format: text

text:
# Properties of Gutzwiller wave-functions for multi-band models

## I Introduction

After so many years since its proposalMott , the Mott transition maintains intact its scientific interest, continually kept alive by the growing number of materials which, on the track of becoming Mott insulators, show rich and fascinating physical properties, among which high Tc superconductivity is likely the most spectacular example.

Much of the theoretical effort has been till now devoted towards cuprate-inspired mainly single-band models, like the standard Hubbard model or its strongly correlated counterpart, the $`tJ`$ model. Although unsolved questions still remain, a lot has been unveiled about the Mott transition in the single-band Hubbard model, especially after the development of the so-called Dynamical Mean Field Theory (DMFT)DMFT . DMFT, which is exact in the limit of infinite coordination lattices, represents a quite reliable approach to investigate the on-site dynamical behavior across the Mott transition. Yet, in its original version, DMFT is not able to yield accurate results concerning inter-site correlations, which are treated in a mean-field like fashion, although several extensions has been proposedClusterDMFT . For that reason, other approaches has been often adopted to investigate the possible occurrence of $`d`$-wave superconductivity in proximity to the Mott insulating phase of single-band Hubbard and $`tJ`$ models, most of which are more or less explicitly related to the Gutzwiller variational techniqueGutzwiller . That amounts to study a variational wave-function consisting of a simple Slater determinant where doubly occupied sites are partially or completely projected out. The Gutzwiller wave-function can be rigorously handled only numerically in finite dimensions, although an approximate analytical scheme to evaluate average values was proposed by Gutzwiller himselfGutzwiller , thereafter called Gutzwiller approximation formula (GAF). Later on, it was realized that the GAF is exact in the case of infinite coordination latticesMetzner  and that it provides similar results to the slave-boson method within the saddle-point approximationMetzner ; Gebhard1 . Moreover, the comparison with the numerical treatment of the Gutzwiller wave-function has shown that the GAF is quite accurate also in finite dimensionsGros ; Shiba ; Gebhard1 .

In spite of the interesting developments achieved by means of the Gutzwiller variational wave-function in single-band Hamiltonians, especially for what it concerns $`d`$-wave superconductivity in the $`t`$-$`J`$ modelGrosSC ; Randeria ; SorellaPRL , there have been not many attempts to apply the same technique to multi-band models, even though there are many interesting strongly correlated materials where orbital degrees of freedom play an important role. Indeed there exist several analyzes based on multi-band extensions of the GAFLu ; Okabe ; Gebhard  and of slave-boson techniqueHasegawa , but a direct comparison with more exact results to test the quality of the Gutzwiller wave-function is still missing, even though DMFT results for multi-band Hubbard models are by now availableRozenberg ; Bulla ; Florens ; Koga ; Massimo .

In this paper, we review the multi-band extension of the Gutzwiller wave-function (GWF) proposed by Ref. Gebhard,  which has the big advantage of being analytically treatable in the limit of infinite dimensions. Moreover, we propose a further extension to broken-symmetry phases, which maintains the same property of being analytically accessible in infinite dimensions without loosing any variational freedom. In this limit, we single out some peculiar properties of the GWF close to the Mott transition. The comparison with exact DMFT results reveals some drawbacks of the GWF which may lead to even qualitatively incorrect results.

The plan of the paper is as follows. In Section II we present the generalized GWF introduced in Ref. Gebhard,  for multi-band models and calculate the general expression of the variational energy in infinite dimensions. In Section III we study the infinite-dimension Mott transition both in the absence and in the presence of multiplet exchange splittings, which we discuss more in detail in Section IV. In Section V we present an extension of the GWF to symmetry broken phases which we use for a specific two-band model in Section VI. Finally, in Section VII we draw some conclusions.

## II The Model

We consider a multi-band lattice model where, on each lattice site, $`N`$ valence orbitals are available for being occupied by the conduction electrons. Since our purpose is to discuss some general features of the Gutzwiller wave-function, we will assume the most simple form of tight-binding energy with nearest neighbor hopping matrix elements diagonal in the orbital index, namely

$$\widehat{H}_0=\frac{t}{\sqrt{z}}\underset{ij}{}\underset{a=1}{\overset{N}{}}\underset{\sigma }{}c_{i,a\sigma }^{}c_{j,a\sigma }^{}+H.c.,$$
(1)
where $`z`$ is the lattice coordination, $`c_{i,a\sigma }^{}`$ and $`c_{i,a\sigma }^{}`$ are, respectively, the annihilation and creation operators for an electron at site $`i`$ with orbital index $`a=1,\mathrm{},N`$ and spin $`\sigma =,`$. The correlation among the electrons is introduced via two local interaction terms: an on-site Hubbard repulsion

$$\widehat{H}_U=\frac{U}{2}\underset{i}{}n_i^2,$$
(2)
where $`n_i=_{a=1}^N_\sigma c_{i,a\sigma }^{}c_{i,a\sigma }^{}`$ is the electron occupation number at site $`i`$, as well as an exchange splitting term

$$\widehat{H}_J=\underset{i}{}\underset{n=0}{\overset{2N}{}}\underset{\mathrm{\Gamma }_n}{}J_{\mathrm{\Gamma }_n}|i,n,\mathrm{\Gamma }_ni,n,\mathrm{\Gamma }_n|.$$
(3)
Here $`|i,n,\mathrm{\Gamma }_n`$ denotes a multiplet of $`n`$-electron states at site $`i`$. $`\mathrm{\Gamma }_n`$ has multiplicity $`g_{\mathrm{\Gamma }_n}`$ such that

$$\underset{\mathrm{\Gamma }_n}{}g_{\mathrm{\Gamma }_n}g_n=\left(\begin{array}{c}2N\\ n\end{array}\right),$$

being the binomial on the right hand side the total number $`g_n`$ of available on-site $`n`$-electron states. Without loss of generality, we assume that $`\widehat{H}_J`$ only splits multiplets at fixed $`n`$ without affecting their center-of-gravity energy, implying $`_{\mathrm{\Gamma }_n}g_{\mathrm{\Gamma }_n}J_{\mathrm{\Gamma }_n}=0`$.

In the absence of interaction, the ground state $`|\mathrm{\Phi }_0`$ of (1) with an average electron number per site $`n_0`$ consists of $`N`$ degenerate bands each one at filling $`n_0/2N`$, namely

$$|\mathrm{\Phi }_0=\underset{a=1}{\overset{N}{}}\underset{|\stackrel{}{k}|<k_F}{}c_{k,a}^{}c_{k,a}^{}|0,$$
(4)
where

$$\frac{2}{V}\underset{|\stackrel{}{k}|<k_F}{}=\frac{n_0}{N}.$$

If $`0<n_0<2N`$ is integer, we expect that, for sufficiently large $`U`$, the correlated ground state of $`\widehat{H}=\widehat{H}_0+\widehat{H}_U+\widehat{H}_J`$ describes a Mott insulator. Analogously to what has been done for the single-band Hubbard model, we would like to have a system which allows to study a metal-to-insulator transition into an ideal paramagnetic Mott insulator, in contrast to a more conventional metal-to-band-insulator transition. For that reason, we further assume that the multiplet exchange splitting (3) is such that, within a perturbation expansion upon the uncorrelated ground state $`|\mathrm{\Phi }_0`$, all self-energy diagrams are diagonal in spin and orbital indices and independent of them. This does not exclude that spin and/or orbital symmetry may be spontaneously broken especially close to the Mott transition, as it is known to occur for the single band Hubbard model on a bipartite lattice. Yet, in what follows, we will discard such a possibility and just discuss an idealized MIT, where neither the metal nor the Mott insulator break any of the symmetries of the Hamiltonian. We postpone to Section V the analysis of spontaneous symmetry breaking nearby the Mott transition.

Having this in mind, we start analyzing the role of the interaction by a Gutzwiller variational approach. Namely we search for the best variational wave-function of the form

$$|\mathrm{\Psi }_G=\underset{i}{}\widehat{P}_{i,G}|\mathrm{\Phi }_0,$$
(5)
where the uncorrelated wave-function is the Slater determinant of Eq. (4) and the on-site Gutzwiller projector Gebhard ; Okabe ,

$$\widehat{P}_{i,G}=\underset{n}{}\underset{\mathrm{\Gamma }_n}{}\lambda _{n\mathrm{\Gamma }_n}|i,n,\mathrm{\Gamma }_ni,n,\mathrm{\Gamma }_n|,$$
(6)
tends to go along with the local interaction terms $`\widehat{H}_U`$ and $`\widehat{H}_J`$ in modifying the relative weights of on-site electronic configurations. As shown by Refs. Gebhard1, ; Gebhard, , there is a certain arbitrariness in the choice of the variational parameters $`\lambda `$’s, related to the fact that any transformation acting on $`|\mathrm{\Phi }_0`$ and involving operators of which $`|\mathrm{\Phi }_0`$ is an eigenstate amounts simply to a multiplicative factor. In our case, that arbitrariness allows to impose, without losing generality, the normalization condition

$$\mathrm{\Phi }_0|\widehat{P}_{i,G}^2|\mathrm{\Phi }_0,$$
(7)
as well as an additional constraint on the single-particle density matrix

$$\mathrm{\Phi }_0|\widehat{P}_{i,G}c_{i,a\sigma }^{}c_{i,b\sigma ^{}}^{}\widehat{P}_{i,G}|\mathrm{\Phi }_0=\mathrm{\Phi }_0|c_{i,a\sigma }^{}c_{i,b\sigma ^{}}^{}|\mathrm{\Phi }_0=\frac{n_0}{2N}\delta _{ab}\delta _{\sigma \sigma ^{}},$$
(8)
where the last equality derives from our choice of $`|\mathrm{\Phi }_0`$, Eq. (4).

A formal solution of (7) and (8) is obtainedGebhard  by writing

$$\lambda _{n\mathrm{\Gamma }_n}^2=\frac{P(n,\mathrm{\Gamma }_n)}{P^{(0)}(n,\mathrm{\Gamma }_n)},$$
(9)
where $`P^{(0)}(n,\mathrm{\Gamma }_n)`$ and $`P(n,\mathrm{\Gamma }_n)`$ represent the occupation probabilities of the $`n`$-electron multiplet $`\mathrm{\Gamma }_n`$ in the uncorrelated, $`|\mathrm{\Phi }_0`$, and correlated, $`|\mathrm{\Psi }_G`$, wave-functions. Upon inserting (9) into (7) and (8) one obtains

$$\underset{n}{}\underset{\mathrm{\Gamma }_n}{}P(n,\mathrm{\Gamma }_n)=1,$$
(10)
the correct normalization for $`P(n,\mathrm{\Gamma }_n)`$, as well as

$$\underset{n}{}\underset{\mathrm{\Gamma }_n}{}nP(n,\mathrm{\Gamma }_n)=n_0,$$
(11)
namely the condition that the average number of electrons per site coincides with the uncorrelated value $`n_0`$. In this representation, the correlated probability distribution is the variational quantity which has to be optimized.

