{'fraction_non_alphanumeric': 0.03965330444203684, 'fraction_numerical': 0.010698808234019502, 'mean_word_length': 4.397017979827511, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gutzwiller description of non-magnetic Mott insulators: a dimer lattice model

We introduce a novel extension of the Gutzwiller variational wavefunction able to deal with insulators that escape any mean-field like description, as for instance non-magnetic insulators. As an application, we study the Mott transition from a paramagnetic metal into a non-magnetic Peierls, or valence-bond, Mott insulator. We analyze this model by means of our Gutzwiller wavefunction analytically in the limit of large coordination lattices, where we find that: (1) the Mott transition is first order; (2) the Peierls gap is large in the Mott insulator, although it is mainly contributed by the electron repulsion; (3) singlet-superconductivity arises around the transition.


I. INTRODUCTION
Among the theoretical tools devised to deal with strongly correlated metals close to a Mott metal-toinsulator transition (MIT), the simplest one likely is the variational approach introduced in the 60ths' by Gutzwiller 1,2 to describe itinerant ferromagnetism and narrow band conductors. In its original version, the Gutzwiller variational wavefunction has the form where |φ is an uncorrelated wavefunction for which Wick's theorem holds, P R an operator at site R, and both |φ and P R have to be determined by minimizing the variational energy. The role of the operator P R is to modify, according to the on-site interaction, the weights of the local electronic configurations with respect to their values in the uncorrelated wavefunction. In spite of its simplicity, the Gutzwiller wavefunction is quite effective in capturing physical properties that supposedly identify strongly correlated metals, as for instance the large increase of the effective mass. 3 However, since the dependence upon the distance |R − R ′ | of intersite correlations are still determined by the uncorrelated wavefunction, while the local operators P R just affect the amplitudes, the Gutzwiller wavefunction can describe a Mott insulator either if P R suppresses completely charge fluctuations, that provides a very poor description of an insulator, or if |φ itself is insulating. The latter case can be stabilized within the original Gutzwiller approach only when |φ is an admissible Hartree-Fock solution of the Hamiltonian. As an example let us consider a single band model at half-filling, for instance the Hubbard model where c † R,σ and c R,σ creates and annihilates, respectively, an electron with spin σ =↑, ↓ at site R and n R,σ = c † R,σ c R,σ . This Hamiltonian admits at the mean-field level two possible phases, one paramagnetic, n R,↑ = n R,↓ , and the other magnetic, n R,↑ = n R,↓ . The latter is the only one that can eventually describe an insulator. In the Hubbard model the action of the operator P is to increase the weight of singly occupied sites at expenses of doubly occupied and empty sites, in order to minimize the Coulomb repulsion U . Evidently, even when the repulsion is very strong, hence the model is a Mott insulator, a realistic wavefunction should still allow for charge fluctuations responsible for the superexchange, that survives even deep inside the Mott phase. However, any paramagnetic uncorrelated wavefunction, for instance the Fermi sea, is unable to generate any super-exchange and necessarily leads to a non-realistic Mott insulator where configurations with empty or doubly occupied sites are fully suppressed. 3 The only way to generate super-exchange is to assume a magnetically ordered |φ , which is also the only insulating wavefunction accessible within Hartree-Fock. However a magnetic state might not always be the right choice, especially if magnetism is sufficiently frustrated.
Recently, an improved version of the Gutzwiller wavefunction has been proposed, 4 in which additional intersite correlations are provided by density-density Jastrow factors, namely where n R is the site R occupation number and v R,R ′ variational parameters. This novel class of wavefunctions has the capability to disentangle charge from other degrees of freedom, hence is more suitable to capture Mott localization, as it has indeed been shown. 4,5,6 However, unlike the conventional Gutzwiller wavefunction (1), the Gutzwiller-Jastrow wavefunction (2) can only be dealt with numerically by variational Monte Carlo, which is inherently limited to finite-size systems, albeit quite large. 7 An alternative approach, that is closely related to recently proposed extensions of Dynamical Mean Field Theory (DMFT) from the original single-site formula-tion 8 to a cluster one 9,10,11,12,13 , is to consider a variational wavefunction of the same form as (1) but defined on a lattice with non-primitive unit cells. In this case, the operator P R acts on all the available electronic configurations of the lattice sites belonging to the nonprimitive cell. The advantage is that in this way one may include additional short-range correlations without losing the property of the wavefunction to be analytically manageable, at least in infinite-coordination lattices. The obvious disadvantage is that this wavefunction could bias the variational solution towards translational-symmetry breaking.
Within this scheme, the variational problem becomes generically equivalent to optimize a Gutzwiller wavefunction for a multi-band Hamiltonian. There have been recently an amount of attempts to extend the Gutzwiller wavefunction to multi-orbital models that include further complications like for instance Coulomb exchange 14,15,16,17 . In this paper we introduce a further extension that is capable to generate inter-site correlations as the super-exchange for paramagnetic wavefunctions, otherwise missed by the conventional Gutzwiller approach. This novel class of wavefunctions also allows to explore new kinds of variational solutions. Specifically, there are interesting examples of correlated models where the Mott insulating phase escapes any Hartree-Fock mean-field treatment, in other words can not be represented by a single Slater determinant. A very simple case, that we will explicitly consider throughout this work, is a Peierls insulator, namely a short-range valencebond crystal, in which pairs of nearest neighbor sites are strongly bound into a singlet configuration, leading to a state that is simply a collection of spin-singlets. Such a Mott insulating state is not accessible by Hartree-Fock theory, just because each singlet is itself not expressible as a Slater determinant, nor by the conventional Gutzwiller approach, which, as mentioned, gives a poor description of paramagnetic insulators.
The paper is organized as follows. In Section II we present the variational wavefunction and discuss under which conditions it can be deal with analytically. In Section III we discuss how to build up the wavefunction in the case in which the basic unit of the lattice model is a dimer. Next, in Section IV, we solve the variational problem for a specific lattice model of dimers. Conclusions are given in Section V.

II. THE VARIATIONAL WAVEFUNCTION
In this Section, we introduce an extension of the Gutzwiller wavefunction (1) which is particularly convenient to perform analytical calculation in the limit of infinite-coordination lattices. 14 Let us consider a generic multi-band Hamiltonian. Each lattice site R contains several orbitals that give rise to a bunch of electronic configurations which we denote individually as |Γ; R . The most general operator P R can be chosen of the form: where λ(R) ΓΓ ′ are variational parameters. In general P R needs not to be hermitean, namely for Γ = Γ ′ it is not required that λ(R) * ΓΓ ′ = λ(R) Γ ′ Γ . Indeed, as we shall see, the non-hermitean character plays a very important role. We further assume that the Wick's theorem holds for the uncorrelated wavefunction, hence that |φ is either a Slater determinant or a BCS wavefunction.
It was realized by Bünemann, Weber and Gebhard 14 that average values of operators on the Gutzwiller wavefunction (1) can be analytically computed in infinite coordination lattices provided the following two constraints are imposed on P R : where C R is the local single-particle density-matrix operator, with elements c † R,α c R,β and c † R,α c † R,β , α labeling single-particle states, while c † R,α and c R,α create and annihilate, respectively, an electron at site R in state α.
The first constraint, Eq. (4), does not actually limit the variational freedom, since P R is defined up to a normalization factor. On the contrary, the latter constraint, Eq. (5), may reduce the variational freedom, although it seems not in a relevant manner, at least in all cases that we have so far investigated. We notice that Eq. (5) is not the same as imposing unless P R commutes with C R , which is a further constraint to be imposed on P R . This actually is the only case that has been hitherto considered, see e.g. Refs. 14 and 15. However, as we shall see, there are interesting models which force to abandon the supplementary condition (6), which is anyway unnecessary. 17 By means of Wick's theorem, the left-hand side of (5) includes a disconnected term φ| P † R P R |φ φ| C R |φ = φ| C R |φ , where the right-hand side follows from (4), plus connected terms that are obtained by selecting in all possible ways a pair of single-fermion operators from P † R P R , averaging on |φ what remains, and finally averaging the two single-fermion operators with those of C R . Therefore, imposing (5) means that the sum of all connected terms vanishes, whatever is the element of the singleparticle density-matrix. In other words, the operator that is left after taking out from P † R P R any pair of singlefermion operators has null average on |φ . In turns, this also implies that, when averaging on |φ P † R P R with multi-particle operators at different sites, the only connected terms that survive are those that involve four or Graphical representation of the average on |φ of P † R P R , drawn as a box, times a generic multi-particle operator at site R ′ , drawn as a circle. Lines that join the two operators represent the average of two single-fermion operators, one at R and the other at R ′ . The dots include all terms where the two sites are joined by more than four lines. The important thing to notice is the absence of terms in which the two sites are connected by two lines. more single-fermion operators of P † R P R , that are represented graphically in Fig. 1 as lines coming out of P † R P R . This property of P R turns out to be extremely useful in infinite-coordination lattices. In this limit, the contribution to the average value on |φ of terms in which more than two fermionic lines come out of P † R P R can be shown to vanish 14 , which simplifies considerably all calculations. For instance, the average value on (1) of any local operator O R becomes which also implies, taking O R = 1, that the variational wavefunction (1) is normalized. In addition, the average value of the inter-site single-particle density matrix turns out to be where the matrices Z and ∆ are determined by inverting the following set of equations Naïvely speaking, it is as if, when calculating the intersite density matrix, a fermionic operator transforms effectively into namely that a particle turns into a particle or a hole with probabilities Z and ∆, respectively. Although all the above expressions are strictly valid only in infinitecoordination lattices, it is quite common to use the same formulas also to evaluate average values on the Gutzwiller wavefunction in finite-coordination lattices. This approximation is refereed to as the Gutzwiller approximation 1,2,18,19,20 , and is known to be equivalent to the saddle point solution within the slave-boson technique. 21 We conclude by noting that the constraint (5) turns out to be useful also when the variational wavefunction (1) is applied to Anderson impurity models. In this case the operator P acts only on the electronic configurations |Γ of the impurity, namely where C imp is the single-particle density matrix of the impurity, then, for any operator of the conduction bath, O bath , and because of (4) and (5), the following result holds A. Some formal definitions In order to perform actual calculations, it is convenient to introduce some notations. We define a matrix F α with elements as well as its hermitean conjugate, F † α , where we assumed that the definition of the local configurations is the same for all sites. It follows that where I is the identity. Next, we introduce the uncorrelated occupation-probability matrix, P 0 (R), with elements that satisfies Analogously, the variational parameters that define P R , λ(R) Γ1Γ2 , are interpreted as elements of a matrix λ(R).
With these definitions, Eqs. (4) and (5) become that suggests to introduce a variational occupationprobability matrix P (R) = P 0 (R) λ(R) † λ(R) with matrix elements that must satisfy Eqs. (20), (21) and (22) replace the constraints (4) and (5). With these definitions, the matrices Z and ∆, see Eqs. (11) and (13), are obtained by solving The above equations simplify if one uses the natural basis, namely the single-particle basis that diagonalizes the density-matrix, Moreover, if one constructs the states |Γ; R so that P 0 (R) is diagonal

B. The variational problem
We are now in position to settle up the variational problem. We consider a generic tight-binding Hamiltonian where α and β stem for spin, orbital and lattice site in the chosen unit cell, and the hermitean matrix E(R) with elements E(R) ΓΓ ′ may be also unit-cell dependent. The average value of this Hamiltonian on the Gutzwiller wavefunction (1) in the limit of infinite coordination lattices or, in finite coordination ones, within the Gutzwiller approximation, is The last term depends only on the local properties of the uncorrelated wavefunction |φ , specifically on the occupation probabilities P 0 (R). Therefore, for any given choice of P 0 (R), the optimal |φ that minimizes the variational energy is the ground state of the Hamiltonian where the parameters µ(R) αβ and ν(R) αβ are Lagrange multipliers to be determined by imposing that the ground state has indeed the chosen P 0 (R). The last task is to find the values of the variational parameters λ(R) ΓΓ ′ as well as of P 0 (R) for which the variational energy (29) is minimum. We note that the variational Hamiltonian (30) that has to be solved may include also inter-site pairing terms, which are absent in the original Hamiltonian (28). Analogously to other more conventional variational approaches, like Hartree-Fock theory, it is common to interpret the single-particle spectrum of the variational Hamiltonian (30) as an approximation of the true coherent spectrum of quasi-particles. 22

III. THE MODEL
Let us now apply the variational wavefunction to specific models that are inspired by the valence-bond crystal example we mentioned in the introduction, and where the off-diagonal elements of the operator P R as well as its non-hermitean character do play an important role. Since the operator P R is built out of purely local properties, namely the available on-site electronic configurations plus a variational guess for the uncorrelated on-site single-particle density-matrix, a lot of preliminary results can be extracted without even specifying how lattice-sites are coupled together. Therefore we start our analysis from defining some local properties and later we will consider a specific lattice model.

A. The isolated dimer
The basic unit of the model we are going to investigate consists of a dimer with Hamiltonian where 1 and 2 refer to the two sites of the dimer and n i , i = 1, 2, is the on-site occupation number. It is more convenient to work in the basis of the even (bonding) and odd (anti-bonding) combinations, defined through and use this basis to built the available electronic configurations, which we will denote as |n, Γ , with n that refers to the number of electrons. The empty and the fourfold occupied dimer states are denoted as |0 and |4 , respectively, while the singly-occupied states as and the states with 3 electrons as There are six doubly-occupied configurations. Two are spin-singlets with two electrons in the even or in the odd orbital, |2, e and |2, o , respectively. When each orbital is singly occupied, the two electrons form either a spin triplet, |2, 1, S z with S z = −1, 0, 1, or a spin singlet, |2, 0 . Since we are not going to consider variational solutions that break spin-SU (2) symmetry, it is convenient to define the projector operators The isolated-dimer ground state in the subspace with two electrons is where tan 2θ = 4t ⊥ /U and has energy |Ψ can be always rewritten in the form of a Gutzwiller wavefunction. First of all, we needs to choose an uncorrelated wavefunction |φ . A natural choice might be the ground state at U = 0, namely |2, e . Indeed |Ψ can be written as and obviously satisfies both (4) and (5).
Another possibility, that we are also going to consider in what follows, is to use an uncorrelated wavefunction that corresponds to a dimer in which the two sites are only coupled by an intersite singlet-Cooper pairing, and, once again, the true ground state can be written as |Ψ = |Ψ φ| |φ ≡ P |φ .
Already at this stage one can appreciate how important is the role of the off-diagonal elements in P, especially for large U/t ⊥ .

B. The non-isolated dimer: variational density matrix
When the dimer is coupled to the rest of the system, in order to built the operator P we need to specify an uncorrelated local single-particle density matrix based on a variational guess of the uncorrelated wavefunction |φ . A simple guess would be a magnetic wavefunction in which the two sites of each dimer have opposite magnetization. This choice is also the only one admitted by an Hartree-Fock decomposition of the interaction term H U . However, a magnetic wavefunction is not the most suitable choice to reproduce the limit of isolated dimers, which is a collection of singlets.
Alternatively, one can consider a paramagnetic |φ that has built in the tendency of each dimer to lock into a spin-singlet. This can be accomplished in two ways that do not exclude each other. The first is to assume an uncorrelated wavefunction with a huge splitting between even and odd orbitals, namely with This implies that, among the doubly-occupied configurations of each dimer, mainly the spin-singlet |2, e survives in the uncorrelated wavefunction. The latter can then be turned into the isolated dimer configuration by an appropriate Gutzwiller operator P, as shown before. The other possibility is to include Cooper pairing correlations in the singlet channel |φ . In this case, the isolated dimer can be recovered by assuming a very strong pairing ∆ SC ≃ 1 and suppressing, through P, configurations with none or two singlet-pairs. Note that both n e − n o and ∆ SC do not appear by a mean-field decoupling of H U , so that a variational wavefunction with such correlations built in can not be stabilized within Hartree-Fock theory. Here the role of P becomes crucial.
Therefore, let us assume for |φ a BCS-wavefunction defined such that with real ∆ SC . In the even/odd basis this translates into where n e + n o = n. As previously mentioned, the calculations simplify considerably in the natural basis, that is derived in the Appendix for this particular choice of density matrices. As a particular application, we assume hereafter that the model is half-filled, namely n e + n 0 = 2. The density matrix of the operators in the natural basis, The two angles θ e and θ o , that are defined by Eq. (A.2) and identify the unitary transformation from the original to the natural basis, are given by θ e = θ and θ o = θ−π/2, where tan 2θ = ∆ SC δ .
We note that, for q → 1/2, the uncorrelated wavefunction describes an insulator where charge fluctuations are completely suppressed since each natural orbital is fully occupied. It is obvious that, if our choice of the variational wavefunction is correct, then the optimal uncorrelated wavefunction must asymptotically acquire q = 1/2 for U → ∞.
The expression in the natural basis of the hopping, Eq. (A.6), and interaction, Eq. (A.7), operators can be derived through (A.5) and have a relatively simple expression at half-filling: where we have defined and denoted the local configurations in the natural basis as | n, Γ to distinguish them from the analogous ones in the original representation.

C. The Gutzwiller operator P
The most general Gutzwiller operator P should include at least all the projectors | n, Γ n, Γ| as well as all the off-diagonal operators | n, Γ n ′ , Γ ′ | that appear in the local Hamiltonian, Eqs. (42) and (43). As we mentioned before, our expectation is that the uncorrelated wavefunction which better connects to the large-U Mott insulator should be identified by q → 1/2, in which locally only the configurations | 3 and | 4 are occupied with nonnegligible probability. This suggests that P must include at least those off-diagonal operators that would turn | 4 into the isolated dimer ground state, namely | 0 4| and | 2, + 4|. The latter forces to include also | 1 3|, as we shall see.
Therefore we assume for P the following variational ansatz: P = e nΓ λ nΓ | n, Γ n, Γ| + λ 13 | 1 3| with real λ's. We define P (n, Γ) = λ 2 nΓ P 0 ( n, Γ), for all n = 3, 4, while, for n = 3, 4, Then the conditions Eqs. (4) and (5) Here P 0 ( n, Γ) are the occupation probabilities in the natural basis of the uncorrelated wavefunction. Specifically where g e n,Γ is the degeneracy of the configuration. Eq. (50) guarantees that the anomalous averages vanish in the natural basis, and explains why we have included | 1 3| in (44). It is convenient to rewrite λ i = P (4) P 0 ( 4) u i for i = 4, 04, 2 + 4, where u 2 3 + u 2 13 = 1, which can be satisfied by choosing u 3 = cos ψ and u 13 = sin ψ, and u 2 4 +u 2 04 +u 2 2+ 4 = 1. The latter parameters can be expressed by means of another In terms of all the variational parameters, the P (n, Γ)'s, θ, q, ψ and v, the average values per dimer of the interaction, H U , and intra-dimer hopping, H ⊥ , are readily found to be where the actual correlated values of the hybridization and of the anomalous average are We note that δ * and ∆ * are mutually exclusive, namely the choice of parameters that maximizes one of the two, makes the other vanishing. Upon the action of P, the single fermion operators in the Nambu spinor representation transform effectively into where τ i , i = 1, 2, 3, are the Pauli matrices that act on the Nambu spinor components, and, finally, with real √ Z and √ ∆. Therefore, if the dimers are coupled one to another by the single particle hopping term where Ψ † R,i = (c † R,i↑ , c R,i↓ ), and Ψ its hermitean conjugate, the uncorrelated wave function |φ minimizes the effective hopping under the condition that the local density matrix satisfies Eqs.(37)-(39). One can readily show that this amounts to find the ground state |φ of the variational Hamiltonian with µ 3 and µ 1 such that φ| H var |φ + 4q µ 3 cos(2θ + 2β) + 4q µ 1 sin(2θ + 2β), is maximum.
Before we consider specific lattice models, it is worth re-deriving within this variational scheme the isolateddimer ground-state energy (32) at half-filling. For that purpose, we take all P (n, Γ) zero but P (4) = 1. The variational energy is simply The minimum under the constraint v · v = 1 is obtained for v 3 = 0 and exactly reproduces (32). We note that the minimum energy is independent on θ, namely there exists a continuous family of variational solutions with equal energy parametrized by θ. However, in spite of the fact that the uncorrelated wavefunction may describe a superconductor, the actual value of the anomalous average ∆ * = 0.

IV. A LATTICE MODEL OF DIMERS
As a particular application, let us consider the following lattice model where n R,i = σ c † R,iσ c R,iσ and ǫ k is the band dispersion induced by t RR ′ , with half-bandwidth D. The Hamiltonian (64) represents two Hubbard models coupled by a single-particle hopping t ⊥ , each model being defined on a lattice with coordination number z. As we mentioned, the variational results that we have so far derived are rigorous strictly speaking only if z → ∞, although, in the spirit of the Gutzwiller approximation, they can be used for generic z as well.
If U ≫ D, t ⊥ , (64) describes at half-filling a Mott insulator which may be magnetic at t ⊥ ≪ D, but is certainly non-magnetic at t ⊥ ≫ D, where the ground state reduces essentially to a collection of singlets. For instance, in the case of a Bethe lattice with nearest neighbor hopping, the transition is at t ⊥ = D/ √ 8, value that is going to decrease if frustration is included. If U is small and the Fermi surface is not nested, then the model is metallic for t ⊥ ≤ D and is a band insulator otherwise. In fact, in the absence of nesting there is generically a finite window of t ⊥ values in which, upon increasing U , the model undergoes a transition from a paramagnetic metal into a non-magnetic Mott insulator, and this is just the case we are going to consider in what follows. The same model have been recently studied by Fuhrmann, Heilmann and Monien using DMFT 23 and by Kancharla and Okamoto 24 using DMFT and cluster DMFT, respectively, that gives us the opportunity to directly check the accuracy of our wavefunction.
The variational Hamiltonian (63) of the model (64) has a very simple expression, The variational single-particle spectrum has the conventional BCS form with eigenvalues hence, for any µ 1 = 0, has a gap equal to 2µ 1 . On the contrary, when µ 1 = 0, the spectrum is gapless for |µ 3 | ≤ D, otherwise is gaped. The Lagrange multipliers µ 1 and µ 2 are obtaining by maximizing In terms of (66), (52) and (53) the variational energy is and depends on eight independent variational parameters.
We have solved numerically the variational problem at fixed t ⊥ /D = 0.5 as function of U/D. To simplify calculations, we have assumed for the band dispersion ǫ k either a flat or a semi-circular density of states, although both would give rise to nesting that could stabilize magnetic phases, which we do not take into account. However, from the point of view of the paramagnetic-metal to paramagnetic-insulator transition, this choice is not influential.
We find that the variational solution displays a first order phase transition at U c ≃ 2.05 D for a flat density of states, as shown by the behavior of the variational energy in Fig. 2. This result agrees almost quantitatively with the DMFT calculation 23 obtained with a semicircular density of states, that also predicts a first order transition with a coexistence region between U ≃ 1.5 D and 1.8 D at the same value of t ⊥ = 0.5 D. We note that the energy is everywhere finite and vanishes like 1/U for large U , see the inset of Fig. 2. The asymptotic behavior U E var /D 2 ∼ −7/6 is compatible with second order perturbation theory in t and t ⊥ using as zeroth-order state a collection of dimers, as explained below. In Fig. 3 we show the behavior across the transition of the three contribution to the energy, namely E U , E ⊥ and E hop . We find that the transition is accompanied by an energy loss in E hop , but a gain in both E ⊥ and E U .
In order to characterize physically the two phases, in Fig. 4 we plot the values of µ 1 and µ 3 across the transition. Since µ 1 = 0, within our numerical accuracy, and |µ 3 | < D, the phase at U < U c is gapless hence metallic, see the behavior of the density of states (DOS) drawn in Fig. 5. On the contrary, on the U > U c side of the transition, µ 1 = 0, that implies a finite gap in the single-particle variational spectrum, see Fig. 5. In the gaped phase at U > U c the spectrum looks like the one of a Peierls insulator with a very large hybridization gap, not consistent with the bare value of t ⊥ . In reality this gap is, more properly, the Mott-Hubbard gap. Indeed the DOS has weight both below and above the chemical potential, suggestive of asymmetric Mott-Hubbard sidebands. Moreover, as we are going to discuss below, the actual difference between the occupations of the bonding and anti-bonding bands, which we denoted as 2δ * in Eq. (54), decreases with U , unlike the single-particle gap, see Fig. 7. This behavior is reminiscent of what has been found by Biermann et al. 25 as an attempt to understand the physics of VO 2 .
The other quantities that identify the variational spectrum are Z and ∆, shown in Fig. 6. We see that Z is decreasing with U but reaches a finite value Z = 1/4 for U → ∞. On the contrary, ∆ = 0 for U < U c , while ∆ = 0 for U > U c and increases monotonically to reach asymptotically the same value 1/4 for large U .
Further insights can be gained by the average values of the intra-dimer hopping and pairing, Eqs. (54) and (55), drawn in Fig. 7. As we mentioned, the intra-dimer hybridization 2δ * is monotonically decreasing with U , apart from the jump at the first order transition. More interestingly, around the transition the variational solution is characterized by a sizeable BCS order parameter ∆ * , which does not follow the behavior of the BCS coupling µ 1 present in the variational Hamiltonian. Indeed, while µ 1 is zero within our numerical accuracy for U < U c , yet a non negligible ∆ * develops just before the transition, see large for U > U c and increases monotonically with U , see Fig. 4, the actual order parameter ∆ * is appreciable only near the transition and fastly decreases with U to very tiny values. Also interesting is that, besides the intra-dimer superconducting order parameter, also an inter-dimer one arises, that can be for instance defined through where V is the number of sites. We find that ∆ * has actually the opposite sign of ∆ * , and both closely follow each other, rapidly decreasing with U , see Fig. 8. Even though ∆ * and ∆ * are everywhere finite for U > U c , suggestive of a superconducting phase that survives up to very large U , we believe that superconductivity, it it occurs at all, may appear only very close to the transition, where the value of the order parameter is larger. Indeed, the optimal solution with finite ∆ * and ∆ * and another solution in which both are forced to be zero are practically degenerate within our numerical accuracy for large U . Moreover, since the variational values of the order parameter are extremely small but close to U c , the inclusion of quantum fluctuations, for instance in the form of a Jastrow factor as in Eq. (2), would likely suppress superconductivity leading to a bona fide insulating wavefunction. Unfortunately, we can not compare this result with the DMFT analyses of Refs. 23,24, where superconductivity has not been looked for. 26

A. Large U limit
In order to appreciate qualities and also single out defects of the variational wavefunction, it is worth discussing the large-U solution. To leading order in 1/U , one can assume all P (n, Γ) = 0 but P (3) and P (4). The two constraints Eqs. (48) and (49) can be solved by defining with d ≪ 1, namely q → 1/2, while Eq. (50) is already satisfied since P (1) = P (2, +) = 0. Moreover, one readily recognizes that the variational solution asymptotically tends to acquire ψ ≃ β → π/4, θ → 0 and v 1 ≫ v 2 , v 3 , which is indeed what we find by numerical minimization. It then follows that This implies that µ 3 → 0, hence that µ 1 is determined by maximizing which leads to µ 1 = ǫ 2 /4d and where At leading order, the variational energy per dimer is therefore and v 2 ≃ 2t ⊥ /U , and takes the value In the case of a flat density of states and with t ⊥ = 0.5 D we recover the numerical result E var ≃ −7/6 D 2 /U , see Fig. 2. We note that, in spite of the hybridization δ * ∼ t ⊥ /U being small, the single-particle gap of the variational spectrum 2µ 1 ≃ U is large, as one should expect in a Mott insulator.
Coming back to the large-U value of the variational energy (69), one can readily realize that it coincides with the second order correction in t RR ′ to the energy of the state which is just a collection of singlets. In other words, in spite of being non-magnetic, our variational wavefunction is able to reproduce the correct super-exchange between dimers. This is a remarkable property that actually derives from the square-root dependence upon d of E hop , see Eq. (68). If we considered a more conventional Gutzwiller operator P commuting with the singleparticle density matrix, that amounts to further impose 2q = δ 2 * + ∆ 2 * , we would find E hop ∝ d, implying a transition into an unrealistic insulator with d = 0 above a critical U . The obvious defect of the wavefunction is that, since it emphasizes strongly the role of individual dimers, the hopping among dimers, although finite for any U , is under-estimated with respect to the intra-dimer one. Therefore we do not expect the wavefunction to be particularly accurate for small t ⊥ /D.

V. CONCLUSIONS
In this work we have proposed an extension of the Gutzwiller variational approach to account for correlated models which display metal-insulator transitions into Mott insulators that escape any simple single-particle descriptions, like the Hartree-Fock approximation. The wavefunction has still the same form as the conventional Gutzwiller wavefunction, with R identifying unit cells that may also be nonprimitive ones, with the novel feature that the operator P R is non-hermitean and does not commute with the local single-particle density matrix. In essence, this property realizes a variational implementation of a Schrieffer-Wolff transformation 27 , although only restricted to the lattice sites within each unit cell. We have shown that, by slightly reducing the variational freedom, this wavefunction, like the conventional Gutzwiller wavefunction, 14 allows for an extension of the Gutzwiller approximation to evaluate average values, approximation that becomes exact in the limit of infinite coordination lattices.
As an application, we have considered the Mott transition into a Peierls, or valence-bond, insulator, namely an insulator that is adiabatically connected to a collection of independent dimers. Such an insulator can not be described by Hartree-Fock simply because the singlet configuration of each dimer is not a Slater determinant. Specifically, we have considered the hypothetical situation shown in Fig. 9, where the splitting between the bonding and antibonding orbitals of each dimer is assumed not to be sufficient to lead, in the absence of in-teraction, to a band insulator. When interaction is taken into account, in the form of an on-site repulsion U , one expects, above a critical U , a transition from the metal into a Mott insulator. If magnetism is prevented, for instance by a sufficiently large splitting between bonding and anti-bonding orbitals and/or by frustration, the Mott insulator is non-magnetic. We have shown that our wavefunction overcomes the difficulties of Hartree-Fock theory and allows to study, albeit variationally, this transition. In particular we find that: (i) at the variational level the Mott transition is first order; (ii) the variational spectrum inside the Mott insulator looks similar to that of a Peierls insulator with a large hybridization gap, namely a large splitting between bonding and anti-bonding bands. In reality the gap is the Mott-Hubbard gap and the actual difference between the occupations of the bonding and anti-bonding bands, is small; (iii) inter-site singlet-superconductivity appears around the transition.
While (i) and (ii) are presumably true, as they have been also found by more rigorous calculations 23,25,28 , the emergence of superconductivity might be an artifact of the variational wavefunction. 26 Nevertheless, the possible occurrence of superconductivity is quite suggestive. It is known for instance that two-leg Hubbard ladders with nearest neighbor hopping display dominant superconducting fluctuations with the same symmetry that we find variationally 29 , although at half-filling they always describe non-magnetic spin-gaped insulators 30,31 because of nesting. Moreover, the uncorrelated wavefunction |φ is quite similar to the wavefunctions used in Refs. 32,33 to simulate t-J ladders. It would be surprising and interesting if this tendency towards superconductivity turned into a true symmetry breaking instability in higher dimensionality, as suggested by our analysis, which we think it is worth deserving further investigations. Note added: During the completion of this work, we became aware of a recent extension of slave-boson technique whose saddle-point solution closely resembles our variational approach. 34 Indeed the two conditions we impose on the Gutzwiller operator, Eqs. (4) and (5), are in oneto-one correspondence with the constraints identified in Ref. 34 within the slave-boson formalism.