The Eqs. (7) and (8) imply that, within a perturbation expansion in the parameters $`\left(1\lambda _{n\mathrm{\Gamma }_n}\right)`$’s, only more than two fermionic lines can exit from any vertex $`\widehat{P}_{i,G}^2`$ at site $`i`$, see Ref. Gebhard, . This property simplifies considerably the calculations in the limit of infinite coordination lattices, $`z\mathrm{}`$ in Eq. (1), where one can show that the average variational energy per site Gebhard

$`E_G`$ $`=`$ $`{\displaystyle \frac{1}{V}}\mathrm{\Psi }_G|\widehat{H}_0+\widehat{H}_U+\widehat{H}_J|\mathrm{\Psi }_G`$ (12)
$`=`$ $`ZT_0+{\displaystyle \frac{U}{2}}{\displaystyle \underset{n}{}}{\displaystyle \underset{\mathrm{\Gamma }_n}{}}(nn_0)^2P(n,\mathrm{\Gamma }_n)+{\displaystyle \underset{n}{}}{\displaystyle \underset{\mathrm{\Gamma }_n}{}}J_{\mathrm{\Gamma }_n}P(n,\mathrm{\Gamma }_n).`$
$`T_0`$ is the average value per site of the non-interacting tight-binding Hamiltonian

$$T_0=\underset{z\mathrm{}}{lim}\frac{1}{V}\mathrm{\Phi }_0|\widehat{H}_0|\mathrm{\Phi }_0,$$

which gets reduced by a factor $`Z`$ through the Gutzwiller projection, and we have introduced a chemical potential term which makes the Hubbard interaction minimum at the average $`n=n_0`$. Our choice of $`\widehat{H}_J`$ leads to the following expression

$`\sqrt{Z}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{2N}{}}}{\displaystyle \underset{\mathrm{\Gamma }_n,\mathrm{\Gamma }_{n1}}{}}{\displaystyle \frac{n}{2N}}{\displaystyle \frac{g_{\mathrm{\Gamma }_n}g_{\mathrm{\Gamma }_{n1}}}{g_{n1}}}\left({\displaystyle \frac{n_0}{2N}}\right)^{n1}\left(1{\displaystyle \frac{n_0}{2N}}\right)^{2Nn}\lambda _{n\mathrm{\Gamma }_n}\lambda _{n1\mathrm{\Gamma }_{n1}}`$ (13)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{2N}{}}}{\displaystyle \underset{\mathrm{\Gamma }_n,\mathrm{\Gamma }_{n1}}{}}{\displaystyle \frac{n}{g_{n1}}}\sqrt{{\displaystyle \frac{g_{\mathrm{\Gamma }_n}g_{\mathrm{\Gamma }_{n1}}}{n_0\left(2Nn_0\right)}}}\sqrt{P(n,\mathrm{\Gamma }_n)P(n1,\mathrm{\Gamma }_{n1})}.`$
When $`\widehat{H}_J=0`$, namely if nothing in the Hamiltonian splits the degeneracy among states at fixed $`n`$, there should exist a variational solution in which all configurations at fixed $`n`$ are equally probable, a property owned by the Slater determinant $`|\mathrm{\Phi }_0`$. In that case, a multiplet with degeneracy $`g_{\mathrm{\Gamma }_n}`$ occurs with probability

$$P(n,\mathrm{\Gamma }_n)=\frac{g_{\mathrm{\Gamma }_n}}{g_n}P(n),$$
(14)
being $`P(n)`$ the probability of a site to be occupied by $`n`$ electrons irrespective of the configuration. Eq. (14) inserted into (13) leads to the following simple expression of $`Z`$ valid for $`\widehat{H}_J=0`$:

$$\sqrt{Z}=\underset{n=1}{\overset{2N}{}}\frac{n}{\sqrt{n_0\left(2Nn_0\right)}}\sqrt{\frac{g_n}{g_{n1}}}\sqrt{P(n)P(n1)}.$$
(15)
## III Mott transition within the Gutzwiller approach

The search for the optimal $`P(n,\mathrm{\Gamma })`$’s with arbitrary values of the interaction parameters is not straightforward. However it is still possible to obtain simple analytical results close to the Mott transition, which is signaled by a vanishing quasiparticle residue $`Z`$. Indeed, when $`Z1`$, one realizes Lu ; nota  by inspection of Eq. (12) that

$$P(n,\mathrm{\Gamma }_n)Z^{|nn_0|}.$$
(16)
It is therefore justified to assume $`P(n,\mathrm{\Gamma }_n)=0`$ for all $`n`$’s but $`n=n_0`$ and $`n=n_0\pm 1`$, which simplifies remarkably all calculations.

In the simple case $`\widehat{H}_J=0`$, which we denote hereafter as $`J=0`$, this approximation leads through Eq. (15) to

$`\sqrt{Z}`$ $``$ $`{\displaystyle \frac{n_0}{\sqrt{n_0\left(2Nn_0\right)}}}\sqrt{{\displaystyle \frac{g_{n_0}}{g_{n_01}}}}\sqrt{P(n_0)P(n_01)}`$ (17)
$`+{\displaystyle \frac{n_0+1}{\sqrt{n_0\left(2Nn_0\right)}}}\sqrt{{\displaystyle \frac{g_{n_0+1}}{g_{n_0}}}}\sqrt{P(n_0+1)P(n_0)}.`$
The distribution probabilities have to satisfy

$$P(n_01)=P(n_0+1)d,$$

and consequently

$$P(n_0)=12d.$$

By inserting those expressions into (17), one obtains

$`Z`$ $`=`$ $`{\displaystyle \frac{d(12d)}{n_0(2Nn_0)}}\left[\sqrt{n_0(2Nn_0+1)}+\sqrt{(n_0+1)(2Nn_0)}\right]^2`$ (18)
$`\gamma (N,n_0)d(12d),`$
hence

$$E_G=\gamma (N,n_0)d(12d)T_0+Ud.$$
(19)
The optimal $`d`$ which minimizes (19) is readily found:

$$d=\frac{T_0\gamma (N,n_0)+U}{4T_0\gamma (N,n_0)}\frac{U_c(J=0)U}{4U_c(J=0)},$$
(20)
where

$$U_c(J=0)=\gamma (N,n_0)T_0=\frac{1}{2Nn_0}\left[\sqrt{n_0(2Nn_0+1)}+\sqrt{(n_0+1)(2Nn_0)}\right]^2\frac{T_0}{n_0},$$
(21)
is the value of the interaction at the Mott transition Lu , when the optimal $`d=0`$. The variational energy is therefore

$$E_G(J=0)=2T_0\gamma (N,n_0)d^2=\frac{1}{8}\frac{\left[U_c(J=0)U\right]^2}{U_c(J=0)}.$$
(22)
Let us now consider the case $`H_J0`$. We assume that the exchange splitting favors at any given $`n`$ a particular multiplet of states, which we denote as $`\mathrm{\Gamma }_n^{}`$. The Mott insulator described by the Gutzwiller wave-function is therefore characterized by $`P(n_0,\mathrm{\Gamma }_{n_0}^{})=1`$ and has energy $`E_{ins}=J_{\mathrm{\Gamma }_{n_0}^{}}<0`$. In order to better understand how the Mott transition occurs when $`J0`$, it is convenient to consider separately two extreme cases which do not require any numerical calculation.

Let us start by considering an hypothetical metallic solution able to smoothly transform into the Mott insulator, in which therefore only the multiplets $`\mathrm{\Gamma }_n^{}`$ favored by $`H_J`$ are occupied close to the MIT. Since only $`n_0\pm 1`$ and $`n_0`$ are relevant, one obtains an expression of $`Z`$ similar to (18) with the only difference that

$`\gamma (N,n_0)\gamma _{eff}(N,n_0){\displaystyle \frac{1}{n_0(2Nn_0)}}[\sqrt{n_0(2Nn_0+1)}\sqrt{{\displaystyle \frac{g_{\mathrm{\Gamma }_{n_01}^{}}}{g_{n_01}}}{\displaystyle \frac{g_{\mathrm{\Gamma }_{n_0}^{}}}{g_{n_0}}}}`$
$`+\sqrt{(n_0+1)(2Nn_0)}\sqrt{{\displaystyle \frac{g_{\mathrm{\Gamma }_{n_0+1}^{}}}{g_{n_0+1}}}{\displaystyle \frac{g_{\mathrm{\Gamma }_{n_0}^{}}}{g_{n_0}}}}]^2<\gamma (N,n_0).`$
The variational energy reads now

$`E_G(J0)`$ $`=`$ $`\gamma _{eff}(N,n_0)d(12d)T_0+Ud`$
$`+\left(J_{\mathrm{\Gamma }_{n_01}^{}}+J_{\mathrm{\Gamma }_{n_0+1}^{}}2J_{\mathrm{\Gamma }_{n_0}^{}}\right)d+J_{\mathrm{\Gamma }_{n_01}^{}}.`$
Provided we substitute $`\gamma \gamma _{eff}`$ and

$$UU_{eff}U+\left(J_{\mathrm{\Gamma }_{n_01}^{}}+J_{\mathrm{\Gamma }_{n_0+1}^{}}2J_{\mathrm{\Gamma }_{n_0}^{}}\right),$$

which is the actual Hubbard repulsion measured with respect to the energies of the lowest multiplets and not from the centers of gravity, the formal solution has the same expression as before. In particular the critical interaction at which that metallic phase becomes unstable is now

$`U_c(J0)`$ $`=`$ $`\gamma _{eff}(N,n_0)T_0\left(J_{\mathrm{\Gamma }_{n_01}^{}}+J_{\mathrm{\Gamma }_{n_0+1}^{}}2J_{\mathrm{\Gamma }_{n_0}^{}}\right)`$ (23)
$`=`$ $`{\displaystyle \frac{\gamma _{eff}(N,n_0)}{\gamma (N,n_0)}}U_c(J=0)\left(J_{\mathrm{\Gamma }_{n_01}^{}}+J_{\mathrm{\Gamma }_{n_0+1}^{}}2J_{\mathrm{\Gamma }_{n_0}^{}}\right).`$
For very small $`J`$’s, $`U_c(J0)`$ is shifted down with respect to $`U_c(J=0)`$ by terms of order $`|T_0|`$, which already suggests that the above solution is not the most energetically favorable, although it has the merit to merge smoothly into the insulator.

Indeed, the best variational solution is actually different from the above one. If the $`J`$’s are very small compared with $`T_0`$, we rather expect in the metallic phase that the probability distributions of the multiplets at any $`n`$ are only slightly modified with respect to the $`J=0`$ case. If this were true, we could still search for a variational solution of the same form as for $`J=0`$, which has an energy given by (20) with $`d`$ as in Eq.(22). This solution, however, does not converge into the insulating one, which has an energy $`E_{ins}=J_{\mathrm{\Gamma }_{n_0}^{}}`$. The two energies indeed cross when

$$\frac{1}{8}\frac{\left[U_c(J=0)U\right]^2}{U_c(J=0)}=J_{\mathrm{\Gamma }_{n_0}^{}},$$

namely when $`U=U_{}`$,

$$U_{}=U_c(J=0)\sqrt{8U_c(J=0)J_{\mathrm{\Gamma }_{n_0}^{}}}.$$
(24)
For small $`J`$, $`U_{}`$ is larger than $`U_c(J0)`$ given in Eq. (23), which suggests not only that the $`J=0`$ metallic solution has lower energy but also that the Mott transition is first order. The explicit solution of the variational equations shows that, in the optimal metal at small $`J`$, the probability distributions are indeed only slightly modified with respect to the $`J=0`$ case by terms of order $`\sqrt{|J/T_0|}`$ hence that the Mott transition becomes first order as soon as a finite $`J`$ is introduced.

In conclusion, if the multiplet exchange splitting term $`|J|`$ is much smaller that the uncorrelated bandwidth $`W`$, then the Gutzwiller variational approach leads to a first order phase transition from a metal, slightly modified with respect to the $`J=0`$ case, into a Mott insulator instead dominated by $`J`$. This transition is predicted to occur when the quasiparticle residue

$$Z\sqrt{\gamma (N,n_0)\left|\frac{J}{T_0}\right|},$$
(25)
and arises because the metallic solution has to pay too much hopping-energy to modify the relative weights of the multiplets, a cost which overcomes the exchange splitting energy gain.

A first order phase transition has indeed been found in Ref. Gebhard,  by explicit numerical minimization of the variational energy in a two-band model, and agrees with linearized-DMFT results obtained in Ref. Bulla,  on a similar model. Moreover, in the same Ref. Bulla,  it is shown that the first order character reinforces with increasing exchange splitting strength $`J`$ from $`J=0`$, while it weakens with further increasing $`J`$ above some intermediate value. This also agrees with our above results. Indeed, when $`J`$ gets so large to make $`U_c(J0)`$, Eq. (23), greater than $`U_{}`$, Eq. (24), we expect the Mott transition to turn again into a second order one. That anyway requires a substantial $`|J||T_0|`$.

## IV Drawbacks of the Gutzwiller wave-function

In the previous Section we have identified the qualitative behavior across the Mott transition for multi-band Hubbard models within the Gutzwiller variational approach. In particular our analysis has been performed in the limit of infinite coordination lattices, where the Gutzwiller wave-function is able to describe a Mott transition. In the same limit, exact results can be obtained by DMFTDMFT , which allows a direct comparison hence a test on the quality of the Gutzwiller wave-function.

As we previously said, upon approaching the Mott transition with a Gutzwiller wave-function for multi-band models, the on-site charge configurations with a number of electrons $`n`$ different from the average one $`n_0`$ get suppressed approximately like $`Z^{|nn_0|}`$, where $`Z`$ is the quasiparticle residue. It is usually assumed that the GWF in the metallic phase of a single band model provides a faithful description of what Landau quasiparticles do, while it gives a very poor characterization of the underlying insulating component. If this property holds in multi-band models too, then the quasiparticle gas should be characterized by a $`P(n)Z^{|nn_0|}`$. This behavior is not displayed by exact DMFT results, which rather indicate a quasiparticle probability distribution $`P(n)Z`$ (see Ref. Massimo, ). The origin of this disagreement can be easily traced back.

The main reason why the Gutzwiller wave-function gives a poor description of the Mott insulator in the single band Hubbard model at half-filling is that it can not account for virtual excitations into unfavorable charge configurations. Those virtual processes are responsible of the finite average value of the double occupancy $`n_{}n_{}`$ in the insulating phase even in infinite dimensions. Indeed $`n_{}n_{}`$ as function of $`U`$ displays a discontinuity in the slope across the Mott transition in infinite dimensionsDMFT . The singular part, which vanishes at the MIT, is attributed to the quasiparticles and it is qualitatively reproduced by the behavior of the double occupancy as obtained through the Gutzwiller wave-function, supporting the common belief that this wave-function does correctly capture quasiparticle properties.

In a multi-band Hubbard model the situation is different. If we just consider the occupation probability $`P(n_0\pm 1)`$ of $`(n_0\pm 1)`$-charge configurations, we expect across the MIT a behavior similar to the double occupancy in the single-band Hubbard model, still compatible with the GWF. Let us instead consider the occupation probability $`P(n_0\pm 2)`$. In the GWF close to the Mott transition, those configurations get suppressed like $`Z^2`$. In reality, virtual processes from the more advantageous $`n_0`$ and $`(n_0\pm 1)`$ charge configurations imply first of all that $`P(n_0\pm 2)`$ is finite in the insulating phase too, and secondly that the singular $`(n_0\pm 2)`$ quasiparticle contribution still linearly vanishes across the MIT, unlike what it is found by the GWF. That disagreement is more profound than what it would seem to be. A quasiparticle probability distribution of the Gutzwiller type, namely $`P(n)Z^{|nn_0|}`$, suggests that quasiparticles remain more strongly interacting than implied by the true behavior $`P(n)Z`$, even after the Hubbard side-bands are well formed. This indicates that, unlike what happens in the single-band Hubbard model, the multi-band GWF is not fully adequate to capture quasiparticle properties.

The second failure of the Gutzwiller variational approach regards the onset of the first order phase transition, which is predicted to occur when $`Z\sqrt{|J|}`$, for $`|J||T_0|`$. It also originates by the lack of virtual processes. If $`J=0`$, the Gutzwiller wave-function leads to a Mott insulator with a finite entropy, related to the finite number of degenerate on-site electronic configurations with $`n_0`$ electrons. This state has an infinite susceptibility to a term $`\widehat{H}_J`$ which splits that degeneracy, with an energy gain linear in $`J`$. This result is obviously wrong. The super-exchange terms generated by virtual processes into unfavorable charge configurations lead to finite susceptibilities even in the Mott insulator. It implies that the actual energy gain is quadratic in $`J`$ so that the Mott transition is either second order or weakly first order, in this case occurring when $`Z|J|`$.

Indeed, this aspect is not peculiar to a multi-band model but also occurs in a single band model in the presence of a magnetic field $`B`$ which splits spin-up singly occupied sites from spin-down ones. Also in that case the Gutzwiller approach would predict a first order transition when $`Z\sqrt{B}`$, while in reality, being the magnetic susceptibility finiteDMFT , the transition occurs at smaller $`Z`$, see Ref. Georges, , likely when $`ZB`$.

Yet those defects of the Gutzwiller wave-function might not affect qualitatively the physical behavior in the most common situations where the multiplet exchange splitting term leads to the conventional Hund’s rules, namely favors high spin and angular momentum configurations. These states are usually representable by a Slater determinant hence are accessible by a mean-field approach which can be improved by the Gutzwiller projector. However there may be less conventional but still interesting cases where the multiplet exchange splitting favors low degeneracy states, which are not Slater determinants hence unaccessible to mean field theories. Here the Gutzwiller wave-function might be inadequate to describe the Mott transition mainly because it is unable to access the interesting region where the metallic kinetic energy gain $`Z|T_0|`$ competes with the exchange splitting $`J`$.

The inability of the GWF to describe phases which are not accessible by Hartree-Fock, like an ideal Mott insulator, is a well known fact in finite dimensionsShiba . However, it is common opinion that the GWF is instead able to yield a correct description also for Hartree-Fock-unlike physical situations in infinite dimensions. Even more, it is believed that the infinite-dimensional picture provided by the GWF is qualitatively right even in finite dimensions. The above discussion suggests that this common belief is partly wrong and that an improvement of the GWF towards including virtual inter-site processes is necessarynota-b .

## V Gutzwiller wave-functions for symmetry broken phases

In spite of its appealing features, an ideal Mott insulator at zero temperature is unlikely to exist, especially if it has huge degeneracy. Commonly one expects a symmetry broken phase to occur at low temperature. For instance, in a single-band model at half-filling, the ideal Mott insulator has an infinite spin degeneracy which is likely reduced at low temperature by some magnetic ordering. Therefore, even though any mean-field type of approach (including more sophisticated ones based on Density Functional Theory) can only stabilize correlated insulators in broken-symmetry phases, namely can only describe band-insulators thus hiding the basic phenomena leading to a Mott insulator, yet they often provide a faithful description of the low temperature physics.

In this situation, the Gutzwiller variational approach should still be useful to improve Hartree-Fock approximation. That would amount to search for the best wave-function of the form

$$|\mathrm{\Psi }_G(\mathrm{\Delta })=\widehat{P}_G|\mathrm{\Phi }(\mathrm{\Delta }_0)=\underset{i}{}\widehat{P}_{i,G}|\mathrm{\Phi }(\mathrm{\Delta }_0),$$
(26)
with $`\widehat{P}_{i,G}`$ still given by Eq. (6), and where $`|\mathrm{\Phi }(\mathrm{\Delta }_0)`$ is a symmetry-broken uncorrelated trial wave-function with single-particle order parameter $`\mathrm{\Delta }_0`$. In general, after Gutzwiller projection, the correlated wave-function will have a different order parameter $`\mathrm{\Delta }`$. This implies that the average values of the single-particle density matrix over $`|\mathrm{\Psi }_G(\mathrm{\Delta })`$ and $`|\mathrm{\Phi }(\mathrm{\Delta }_0)`$ do not coincide. This does not obviously represent a problem for a numerical treatment, whereas would seem to prevent the use of the method developed in Ref. Gebhard,  for analytically evaluating the average values in infinite dimensions. In fact that method relies on the possibility of constructing a GWF with the same average value of the single particle density matrix as the uncorrelated wave-function. We could still impose that condition for the wave-function defined by (26), but that would reduce the variational freedom.

In this Section we present a simple extension of the GWF to account for symmetry broken phases while leaving the property of being analytically treatable in infinite dimensions without any variational loss. We start by noticing that there always exists a non unitary operator $`\widehat{U}=_i\widehat{U}_i`$ such that its action

$$\widehat{U}|\mathrm{\Phi }(\mathrm{\Delta }_0)=|\mathrm{\Phi }(\mathrm{\Delta }),$$
(27)
leads to a trial wave-function $`|\mathrm{\Phi }(\mathrm{\Delta })`$ of the same form as $`|\mathrm{\Phi }(\mathrm{\Delta }_0)`$ but with the same average value of the single-particle density matrix as the Gutzwiller projected $`|\mathrm{\Psi }_G(\mathrm{\Delta })`$, hence the same order parameter $`\mathrm{\Delta }`$. Therefore,

$$|\mathrm{\Psi }_G(\mathrm{\Delta })=\widehat{P}_G|\mathrm{\Phi }(\mathrm{\Delta }_0)=\widehat{P}_G\widehat{U}|\mathrm{\Phi }(\mathrm{\Delta })\widehat{\mathrm{P}}_G|\mathrm{\Phi }(\mathrm{\Delta }),$$
(28)
which implies that, provided we substitute

$$\widehat{P}_G\widehat{\mathrm{P}}_G=\widehat{P}_G\widehat{U},$$
(29)
we can still search without loss of generality for a variational wave-function

$$|\mathrm{\Psi }_G(\mathrm{\Delta })=\widehat{\mathrm{P}}_G|\mathrm{\Phi }(\mathrm{\Delta }),$$

where the average single-particle density matrix stays unchanged after Gutzwiller projection. The cost is that we must work with a Gutzwiller projector as given by (29), which is in general neither diagonal in the multiplets which appears in $`\widehat{H}_J`$, Eq. (3), nor hermitean.

Actually we can identify two distinct situations which may occur, one of them being already included in the formalism developed by Ref. Gebhard, . If the exchange splitting $`J`$ favors a degenerate atomic configuration, then we reasonably expect that in the true ground state of the lattice only one of the degenerate states will be occupied on a given site, eventually changing from site to site. There the order parameter corresponds locally to a conserved quantity of the atomic Hamiltonian, for instance the $`z`$-component of the on-site spin, which leads to a generalized Gutzwiller projector $`\widehat{\mathrm{P}}_G`$ still hermitean and diagonal. Let us consider as a simple example a one band model. The local Gutzwiller projector is in general

$$\widehat{P}_{i,G}=\lambda _0|i,0i,0|+\lambda _2|i,2i,2|+\lambda _1[|i,i,|+|i,i,|],$$

where $`|i,0(2)`$ denote the empty or doubly-occupied site $`i`$, while $`|i,\sigma `$ the singly occupied site with spin $`\sigma `$. We assume that the uncorrelated wave-function is magnetically ordered. Then (27) may be constructed by local operators

$$\widehat{U}_i=\mathrm{e}^{2\alpha _i\widehat{S}_i^z},$$

where $`\widehat{S}_i^z`$ is the $`z`$-component of the spin operator at site $`i`$. We find that

$$\widehat{\mathrm{P}}_{i,G}=\widehat{P}_{i,G}\widehat{U}_i=\lambda _0|i,0i,0|+\lambda _2|i,2i,2|+\lambda _1\mathrm{e}^{\alpha _i}|i,i,|+\lambda _1\mathrm{e}^{\alpha _i}|i,i,|.$$

Indeed the modified Gutzwiller projector is still diagonal and hermitean, although there appear different variational parameters for spin up and down components. This additional degree of freedom gets fixed once we require that the order parameters of the correlated and uncorrelated wave-functions coincide. The above type of GWF is actually a particular case of the generalized GWF’s introduced in Ref. Gebhard, .

We can however envisage a different situation where the uncorrelated wave-function $`|\mathrm{\Phi }_0`$ has an order parameter which does not correspond locally to a conserved quantity. There a modified Gutzwiller projector $`\widehat{\mathrm{P}}_G`$ is unavoidably off-diagonal and non-hermitean. Let us discuss an over-simplified example. We assume that the exchange term leads to a local Gutzwiller projector of the form (we drop the site label)

$$\widehat{P}_G=\lambda _a|aa|+\lambda _b|bb|,$$

and moreover that the uncorrelated wave-function has an order parameter identified by non-zero matrix elements

$$\mathrm{\Phi }_0|ab|\mathrm{\Phi }_0=\mathrm{\Phi }_0|ba|\mathrm{\Phi }_0.$$

The transformation $`\widehat{U}`$ can be now taken of the form

$$\widehat{U}=\alpha \left(|aa|+|bb|\right)+\beta \left(|ab|+|ba|\right).$$

We readily find that

$$\widehat{\mathrm{P}}_G=\widehat{P}_G\widehat{U}=\alpha \lambda _a|aa|+\alpha \lambda _b|bb|+\beta \lambda _a|ab|+\beta \lambda _b|ba|,$$

indeed containing off-diagonal terms and evidently non-hermitean. This is a novel situation which we are going to discuss more in detail in a particular example.

We conclude this Section by pointing out that the non-hermitean character plays a crucial role only if the ordering involves just the quasiparticles, while it is essentially irrelevant when both the quasiparticles and the Mott-Hubbard side-bands contribute to the order parameter. In the following section we analyse a two-band model where both cases may appear.

## VI A two-band model study

We consider a two-band Hubbard model described by the Hamiltonian (1) and (2) where the orbital index $`a=1,2`$. Besides the local spin-density operators

$$\stackrel{}{S}_i=\frac{1}{2}\underset{a=1}{\overset{2}{}}\underset{\alpha \beta }{}c_{i,a\alpha }^{}\stackrel{}{\sigma }_{\alpha \beta }c_{i,a\beta }^{},$$

where $`\widehat{\sigma }_x`$, $`\widehat{\sigma }_y`$ and $`\widehat{\sigma }_z`$ are the Pauli matrices, we introduce orbital pseudo-spin operators

$$\stackrel{}{T}_i=\frac{1}{2}\underset{a,b=1}{\overset{2}{}}\underset{\sigma }{}c_{i,a\sigma }^{}\stackrel{}{\tau }_{ab}c_{i,b\sigma }^{},$$
(30)
where the Pauli matrices $`\widehat{\tau }`$’s act on the orbital indices. The hopping and Hubbard terms, (1) and (2), have a very large SU(4) symmetry. Having in mind common physical realizations, like e.g. $`d`$-orbitals of $`e_g`$ symmetry, we assume that SU(4) is lowered down to a spin SU(2) times an orbital O(2) by the exchange term, which can therefore be written as

$$\widehat{H}_J=\underset{i}{}J_S\stackrel{}{S}_i\stackrel{}{S}_i+J_T\stackrel{}{T}_i\stackrel{}{T}_i3\left(J_S+J_T\right)\left(T_i^z\right)^2.$$
(31)
The above $`\widehat{H}_J`$ just splits the on-site configurations with two-electrons. There are six of those states on a given site $`i`$, a spin-triplet orbital-singlet $`|i,n=2;S=1,S_z;T=0`$, which we denote hereafter as $`|i,2,t`$, and a spin-singlet orbital-triplet, for which we use the short notations $`|i,n=2;S=0;T=1,T_z=0|i,2,0`$ and $`|i,n=2;S=0;T=1,T_z=\pm 1|i,2,\pm `$. In this subspace, $`\widehat{H}_J`$ has the form

$$\widehat{H}_J=\underset{i}{}\mathrm{\hspace{0.17em}2}J_S|i,2,ti,2,t|+2J_T|i,2,0i,2,0|(3J_S+J_T)|i,2,\pm i,2,\pm |.$$
(32)
The standard Hund’s rules correspond to $`J_T<J_S<5J_T/6<0`$, when the spin-triplet orbital-singlet configuration has the lowest energy, followed by the spin-singlet orbital doublet with $`T_z=\pm 1`$. In this case the ideal Mott insulator at half-filling, $`n_0=2`$, represents localized spin-1 moments which should order at low enough temperature to freeze out the spin-entropy. On general grounds one expects that the magnetic ordering in the insulator should contaminate the nearby metallic phase so that, as $`U`$ increases from weak coupling, first a transition from a paramagnetic into a magnetic metal should occur, followed by a Mott transition into a magnetically ordered insulator. Even a mean field approach is in principle able to reproduce the above scenario. In this situation, as we discussed before, the Gutzwiller-projected Hartree-Fock wave-function does improve the mean-field solution, providing a better physical description. In Appendix I we analyse in detail the case in which a bipartite lattice stabilizes an antiferromagnetic ordering.

Less conventional is the situation $`J_T<|J_S|`$, where the non-degenerate spin-singlet with $`T_z=0`$, namely

$$|i,2;S=0;T=1,T_z=0=\frac{1}{\sqrt{2}}\left(c_{i,1}^{}c_{i,2}^{}+c_{i,2}^{}c_{i,1}^{}\right)|0,$$
(33)
is the lowest energy configuration. That would be for instance the case of two Hubbard models (two-chains, two-planes, etc…), coupled by an antiferromagnetic exchange. Here the large $`U`$ Mott insulator with $`n_0=2`$ describes a collection of on-site singlets, a local version of a valence-bond (VB) insulator. Since it is not degenerate and fully gaped, we expect the VB insulator to be stable at large $`U`$ against any spin and/or orbital order. Just to avoid unessential complications, we assume that the lattice is sufficiently frustrated to prevent any spin/orbital ordering at any $`U`$. This situation is far less trivial than the previous one. In fact, being the Mott insulator not describable by a single Slater determinant, it is inherently unreachable by any mean-field approach, which necessarily leads to some kind of ordered state. According to our previous discussion, we expect in this case that also the GWF is unable to provide a faithful description of the Mott transition.

As we showed in Section III, the GWF without any symmetry breaking would undergo a first order metal-insulator transition when the quasiparticle residue

$$Z\sqrt{\frac{J_T}{6T_0}},$$

namely when $`U=U_{}`$ being

$$U_{}=6T_04\sqrt{6T_0J_T}.$$

We also argued that this result is wrong since a metallic phase is able to enter the regime in which $`Z|J_T/T_0|`$. Let us now check whether there exists a better broken-symmetry metallic solution. Indeed, even if lattice frustration prevents spin/orbital order, there is still a broken-symmetry GWF which might in principle compete with the above metallic solution.

### VI.1 Superconducting Gutzwiller wave-function

When $`U=0`$, the multiplet exchange term (32) favors a BCS Hartree-Fock wave-function with the $`s`$-wave order parameter

$$\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0)|c_{i,1}^{}c_{i,2}^{}|\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0)=\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0)|c_{i,2}^{}c_{i,1}^{}|\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0)\mathrm{\Delta }_0\frac{1}{2}.$$
(34)
When $`\mathrm{\Delta }_01/2`$, the doubly occupied sites in the spin-triplet, $`|2,t`$, or in the doublet of spin-singlets, $`|2,\pm `$, configurations get suppressed by a factor $`(12\mathrm{\Delta }_0)^2`$ with respect to the $`T_z=0`$ spin-singlet, $`|2,0`$. Similarly, the probability of singly, $`|1`$, or triply-occupied, $`|3`$, sites vanish like $`(12\mathrm{\Delta }_0)`$. This suggests that by Gutzwiller projecting out sites with zero and four electrons, $`|0`$ and $`|4`$, respectively, through the variational wave-function

$$|\mathrm{\Psi }_G(\mathrm{\Delta })=\widehat{P}_G|\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0),$$
(35)
one might indeed smoothly connect to the VB insulating state, with $`(12\mathrm{\Delta }_0)`$ playing the role of $`Z`$. However, even though the uncorrelated BCS wave-function has a large order parameter $`\mathrm{\Delta }_01/2`$, the correlated $`|\mathrm{\Psi }_G`$ should have a much smaller one, $`\mathrm{\Delta }Z\mathrm{\Delta }_0`$, since only quasiparticles are involved in superconductivity. In such a case we are therefore obliged to implement the non-hermitean Gutzwiller projector in order to get analytical results for large coordination lattices. In other words, we shall work with a Gutzwiller wave-function of the same form as (35) and impose that $`|\mathrm{\Psi }_G`$ and $`|\mathrm{\Phi }_{BCS}`$ have the same order parameter $`\mathrm{\Delta }`$ through a non-hermitean $`\widehat{\mathrm{P}}_G`$, see Eq. (29). Since this is a novel situation in the Gutzwiller variational approach, we prefer to describe it in detail.

In order to simplify the analysis at half-filling, we assume that the Hamiltonian has particle-hole symmetry and search for solutions which do not break it. This guarantees that there are still two conditions we can impose without losing variational freedom: an overall normalization, (7), and a particle-hole symmetry constraint.

To accomplish this job, it is convenient to work not in the original electron basis but in the natural basis where the on-site single-particle density matrix is diagonal. This is done by the following unitary transformation, which is valid at half-filling $`n_0=2`$,

$$\left(\begin{array}{c}c_{i,1(2)}^{}\\ c_{i,2(1)}^{}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1\hfill & \hfill 1\\ 1\hfill & \hfill 1\end{array}\right)\left(\begin{array}{c}a_{i,1(2)}^{}\\ a_{i,2(1)}^{}\end{array}\right).$$

In the natural basis the only non-zero on-site average is

$$\mathrm{\Phi }_{BCS}(\mathrm{\Delta })|a_{i,a\sigma }^{}a_{i,a\sigma }^{}|\mathrm{\Phi }_{BCS}(\mathrm{\Delta })=\frac{1}{2}\mathrm{\Delta }.$$

In order to distinguish the local configurations in the natural from the original basis we will denote the former ones as $`|\overline{n},\mathrm{\Gamma }_{\overline{n}}`$. The most general local Gutzwiller projector is in this case

$$\widehat{\mathrm{P}}_{i,G}=\underset{\overline{n}}{}\underset{\mathrm{\Gamma }_{\overline{n}}}{}\lambda _{\overline{n}\mathrm{\Gamma }_{\overline{n}}}|i,\overline{n},\mathrm{\Gamma }_{\overline{n}}i,\overline{n},\mathrm{\Gamma }_{\overline{n}}|+\lambda _{\overline{0}\overline{4}}|i,\overline{0}i,\overline{4}|+\lambda _{\overline{4}\overline{0}}|i,\overline{4}i,\overline{0}|.$$
(36)
The last two terms are the only possible off-diagonal elements at half-filling when particle-hole symmetry holds. The normalization condition and the conservation of the single-particle density-matrix lead to the following parametrization of the $`\lambda `$’s for $`\overline{n}0,4`$

$$\lambda _{\overline{n}\mathrm{\Gamma }_{\overline{n}}}^2=\frac{P(\overline{n},\mathrm{\Gamma }_{\overline{n}})}{P^{(0)}(\overline{n},\mathrm{\Gamma }_{\overline{n}})},$$

where $`P(\overline{n},\mathrm{\Gamma }_{\overline{n}})`$ and $`P^{(0)}(\overline{n},\mathrm{\Gamma }_{\overline{n}})`$ are the correlated and uncorrelated occupation probabilities in the natural basis. For $`\overline{n}=0,4`$ we have instead

$`\lambda _{\overline{0}}^2P^{(0)}(\overline{0})+\lambda _{\overline{0}\overline{4}}^2P^{(0)}(\overline{4})`$ $`=`$ $`P(\overline{0}),`$ (37)
$`\lambda _{\overline{4}\overline{0}}^2P^{(0)}(\overline{0})+\lambda _{\overline{4}}^2P^{(0)}(\overline{4})`$ $`=`$ $`P(\overline{4}),`$ (38)
$`\lambda _{\overline{0}}\lambda _{\overline{4}\overline{0}}P^{(0)}(\overline{0})+\lambda _{\overline{0}\overline{4}}\lambda _{\overline{4}}P^{(0)}(\overline{4})`$ $`=`$ $`A(\overline{0};\overline{4})=A(\overline{4};\overline{0}),`$ (39)
where we introduce the transition amplitudes

$$A(\overline{n},\mathrm{\Gamma }_{\overline{n}};\overline{m},\mathrm{\Gamma }_{\overline{m}})=\mathrm{\Psi }_G|i,\overline{n},\mathrm{\Gamma }_{\overline{n}}i,\overline{m},\mathrm{\Gamma }_{\overline{m}}|\mathrm{\Psi }_G.$$

The occupation probabilities in the natural basis are related to those in the original one through

$$\begin{array}{ccc}P(\overline{0})\hfill & =& \frac{1}{2}\left(P(0)+P(2,0)+A(0;4)\right)+\sqrt{2}A(0;2,0),\hfill \\ P(\overline{1})\hfill & =& P(1)+A(1;3),\hfill \\ P(\overline{2},0)\hfill & =& P(0)A(0;4),\hfill \\ P(\overline{2},\pm )\hfill & =& P(2,\pm ),\hfill \\ P(\overline{2},t)\hfill & =& P(2,t),\hfill \\ P(\overline{3})\hfill & =& P(1)A(1;3),\hfill \\ P(\overline{4})\hfill & =& \frac{1}{2}\left(P(0)+P(2,0)+A(0;4)\right)\sqrt{2}A(0;2,0),\hfill \\ A(\overline{0};\overline{4})\hfill & =& \frac{1}{2}\left(P(0)P(2,0)+A(0;4)\right),\hfill \end{array}$$
(40)
where we have used the fact that, by particle-hole symmetry, $`P(1)=P(3)`$ and $`P(0)=P(4)`$. The order parameter $`\mathrm{\Delta }`$ is given by

$$2\mathrm{\Delta }=2\sqrt{2}A(02+)+A(13),$$

and the following inequalities should be verified

$`2A(0;2,0)^2`$ $``$ $`P(2,0)\left[P(0)+A(0;4)\right],`$
$`A(0;4)`$ $``$ $`P(0),`$
$`A(1;3)`$ $``$ $`P(1).`$
The big advantage in working with natural orbitals is that the $`Z`$ reduction factor has the same expression as in (13), namely

$$\sqrt{Z(\mathrm{\Delta })}=\underset{\overline{n}=1}{\overset{4}{}}\underset{\mathrm{\Gamma }_{\overline{n}},\mathrm{\Gamma }_{\overline{n}1}}{}\frac{\overline{n}}{4}\frac{g_{\mathrm{\Gamma }_{\overline{n}}}g_{\mathrm{\Gamma }_{\overline{n}1}}}{g_{\overline{n}1}}\left(\frac{1}{2}\mathrm{\Delta }\right)^{\overline{n}1}\left(\frac{1}{2}+\mathrm{\Delta }\right)^{4\overline{n}}\lambda _{\overline{n}\mathrm{\Gamma }_{\overline{n}}}\lambda _{\overline{n}1\mathrm{\Gamma }_{\overline{n}1}}.$$
(41)
Once we know how to relate the parameters $`\lambda `$’s and $`Z`$ to the variational occupation probabilities defining the correlated and uncorrelated wave-functions, we can solve the most general variational problem by minimizing the energy functional

$`E_G(\mathrm{\Delta })=Z(\mathrm{\Delta })T_0(\mathrm{\Delta })+U\left[4P(0)+P(1)\right]`$
$`+2J_SP(2,t)+2J_TP(2,0)(3J_S+J_T)P(2,\pm ),`$ (42)
where

$$\mathrm{\Phi }_{BCS}(\mathrm{\Delta })|\widehat{H}_0|\mathrm{\Phi }_{BCS}(\mathrm{\Delta })=T_0(\mathrm{\Delta }).$$

For the sake of clarity, we present here an analysis based on the following parametrization of the $`\lambda `$’s in Eq. (36), although it has less variational freedom:

$$\begin{array}{ccc}\lambda _{\overline{0}}\hfill & =& \frac{1}{2}\left(\lambda _0+\lambda _{20}\right)\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{0})}{P_\mathrm{\Delta }^{(0)}(\overline{0})}},\hfill \\ \lambda _{\overline{4}}\hfill & =& \frac{1}{2}\left(\lambda _0+\lambda _{20}\right)\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{4})}{P_\mathrm{\Delta }^{(0)}(\overline{4})}},\hfill \\ \lambda _{\overline{0}\overline{4}}\hfill & =& \frac{1}{2}\left(\lambda _0\lambda _{20}\right)\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{4})}{P_\mathrm{\Delta }^{(0)}(\overline{4})}},\hfill \\ \lambda _{\overline{4}\overline{0}}\hfill & =& \frac{1}{2}\left(\lambda _0\lambda _{20}\right)\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{0})}{P_\mathrm{\Delta }^{(0)}(\overline{0})}},\hfill \\ \lambda _{\overline{1}}\hfill & =& \lambda _1\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{1})}{P_\mathrm{\Delta }^{(0)}(\overline{1})}},\hfill \\ \lambda _{\overline{3}}\hfill & =& \lambda _1\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{3})}{P_\mathrm{\Delta }^{(0)}(\overline{3})}},\hfill \\ \lambda _{\overline{2}0}\hfill & =& \lambda _0\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{2},0)}{P_\mathrm{\Delta }^{(0)}(\overline{2},0)}},\hfill \\ \lambda _{\overline{2}t}\hfill & =& \lambda _{2t}\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{2},t)}{P_\mathrm{\Delta }^{(0)}(\overline{2},t)}},\hfill \\ \lambda _{\overline{2}\pm }\hfill & =& \lambda _{2\pm }\sqrt{\frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{2},\pm )}{P_\mathrm{\Delta }^{(0)}(\overline{2},\pm )}}.\hfill \end{array}$$
(43)
As before $`\mathrm{\Delta }_0`$ is the order parameter of the uncorrelated BCS wave-function while $`\mathrm{\Delta }`$ is the true order parameter after Gutzwiller projection, see Eq. (35). $`P_{\mathrm{\Delta }_0}^{(0)}(\overline{n},\mathrm{\Gamma }_{\overline{n}})`$ and $`P_\mathrm{\Delta }^{(0)}(\overline{n},\mathrm{\Gamma }_{\overline{n}})`$ are the distribution probabilities in the natural basis for the BCS wave-functions with order parameter $`\mathrm{\Delta }_0`$ and $`\mathrm{\Delta }`$, respectively. They are explicitly written in the Appendix, Eq. (55).

The normalization condition as well as the conservation of the single particle density matrix imply that

$$\lambda _{n\mathrm{\Gamma }_n}^2=\frac{P(n,\mathrm{\Gamma }_n)}{P_{\mathrm{\Delta }_0}^{(0)}(n,\mathrm{\Gamma }_n)},$$
(44)
where $`P_{\mathrm{\Delta }_0}^{(0)}(n,\mathrm{\Gamma }_n)=\mathrm{\Phi }_{\mathrm{\Delta }_0}|i;n,\mathrm{\Gamma }_ni;n,\mathrm{\Gamma }_n|\mathrm{\Phi }_{\mathrm{\Delta }_0}`$ is the occupation probability for configurations in the original electronic basis within the uncorrelated BCS wave-function with large order parameter $`\mathrm{\Delta }_0`$, see Eq. (56), while $`P(n,\mathrm{\Gamma }_n)`$ is the same quantity for the correlated wave-function. Eq. (44) is the most natural generalization of (9) to a broken-symmetry phase, which is the reason why we have chosen the above parametrization. The true order parameter $`\mathrm{\Delta }`$ is defined through

$`2\mathrm{\Delta }`$ $`=`$ $`P_\mathrm{\Delta }^{(0)}(\overline{0})+{\displaystyle \frac{1}{2}}P_\mathrm{\Delta }^{(0)}(\overline{1}){\displaystyle \frac{1}{2}}P_\mathrm{\Delta }^{(0)}(\overline{3})+P_\mathrm{\Delta }^{(0)}(\overline{4})`$ (45)
$`=`$ $`\mathrm{\Delta }_0\left(1+4\mathrm{\Delta }_0^2\right)\sqrt{{\displaystyle \frac{P(2,0)P(0)}{P_{\mathrm{\Delta }_0}^{(0)}(2,0)P_{\mathrm{\Delta }_0}^{(0)}(0)}}}+\mathrm{\Delta }_0\left(14\mathrm{\Delta }_0^2\right){\displaystyle \frac{P(1)}{P_{\mathrm{\Delta }_0}^{(0)}(1)}},`$
which indeed is of order $`Z`$ when $`1/2\mathrm{\Delta }_0Z`$, $`P(2,0)1`$, $`P(1)Z`$ and $`P(0)Z^2`$. The explicit evaluation of the $`Z`$ reduction factor is presented in the Appendix, see Eq. (53).

Let us now compare the variational energy as given by Eq. (42) with $`Z`$ of Eq. (57), valid for $`1/2\mathrm{\Delta }_0=\delta 1`$, to the energy of a non superconducting paramagnetic solution, Eq. (42) with $`\mathrm{\Delta }=0`$, $`Z`$ being given by Eq. (54). We find that the Gutzwiller projected BCS wave-function has always higher energy by terms roughly of order $`Z|T_0|`$.

Therefore, even though the Gutzwiller projector is quite efficient to transform the huge Hartree-Fock energy cost, namely

$`\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0)|\widehat{H}|\mathrm{\Phi }_{BCS}(\mathrm{\Delta }_0)\mathrm{\Phi }_{BCS}(0)|\widehat{H}|\mathrm{\Phi }_{BCS}(0)`$ $`=`$ $`T_0(\mathrm{\Delta }_0)T_0(0)+2U\mathrm{\Delta }_0^2+4J_T\mathrm{\Delta }_0^2`$
$`T_0(0)+{\displaystyle \frac{U}{2}}+J_T,`$
into a much smaller one of order $`Z|T_0|ZU`$ close to the Mott transition, yet it is not able to make superconductivity favorable. Namely, the best variational metallic solution remains the one described in Section III, with all the drawbacks discussed in Section IV. In conclusion, as we anticipated, the Gutzwiller variational approach does not properly describe a mean-field-unlike Mott transition.

In reality we may expect a superconducting phase just before the VB Mott insulator. Recently it has been studied a model which shares many common features with the present one, namely a three-band Hubbard model with inverted Hund’s rules Massimo , mimicking a strongly dynamical Jahn-Teller effect. For an average number of electrons per site $`n_0=2`$, the inverted Hund’s rules favor, like in our example, a non degenerate singlet on-site configuration. By a DMFT calculation, a superconducting instability was discovered just before the singlet Mott insulator. However, that instability was found to appear when the quasiparticle residue $`Z|J|`$ (see Ref. Massimo, ). As we discussed at length previously, the simplest metallic Gutzwiller wave-function which we have so far considered is unable to reach $`ZJ`$, since it becomes disadvantageous with respect to the insulating one already at $`Z\sqrt{J}`$. Therefore we can not exclude that superconductivity may occur even in the two-band model we have considered.

## VII Conclusions

In this paper we have analysed some peculiar features of the Mott transition displayed by a multi-band Gutzwiller variational wave-function (GWF) in infinite dimensions. The analysis has been carried out by using the generalized GWF introduced in Ref. Gebhard, , which allows a simple analytical treatment in infinite dimensions. Moreover, we have extended that wave-function to account for broken-symmetry phases.

It is usually assumed that the GWF in infinite dimensions gives a faithful description of the quasiparticle behavior around the Mott metal-insulator transition. We have shown that while this belief is partly true for single-band models, it is incorrect for multi-band models. In particular, we have identified at least two major failures of the GWF across the Mott transition. The first concerns the occupation probability $`P(n)`$ of on-site charge configurations with $`n`$ electrons different from the integer average one $`n_0`$, which is believed to represent just the quasiparticle occupation probability normalized to the quasiparticle residue $`Z`$. The GWF in infinite dimensions predicts $`P(n)Z^{|nn_0|}`$ close to the Mott transition, while both physical arguments as well as Dynamical Mean Field Theory results suggest a $`P(nn_0)Z`$, even in infinite dimensions. This apparently innocuous disagreement is instead profound. In fact, the GWF results imply that the quasiparticles remain much more strongly interacting than what the correct $`P(n)Z`$ behavior suggests.

Another drawback concerns the Mott transition in the presence of a weak multiplet exchange splitting term $`J`$. Within the GWF, the Mott transition turns into a first order one and occurs when the quasiparticle residue $`Z\sqrt{|J|/W}`$, $`W`$ being the bare bandwidth, much before the quasiparticle gas has had the time to react against $`J`$. This happens because the susceptibility to an infinitesimal exchange splitting $`J`$ diverges at the Mott transition for a GWF. In reality, that susceptibility is finite so that the Mott transition is either second order or weakly first order, in that case occurring when $`Z|J|/W`$. The main consequence is that the interesting region where the metallic hopping-energy gain $`ZW`$ competes against the exchange $`J`$ is not even accessible by a Gutzwiller wave-function. Both the above mentioned shortcomings have the same origin: the inability of the GWF to account for virtual processes into unfavorable charge configurations.

We have then argued, on the basis of a two-band model study, that the GWF is still a good variational wave-function in all cases which can be qualitatively described by a mean-field theory, but it fails otherwise, as for instance in the case we have explicitly analysed where the Mott insulator is a local version of a valence bond insulator. There an improvement of the GWF is necessary.

###### Acknowledgements.

We acknowledge helpful discussions with Joerg Büenemann, Sandro Sorella and Federico Becca. We are particularly grateful to Claudio Castellani for his un-valuable comments and suggestions. This work was partly supported by MIUR COFIN-2001.

## Appendix A Antiferromagnetic Gutzwiller wave-function for a two-band model

The best Hartree-Fock wave-function $`|\mathrm{\Phi }_0(m)`$ on a bipartite lattice when the multiplet exchange term favors the spin-triplet configuration is the ground state of an Hamiltonian

$$\widehat{H}_{HF}=\widehat{H}_0+\mu \underset{i}{}\underset{a=1,2}{}(1)^i\left(n_{i,a}n_{i,a}\right),$$

which describes an antiferromagnetic insulator with order parameter

$$\mathrm{\Phi }_0|n_{i,1}n_{i,1}|\mathrm{\Phi }_0=\mathrm{\Phi }_0|n_{i,2}n_{i,2}|\mathrm{\Phi }_0=2(1)^im.$$

We search for the optimal Gutzwiller wave-function

$$|\mathrm{\Psi }_G=\underset{i}{}\widehat{P}_{i,G}|\mathrm{\Phi }_0(m),$$

where for a given sublattice and making use of particle-hole symmetry,

$`\widehat{P}_{i,G}`$ $`=`$ $`\lambda _0\left[|i,0i,0|+|i,4i,4|\right]`$
$`+\lambda _{1+}\left[|i,1;S_z=1/2i,1;S_z=1/2|+|i,3;S_z=1/2i,3;S_z=1/2|\right]`$
$`+\lambda _1\left[|i,1;S_z=1/2i,1;S_z=1/2|+|i,3;S_z=1/2i,3;S_z=1/2|\right]`$
$`+\lambda _{2\pm }|i,2,\pm i,2,\pm |+\lambda _{20}|i,2,0i,2,0|`$
$`+\lambda _{2t+}|i,2;S=1,S_z=1i,2;S=1,S_z=1|`$
$`+\lambda _{2t}|i,2;S=1,S_z=1i,2;S=1,S_z=1|`$
$`+\lambda _{2t0}|i,2;S=1,S_z=0i,2;S=1,S_z=0|`$
while for the other sublattice $`+`$ with $``$ interchange. In this particular case the Gutzwiller projector remains hermitean since the non-unitary transformation $`\widehat{U}_i`$ is diagonal in the above multiplet basis.

We define the correlated probability distributions for a given sublattice

$`P(1,+)`$ $`=`$ $`P(1,S_z=1/2)=P(3,S_z=1/2),`$
$`P(1,)`$ $`=`$ $`P(1,S_z=1/2)=P(3,S_z=1/2),`$
$`P(2,t+)`$ $`=`$ $`P(2,S=1,S_z=1),`$
$`P(2,t)`$ $`=`$ $`P(2,S_z=1),`$
$`P(2,t0)`$ $`=`$ $`P(2,S_z=0),`$
and analogously for the uncorrelated ones $`P^{(0)}(\mathrm{})`$. For the other sublattice $`S_zS_z`$.

The two conditions (7) and (8) imply that

$`\lambda _0^2`$ $`=`$ $`{\displaystyle \frac{P(0)}{P^{(0)}(0)}}={\displaystyle \frac{P(0)}{\left(\frac{1}{4}m^2\right)^2}},`$
$`\lambda _{1+}^2`$ $`=`$ $`{\displaystyle \frac{P(1,+)}{P^{(0)}(1,+)}}={\displaystyle \frac{P(1,+)}{2\left(\frac{1}{2}+m\right)^3\left(\frac{1}{2}m\right)}},`$
$`\lambda _1^2`$ $`=`$ $`{\displaystyle \frac{P(1,)}{P^{(0)}(1,)}}={\displaystyle \frac{P(1,)}{2\left(\frac{1}{2}m\right)^3\left(\frac{1}{2}+m\right)}},`$
$`\lambda _{2,\pm }^2`$ $`=`$ $`{\displaystyle \frac{P(2,\pm )}{P^{(0)}(2,\pm )}}={\displaystyle \frac{P(2,\pm )}{2\left(\frac{1}{4}m^2\right)^2}},`$
$`\lambda _{2,0}^2`$ $`=`$ $`{\displaystyle \frac{P(2,0)}{P^{(0)}(2,0)}}={\displaystyle \frac{P(2,0)}{\left(\frac{1}{4}m^2\right)^2}},`$
$`\lambda _{2,t0}^2`$ $`=`$ $`{\displaystyle \frac{P(2,t0)}{P^{(0)}(2,t0)}}={\displaystyle \frac{P(2,t0)}{\left(\frac{1}{4}m^2\right)^2}},`$
$`\lambda _{2,t+}^2`$ $`=`$ $`{\displaystyle \frac{P(2,t+)}{P^{(0)}(2,t+)}}={\displaystyle \frac{P(2,t+)}{\left(\frac{1}{2}+m\right)^4}},`$
$`\lambda _{2,t}^2`$ $`=`$ $`{\displaystyle \frac{P(2,t)}{P^{(0)}(2,t)}}={\displaystyle \frac{P(2,t)}{\left(\frac{1}{2}m\right)^4}}.`$
The correlated occupation probabilities satisfy the normalization condition

$$2P(0)+2P(1,+)+2P(1,)+P(2,0)+P(2,\pm )+P(2,t0)+P(2,t+)+P(2,t)=1,$$
(46)
as well as the conservation of the order parameter

$$P(1,+)+P(2,t+)P(1,)P(2,t)=2m.$$
(47)
By all the above defined quantities, the variational energy is found to be

$`E_G(m)`$ $`=`$ $`Z(m)T_0(m)+U\left[4P(0)+P(1,+)+P(1,)\right]+2J_TP(2,0)`$ (48)
$`+2J_S\left[P(2,t+)+P(2,t0)+P(2,t)\right]\left(3J_S+J_T\right)P(2,\pm ),`$
where

$$T_0(m)=\mathrm{\Phi }_0(m)|\widehat{H}_0|\mathrm{\Phi }_0(m),$$

and

$`\sqrt{Z(m)}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{14m^2}}}[\sqrt{P(0)}(\sqrt{{\displaystyle \frac{P(1,+)}{2}}}+\sqrt{{\displaystyle \frac{P(1,)}{2}}})`$ (49)
$`+{\displaystyle \frac{1}{2}}\left(\sqrt{P(1,+)}+\sqrt{P(1,)}\right)\left(\sqrt{P(2,\pm )}+\sqrt{{\displaystyle \frac{P(2,0)}{2}}}+\sqrt{{\displaystyle \frac{P(2,t0)}{2}}}\right)`$
$`+\sqrt{{\displaystyle \frac{P(2,t+)P(1,+)}{2}}}+\sqrt{{\displaystyle \frac{P(2,t)P(1,)}{2}}}]`$
For any finite $`U`$ the optimal solution has always $`m0`$ due to the nesting property. For very large $`U`$ we expect $`m1/2`$. In this limit we can neglect all $`P`$’s but $`P(1,+)`$ and $`P(2,t+)`$ hence, from Eqs. (46) and (47),

$$P(2,t+)=4m1,P(1,+)=12m,$$

which implies that

$$Z(m)2\frac{4m1}{1+2m}.$$

In the same limit the uncorrelated hopping energy has the expression

$$T_0(m)2\sqrt{M_2\left(24m\right)},$$

where

$$M_2=𝑑ϵ\rho (ϵ)ϵ^2,$$

is the second moment of the uncorrelated density of states per spin and orbital, $`\rho (ϵ)`$. Therefore the variational energy as function of the order parameter $`m`$ for $`U|T_0|`$ is

$$E_G(m)4\frac{4m1}{1+2m}\sqrt{M_2\left(24m\right)}+U\left(12m\right)+2J_S\left(4m1\right),$$

and it is optimized by

$$m\frac{1}{2}\frac{M_2}{\left(U4J_S\right)^2},$$
(50)
leading to

$$E_G2J_S\frac{2M_2}{U4J_S}.$$
(51)
## Appendix B Evaluation of the $`Z`$-factor for the two-band model

The explicit expression of the $`Z`$ reduction factor in the two-band model of Section VI allowing for a superconducting order parameter is, through Eq. (41),

$`\sqrt{Z}`$ $`=`$ $`\lambda _{\overline{0}}\lambda _{\overline{1}}\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)^3+{\displaystyle \frac{3}{2}}\lambda _{\overline{1}}\lambda _{\overline{2}t}\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)^2\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)`$ (52)
$`+`$ $`\lambda _{\overline{1}}\lambda _{\overline{2}\pm }\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)^2\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)+{\displaystyle \frac{1}{2}}\lambda _{\overline{1}}\lambda _{\overline{2}0}\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)^2\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)`$
$`+`$ $`{\displaystyle \frac{3}{2}}\lambda _{\overline{3}}\lambda _{\overline{2}t}\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)^2+\lambda _{\overline{3}}\lambda _{\overline{2}\pm }\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)^2`$
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _{\overline{3}}\lambda _{\overline{2}0}\left({\displaystyle \frac{1}{2}}+\mathrm{\Delta }\right)\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)^2+\lambda _{\overline{3}}\lambda _{\overline{4}}\left({\displaystyle \frac{1}{2}}\mathrm{\Delta }\right)^3.`$
If we parametrize the $`\lambda `$’s according to Eq. (43) and make use of (44) we find

$`\sqrt{Z}`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{14\mathrm{\Delta }^2}}}\{{\displaystyle \frac{1}{4}}\sqrt{P(0)P(1)}[\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{0})P_{\mathrm{\Delta }_0}^{(0)}(\overline{1})}{P_{\mathrm{\Delta }_0}^{(0)}(0)P_{\mathrm{\Delta }_0}^{(0)}(1)}}}+\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{4})P_{\mathrm{\Delta }_0}^{(0)}(\overline{3})}{P_{\mathrm{\Delta }_0}^{(0)}(0)P_{\mathrm{\Delta }_0}^{(0)}(1)}}}`$ (53)
$`+`$ $`\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{2},0)P_{\mathrm{\Delta }_0}^{(0)}(\overline{1})}{P_{\mathrm{\Delta }_0}^{(0)}(0)P_{\mathrm{\Delta }_0}^{(0)}(1)}}}+\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{2},0)P_{\mathrm{\Delta }_0}^{(0)}(\overline{3})}{P_{\mathrm{\Delta }_0}^{(0)}(0)P_{\mathrm{\Delta }_0}^{(0)}(1)}}}]`$
$`+`$ $`{\displaystyle \frac{1}{4}}\sqrt{P(2,0)P(1)}\left[\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{0})P_{\mathrm{\Delta }_0}^{(0)}(\overline{1})}{P_{\mathrm{\Delta }_0}^{(0)}(2,0)P_{\mathrm{\Delta }_0}^{(0)}(1)}}}+\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{4})P_{\mathrm{\Delta }_0}^{(0)}(\overline{3})}{P_{\mathrm{\Delta }_0}^{(0)}(2,0)P_{\mathrm{\Delta }_0}^{(0)}(1)}}}\right]`$
$`+`$ $`[{\displaystyle \frac{\sqrt{3}}{4}}\sqrt{P(2,t)P(1)}+{\displaystyle \frac{\sqrt{2}}{4}}\sqrt{P(2,\pm )P(1)}][\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{1})}{P_{\mathrm{\Delta }_0}^{(0)}(1)}}}+\sqrt{{\displaystyle \frac{P_{\mathrm{\Delta }_0}^{(0)}(\overline{3})}{P_{\mathrm{\Delta }_0}^{(0)}(1)}}}]\}.`$
When $`\mathrm{\Delta }=\mathrm{\Delta }_0=0`$, the above expression reduces to the $`Z`$-factor for a paramagnetic non superconducting solution, namely

$`\sqrt{Z}`$ $`=`$ $`2\sqrt{P(0)P(1)}+\sqrt{3}\sqrt{P(1)P(2,t)}`$ (54)
$`+`$ $`\sqrt{2}\sqrt{P(1)P(2,\pm )}+\sqrt{P(1)P(2,0)}.`$
The uncorrelated probabilities distributions in the natural basis with order parameter $`\mathrm{\Delta }_0`$ (analogous expressions hold for $`\mathrm{\Delta }`$) are

$$\begin{array}{ccc}P_{\mathrm{\Delta }_0}^{(0)}(\overline{0})\hfill & =& \left(\frac{1}{2}+\mathrm{\Delta }_0\right)^414\delta +6\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(\overline{1})\hfill & =& 4\left(\frac{1}{2}+\mathrm{\Delta }_0\right)^3\left(\frac{1}{2}\mathrm{\Delta }_0\right)4\delta 12\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(\overline{2};0)\hfill & =& \left(\frac{1}{2}+\mathrm{\Delta }_0\right)^2\left(\frac{1}{2}\mathrm{\Delta }_0\right)^2\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(\overline{2};\pm )\hfill & =& 2\left(\frac{1}{2}+\mathrm{\Delta }_0\right)^2\left(\frac{1}{2}\mathrm{\Delta }_0\right)^22\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(\overline{2};t)\hfill & =& 3\left(\frac{1}{2}+\mathrm{\Delta }_0\right)^2\left(\frac{1}{2}\mathrm{\Delta }_0\right)^23\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(\overline{3})\hfill & =& 4\left(\frac{1}{2}+\mathrm{\Delta }_0\right)\left(\frac{1}{2}\mathrm{\Delta }_0\right)^30,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(\overline{4})\hfill & =& \left(\frac{1}{2}\mathrm{\Delta }_0\right)^40,\hfill \end{array}$$
(55)
where the last expressions on the left hand side correspond to the limit $`\mathrm{\Delta }_0=1/2\delta `$ with $`\delta 1`$. The uncorrelated occupation probabilities in the original electronic basis are readily obtained by the latter upon inverting Eq. (40):

$$\begin{array}{ccc}P_{\mathrm{\Delta }_0}^{(0)}(0)\hfill & =& P_{\mathrm{\Delta }_0}^{(0)}(4)=\left(\frac{1}{4}+\mathrm{\Delta }_0^2\right)^2\frac{1}{4}\delta +2\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(1)\hfill & =& P_{\mathrm{\Delta }_0}^{(0)}(3)=4\left(\frac{1}{16}\mathrm{\Delta }_0^4\right)2\delta 6\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(2,0)\hfill & =& \frac{1}{16}+\frac{3}{2}\mathrm{\Delta }_0^2+\mathrm{\Delta }_0^4\frac{1}{2}2\delta +3\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(2,\pm )\hfill & =& 2\left(\frac{1}{4}\mathrm{\Delta }_0^2\right)^22\delta ^2,\hfill \\ P_{\mathrm{\Delta }_0}^{(0)}(2,t)\hfill & =& 3\left(\frac{1}{4}\mathrm{\Delta }_0^2\right)^23\delta ^2.\hfill \end{array}$$
(56)
In the limit of small $`\delta `$, by inserting (55) and (56) into (53) one finds at leading order (recalling that $`\mathrm{\Delta }Z1`$)

$`\sqrt{Z}`$ $``$ $`\sqrt{2}\sqrt{P(0)P(1)}\left(1+\delta \right)+\sqrt{P(1)P(2,0)}`$ (57)
$`+`$ $`\left[\sqrt{{\displaystyle \frac{3}{2}}}\sqrt{P(1)P(2,t)}\sqrt{P(1)P(2,\pm )}\right]\left(1+\delta \right).`$


text:
# Chapter 0 AN XMM and Chandra view of massive clusters of galaxies to z=1

## Chapter 0 AN XMM and Chandra view of massive clusters of galaxies to z=1

### 1 Introduction

Massive clusters of galaxies are rare objects, forming from the high-sigma tail of the cosmological density distribution. Their properties are powerful probes of cosmology, and give insight into the process of structure formation on large scales. For example, compared to z=0, and in a $`\mathrm{\Lambda }`$CDM Universe, the number density of virialized halos of $`M>5`$x1014 $`M_{}`$ is predicted to be a factor of $``$10 lower at z=0.5, or a factor of $``$100 lower at z=1, from both simulations and the Press-Schechter approximation (eg. Bode et al.2001). Such a large change in number density implies that a significant fraction of the progenitor systems at high redshift should be unvirialized, perhaps containing lower mass systems in the process of merging.

The number of genuinely massive clusters known at high redshifts is very small, and the number of relaxed high redshift massive clusters is even smaller. Thus the fraction which are relaxed, giving direct information on the epoch and mode of assembly of massive clusters, is poorly known. Relaxed, massive clusters at high redshifts also offer the best opportunities for deriving total masses and gas mass fractions from the X-ray data without too many uncertainties due to hydrostatic equilibrium assumptions, and also with reasonable signal-to-noise.

Detailed analyses of 3 of the clusters presented here can be found in Maughan et al.(2003a,b). We use $`\mathrm{\Omega }_m`$=0.3, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7 and H0=70 km s-1 unless otherwise stated.

### 2 The cluster sample

The clusters were discovered in a $`ROSAT`$ serendipitous X-ray survey (WARPS - Scharf et al.1997, Jones et al.1998, Perlman et al.2002). Of 16 clusters at z$`>`$0.6 which form a complete, X-ray selected sample, 10 have been observed with $`XMM`$ and/or $`Chandra`$. For 3 of these, the $`XMM`$ observations were of low quality due to periods of high background. The remaining observations, however, represent a major step forward compared to previous results from $`ROSAT`$ (see also Rosati et al., these proceedings and references therein).

### 3 Analysis of the high redshift clusters

Standard analysis techniques were applied to the $`XMM`$ and $`Chandra`$ data, taking into account the ACIS contamination affecting the low energy efficiency, and using a local background subtraction method to obtain spectra. Images were exposure corrected and adaptively smoothed so that all features are significant at the 99% confidence level.

In general a single temperature was measured (but see CLJ1226.9+3332 below for a temperature profile). Two-dimensional $`\beta `$-profiles were fit to the images, including PSF blurring, to obtain values for the core radii and $`\beta `$. For the clusters which appeared relaxed, and assuming spherical symmetry, the gas mass was derived from an integral of the density profile given by the $`\beta `$-profile. The total gravitating mass was derived by additionally assuming isothermality and hydrostatic equilibrium. The total mass within a radius $`r`$ is then given by

$`M(<r)`$ $`=`$ $`1.13\times 10^{14}\beta {\displaystyle \frac{T}{keV}}{\displaystyle \frac{r}{Mpc}}{\displaystyle \frac{(r/r_c)^2}{1+(r/r_c)^2}}M_{}.`$ (1)
To measure the virial radius, we used eqn (1.1) to find the radius within which the total overdensity was 200 times the critical density at the redshift of each cluster. The quoted fluxes, luminosities and masses were measured via an extrapolation from the typical detection radius of $``$0.3-0.45 $`r_{vir}`$ to the virial radius.

### 4 Descriptions of the clusters

#### 1 A massive, relaxed cluster at z=0.89

The discovery of CLJ1226.9+3332 is described by Ebeling et al.(2001). The $`XMM`$ image of CLJ1226.9+3332 is shown in Fig 1, overlaid on a deep I band image. The X-ray morphology appears generally relaxed, in contrast to cluster MS1054-0321, the other high temperature cluster known at z$`>`$0.8 (Jeltema et al.2001). The $`XMM`$ temperature of 11.5$`{}_{0.9}{}^{}{}_{}{}^{+1.1}`$ keV is consistent with the velocity dispersion of $``$1100 km s-1, based on 15 galaxy redshifts. The bolometric X-ray luminosity is 5.4x1045 erg s-1. We derive a mass of (1.4$`\pm `$0.2)x1015 $`M_{}`$ within the virial radius. Our analysis of a short $`Chandra`$ observation confirms the lack of significant point-source contamination, and the temperature and luminosity, albeit with lower precision (see also Cagnoni et al.2001).

A temperature profile (Fig 1) is consistent with the cluster being isothermal out to 45% of the virial radius. The metal abundance of $`Z`$=0.33$`{}_{0.10}{}^{+0.14}Z_{}^{}`$ and gas mass fraction of 8.6$`{}_{1.0}{}^{}{}_{}{}^{+1.1}`$% are consistent with those of local clusters of the same temperature (see below).

This cluster is unique in being at high redshift, yet massive and generally relaxed. It must have been assembled when the age of the Universe was significantly less than 6 Gyr. Further details of the $`XMM`$ analysis and the cosmological implications are given in Maughan et al.(2003b).

#### 2 A massive cluster in formation at z=0.83

In stark contrast is the highly unrelaxed, massive cluster ClJ0152.7-1357, described in Ebeling et al.(2000), Della Ceca et al.(2000) and Maughan et al.(2003a). The $`Chandra`$ image and a recently acquired deep $`XMM`$ image are shown in Fig 2. Based on the $`Chandra`$ data, we find that the cluster consists of two major sub-clumps at very similar redshifts and of temperatures 5.5$`{}_{0.8}{}^{}{}_{}{}^{+0.9}`$ keV and 5.2$`{}_{0.9}{}^{}{}_{}{}^{+1.1}`$ keV. The total mass of each subcluster is $``$(5-6)x1014 $`M_{}`$ . A dynamical analysis of the system shows that the subclusters are likely to be gravitationally bound. When merged, the system mass will be similar to that of the Coma cluster.

The deeper $`XMM`$ image shows lower surface brightness features. To the east is a low luminosity system also at the cluster redshift (Demarco et al.2003, in preparation). To the NW is a possible filamentary structure with a point source and one or two low luminosity extended sources, probably groups, embedded within it. This structure falls within a radius of 1.6 Mpc (the dashed circle), the estimated virial radius of the final merged system, suggesting that it is part of the same formation process. The faintest filamentary emission is significant at $`>`$5$`\sigma `$ in the soft 0.3-1 keV band and has a hardness ratio consistent with a temperature of 1$`\pm `$0.7 keV. We may be observing the warm baryons predicted to lie within filaments, based on simulations of structure formation.

#### 3 Other clusters at z$`>`$0.6

$`Chandra`$ and $`XMM`$ contours of the remaining 5 clusters are shown in Figs 3, 4 and 5 in order of the most relaxed first. There is a wide range of morphologies, and some (eg ClJ1559) clearly have bright point X-ray sources nearby.

An additional cluster, ClJ1227.3+3333 (z=0.766) is shown in Fig 6. This is an extremely unrelaxed system. Although not in the original $`ROSAT`$ sample (it fell below a window support shadow), it is detected in an $`XMM`$ image and is a factor $``$12 less luminous than ClJ1226.9+3332. There is no clear optical or X-ray core, but rather a system of $``$4 X-ray subclumps separated by $``$0.5 Mpc which are at a very similar redshift. There is a red sequence in the galaxy colours. This is probably a highly non-virialized system seem in formation, yet it is detected in X-rays.

### 5 Metallicities and gas mass fractions

Most of the metals in a cluster are in the X-ray gas, and so the evolution of the metallicity of the gas provides important information on the chemical history of the Universe. In Fig 7 we show that our measurements, and those of others, of cluster metallicities up to z=0.9 are consistent with the canonical value of 0.3 times the solar value found at lower redshifts by Mushotzky & Loewenstein (1997).

The gas mass fractions we measure, within a radius of 0.3$`r_{vir}`$ (where they are the most reliable) are consistent with those found by Sanderson et al.(2002) for local clusters of the same temperature. Thus the gas, and the metals it contains, were in place by z$``$0.8 in at least some massive clusters.

### 6 The evolution of the mass-temperature relation

In Fig 8 we show the mass-temperature relation for those high redshift clusters which appear reasonably relaxed (shown as red squares). The circles are the low redshift data points of Sanderson et al.(2002). In both cases the mass is that measured within the virial radius, $`r_{200}`$, and in both studies $`r_{200}`$ is measured from the derived mass profile as the radius within which the density is 200 times the critical density at the redshift of observation. The low and high redshift samples have different treatments of the temperature profile, however. For the high redshift sample we (necessarily) assume isothermality. For the low redshift sample, the temperature profile is measured and thus a more accurate mass is derived.

The upper dashed line shows the best fit to the low redshift points. The lower (red) dashed line shows the expected change in normalisation assuming clusters observed at z=0.8 were formed at z=0.8, using the relation $`E(z)MT^{3/2}`$, where $`E(z)=(1+z)\sqrt{(1+z\mathrm{\Omega }_M+\mathrm{\Omega }_\mathrm{\Lambda }/(1+z)^2-\mathrm{\Omega }_\mathrm{\Lambda })}`$. It is interesting to note that all the high redshift clusters except the most massive one are consistent with the low redshift relation, rather than the z=0.8 prediction. It is however the most massive cluster (ClJ1226) which has evidence of isothermality and the most reliable mass. The masses of the other high redshift clusters would need to be reduced by $``$40% in order to be consistent with the z=0.8 prediction. Sanderson et al.have shown that an incorrect assumption of isothermality when clusters are not isothermal produces, on average, masses which are too high by 30% at $`r_{200}`$. This would explain most of the offset we find. Further work is in progress, and will be reported in future papers, but more accurate masses at high redshifts, based on temperature profiles, may be required before firm conclusions can be reached.

### 7 Point source contamination of $`ROSAT`$ luminosities

The $`ROSAT`$ detections were obtained at $`ROSAT`$ off-axis angles of up to 15 arcmin, where the PSF degrades to $``$50 arcsec (FWHM). Thus some point-source contamination is expected in all the $`ROSAT`$ cluster surveys. The degree of contamination depends on the signal-to-noise ratio of the detection and the algorithm used for source detection and characterisation.

Based on a preliminary analysis, 2 or 3 out of 9 clusters at z$`>`$0.6 in the WARPS survey were artificially boosted above the survey flux limit by significant point source contamination. Removing these sources from the sample, and using the updated $`XMM`$ and $`Chandra`$ fluxes, we recomputed the X-ray luminosity function at 0.6$`<z<`$1.1. The resultant shifts to lower space densities were within the Poissonian error bars, suggesting that point source contamination can be an important effect for individual $`ROSAT`$ measurements, but the previous broad conclusions regarding cluster evolution are unaffected.

### 8 References

Bode, P., Bahcall, N. A., Ford, E. B., Ostriker, J. P., 2001, ApJ, 551, 15
Cagnoni, I., Elvis, M., Kim, D.-W., Mazzotta, P., Huang, J.-S., Celotti, A., 2001, ApJ, 560, 86
Della Ceca, R., Scaramella, R., Gioia, I.M., Rosati, P., Fiore, F., Squires, G., 2000, A&A, 353, 498
Ebeling, H., Jones, L. R., Perlman, E. S., Scharf, C. A., Horner, D., Wegner, G., Malkan, M., Mullis, C.R. 2000, ApJ, 534, 133.
Ebeling, H., Jones, L. R., Fairley, B. W., Perlman, E., Scharf, C., Horner, D. 2001, ApJ 548, L23.
Jeltema, T. E., Canizares, C. R., Bautz, M. W., Malm, M. R., Donahue, M., Garmire, G. P., 2001, ApJ, 562, 124
Jones, L.R., Scharf, C.A., Ebeling, H., Perlman, E., Wegner, G., Malkan, M., Horner, D., 1998, ApJ, 495, 100
Maughan, B.J., Jones, L.R., Ebeling, H., Perlman, E., Rosati, P., Frye, C., Mullis, C.R., 2003a, ApJ in press (astro-ph/0301218)
Maughan, B.J., Jones, L.R., Ebeling, H., Scharf C., 2003b, in preparation
Mushotzky, R. F., Loewenstein, M., 1997, ApJ, 481, L63
Perlman, E.S., Horner, D., Jones, L.R., Scharf, C., Ebeling, H., Wegner, G., Malkan, M., 2002, ApJS, 140, 265
Scharf, C., Jones, L.R., Ebeling, H., Perlman, E., Wegner, G., Malkan, M., 1997, ApJ, 477, 79

### 9 Acknowledgements

We acknowledge useful discussions with Stefano Ettori, Gus Evrard, Monique Arnuad & Trevor Ponman.