Gutzwiller variational wave function for multi-orbital Hubbard modelsin finite dimensions

Gutzwiller variational wave function for multi-orbital Hubbard models
in finite dimensions

Kevin zu Münster    Jörg Bünemann Fachbereich Physik, Philipps Universität Marburg, D-35032 Marburg, Germany
Abstract

We develop a diagrammatic method for the evaluation of general multi-band Gutzwiller wave functions in finite dimensions. Our approach provides a systematic improvement of the widely used Gutzwiller approximation. As a first application we investigate itinerant ferromagnetism and correlation-induced deformations of the Fermi surface for a two-band Hubbard model on a square lattice.

pacs:
71.10.-w,71.10.Fd,75.10.Lp,71.18.+y

I Introduction

Strongly-correlated electron systems display a variety of intriguing phases, such as superconductivity, (anti)ferromagnetism, or Mott insulating phases. In order to study the fundamental properties of strongly correlated lattice systems, simplifying Hubbard-type models are often employed. Unfortunately, the calculation of ground-state and dynamical properties is notoriously difficult even for these relatively simple models.

In one dimension, the density-matrix renormalization group (DMRG) method permits the numerical investigation of Hubbard-type models for fairly large chains. However, even modern variants of the DMRG such as tensor network approaches, Schollwöck (2011) are not satisfactory when applied to many-orbital models or higher-dimensional systems. In the limit of infinite dimensions, the dynamical mean-field theory (DMFT) Kotliar et al. (2006) maps the problem onto a single-impurity model whose spectral function must be calculated numerically. For multi-band systems, the solution of this task requires sophisticated quantum Monte-Carlo techniques and substantial computational resources. Concomitantly, it is very difficult to go beyond the mean-field limit and access multi-band models in finite dimensions.

In the absence of exact analytical or quasi-exact numerical methods, variational approaches have proven helpful. In this work, we employ the Gutzwiller wave function approach.Gutzwiller (1964) The evaluation of expectation values with the Gutzwiller correlated many-particle wave function poses itself a difficult many-body problem. Therefore, even the Gutzwiller-correlated single-band Fermi sea can be evaluated exactly only in one dimension. Metzner and Vollhardt (1987); Gebhard and Vollhardt (1987); Metzner and Vollhardt (1988); Gebhard and Vollhardt (1988); Metzner and Vollhardt (1989); Kollar and Vollhardt (2002) In the limit of infinite dimensions, the so-called Gutzwiller Approximation (GA) becomes exact for the single-band Hubbard model. Metzner and Vollhardt (1987); Gebhard (1990) Later, the method was extended to the multi-band case. Bünemann and Weber (1997, 1997); Bünemann et al. (1998, 2005) Recently, it has been combined with the density functional theory in a self-consistent manner to describe transition metals and their compounds. Ho et al. (2008); Deng et al. (2008, 2009); Wang et al. (2008, 2010); Weng et al. (2011); Yao et al. (2011); Tian et al. (2011); Lanatà et al. (2013); Dong et al. (2014); Schickling et al. (2014)

Despite many successes of the GA in improving our understanding of correlated metals, there are certain phenomena which it cannot describe properly. For example, in single-band models the Fermi surface is independent of the local Coulomb interaction within the GA, unless a state with broken spin or translational symmetry is considered. This is obviously incorrect, as can be seen already from straightforward perturbation theory for the paramagnetic ground state. Schweitzer and Czycholl (1991); Halboth and Metzner (1997) In order to describe a Fermi-surface deformation one needs to evaluate the Gutzwiller wave function in finite dimensions. A well established way to do this, is the ‘variational Monte Carlo method’ in which the Gutzwiller energy functional is minimized numerically on finite lattices. Tocchio et al. (2008, 2011, 2012) Although numerically less demanding than other techniques, such as quantum Monte-Carlo, this method still has significant finite-size limitations.

An alternative approach, which has first been proposed in Refs. [Gebhard, 1990,Bünemann et al., 2012], constitutes a systematic improvement of the GA for Gutzwiller wave functions on finite dimensional lattices. The method has been used successfully to study Fermi-surface deformations, d-wave superconductivity, and quasi-particle band structures in single band Hubbard models, Bünemann et al. (2012); Kaczmarczyk et al. (2013); Kaczmarczyk (2015a); Kaczmarczyk et al. (2015) t-J models, Kaczmarczyk et al. (2014) and periodic Anderson models. Wysokiński et al. (2015, )

Most transition metals and their compounds cannot be described properly by single-band models. For example, in iron-pnictides, such as LaOFeAs, all five -orbitals are partially occupied and may have to be taken into account in any model study that aims to describe the superconductivity or the antiferromagnetism in these systems.Schickling et al. (2011, 2012) Hence, it is clearly desirable to generalize the method, developed in Ref. [Bünemann et al., 2012] for the single-band model, to the multi-band case. It is the main purpose of this work, to formulate such a generalization. As a first application, we shall study ferromagnetism and Fermi-surface deformations in a two-orbital Hubbard model on a square lattice.

Our work is organized as follows. In Sect. II, we introduce the multi-band Hubbard model and the corresponding Gutzwiller wave function. Moreover, we present a detailed derivation of the diagrammatic expansion for ground-state expectation values. In Sect. III we discuss our results for ferromagnetism and interaction-induced deformations of the Fermi surface in a two-orbital Hubbard model on a square lattice. Finally, Sect. IV summarizes our findings and gives a brief outlook. For some technical details we refer to three appendices. A more detailed derivation of the results presented in this work can be found in Ref. [zu Münster, 2015].

Ii Gutzwiller wave functions

In this work, we employ Gutzwiller variational functions for multi-band Hubbard models. Such wave functions start from an independent-particle picture where the electrons are distributed over all lattice sites to optimize the single-particle energy. This statistical distribution leads to atomic configurations that are energetically unfavorable for finite Hubbard interactions. In the Gutzwiller wave function, the weight of such configurations is reduced with the help of the Gutzwiller ‘correlator’, a product of local operators, see below. Non-local (‘extended’) Gutzwiller correlators for the single-band Hubbard model can be studied analytically, Baeriswyl et al. (2009) and numerically using variational Monte Carlo. Tocchio et al. (2008, 2011, 2012) More recently, a two-dimensional bilayer Hubbard model was studied with the same method. Rüger et al. (2014)

ii.1 Multi-band Hubbard model

In this work, we investigate Gutzwiller-correlated wave functions for general multi-band Hubbard models. Only in the numerical applications in Sec. III, we shall be more specific by considering a two-orbital Hubbard model on a square lattice where the degenerate orbitals obey a - symmetry. The Hubbard Hamilton operator with purely local interactions reads

(1)
(2)
(3)

Here, and are fermionic creation and annihilation operator, respectively. The site index is given by and and the combined spin-orbital index by . The lattice indices run over all lattice sites of the lattice . Periodic boundary conditions apply. The hopping amplitudes and the coefficients of the on-site interaction energy are considered to be free model parameters.

The hopping and Coulomb parameters are restricted by symmetry. Spin conservation and rotational symmetry of the - orbitals reduce the nearest-neighbor and next-nearest-neighbor hopping amplitudes to four independent parameters. Furthermore, the coefficients of the on-site energy can be expressed solely in terms of the Hubbard interaction  and the Hund’s-rule coupling , as shown in appendix A. Note that the symmetry of the - orbitals is the same as that of the pair of - orbitals. Therefore, our two-band Hubbard model applies to - orbitals and to - orbitals equally well.

ii.2 Definition of Gutzwiller variational states

The Gutzwiller correlator is given by the product of the local Gutzwiller correlators for all sites on our lattice ,

(4)

If the context does not lead to any ambiguities, the local index will frequently be dropped in the following. In this work, we restrict ourselves to the homogeneous case where the variational parameters in are the same for all lattice sites. The local Gutzwiller operator is given by

(5)
(6)

with

(7)

where we already assumed that the parameters are real. The operators in (5), (6) that act on the site can be written explicitly as

(8)

where represents the identity operator on site . In our two-band application the local indices run over all local configurations which can contain up to four electrons.

In order to simplify the notation we define a product of local creation or annihilation operators by the introduction of the following symbols

(9)
(10)

where we introduced some arbitrary order of the spin-orbit indices . The multi-particle states

(11)

are uniquely determined by the lexicographical order of their sub-indices in (9).

A single particle product state (SPPS) can always be cast in the form

(12)

in some fermionic basis

(13)

We will assume that the SPPS are normalized, , and that the canonical commutation relations hold, . Now, we define the Gutzwiller wave function as

(14)

In the remaining part of this work, we optimize the Gutzwiller correlator and the SPPS so that the approximate ground state energy

(15)

becomes minimal.

ii.3 Diagrammatic expansion in finite dimensions

We consider the expectation value of some local operator

(16)

Note that the following calculation can equally be performed for expectation values of nonlocal operators such as .

As a first step, we follow the analysis for the single-band case derived in [Gebhard, 1990,Bünemann et al., 2012] and partly worked out for the multi-band case in infinite dimensions in [Bünemann, 2009]. In the numerator of Eq. (16) we pull the Gutzwiller correlators with indices to the right side of and denote the sandwich as ,

(17)

The operator and the squares of the Gutzwiller operator can be written in terms of creation and annihilation operators.

(18)
(19)

where the scalar contribution to could always be set to unity after rescaling the Gutzwiller wave function. However, we will postpone this step to a later stage of our analysis. In (19) we have introduced the number of electrons in a configuration state .

As a next step, we apply Wick’s theorem

(20)

where gives the sum over all possible contractions with respect to the density matrix with the elements

(21)

For example,

We depict the different contributions in a diagrammatic way, as shown in Fig. 1. Each summand of the operator and defines an ‘external node’ with weight or an ‘internal node’ with weights , respectively. Each operator contraction can be represented by a line which is either a ‘self-closing line’ (also denoted as ‘local contractions’ or ‘Hartree bubbles’) or connects two different nodes. In the following, we will simplify this diagrammatic analysis in three steps.

Figure 1: Diagrammatic representation of the numerator and denominator of the expectation value of an local operator on the lattice site . The blue square and the red circles gives the external and internal nodes, respectively. Black lines correspond to the single-particle density matrix.

First, we aim to eliminate all local contractions. Therefore, we map our operators to so called HF-operators, Gebhard (1990); Bünemann (2009) which do not have local contractions by definition. For example,

(23)

The mapping between the normal creation and annihilation operators and the HF-operators depends on the local density matrix as can be seen from the simplest case

An extension of this mapping for operator products is given in appendix B. We write the operator and the square of the Gutzwiller correlator as

(24)
(25)

where we set the coefficient . As mentioned above, this is equal to a rescaling of the Gutzwiller wave function by a factor which is always canceled out by the denominator in Eq. (16).

All operators in Eq. (20) are normal ordered because all site indices are different when we apply Wick’s theorem. We can set all local entries in to zero because we work with the HF-operators so that all local contractions vanish automatically. Therefore, we can carry out all contractions with a new density matrix

(26)

and drop the HF-operator notation at the same time

(27)

Without any nonzero local contraction we get

(28)

Thus, we can replace the fermionic operators , by Graßmann variables , , respectively. These Graßmann variables are nilpotent

(29)

In principle, the introduction of the HF mapping is not a necessary step for the introduction of Graßmann operators as we discuss in appendix C. All local entries of the new density matrix vanish so that the diagrammatic expansion cannot have nodes with self-closing lines, as shown in Fig. 2.

Figure 2: All nodes with internal lines cancel because the local entries of are set to zero.

The coefficients and are not affected by our mapping so that we can write

(30)
(31)

The numerator in Eq. (20) becomes

(32)

whereas the denominator reads

(33)

As a second step, we merge the diagrams of the numerator and denominator with the help of the linked cluster theorem. To this end, the lattice site restrictions on the right hand site of Eq. (32) must be removed. Therefore, we define

(34)
(35)

with

(36)
(37)

The exponential series expansion stays finite due to the nilpotency of the Graßmann variables. Therefore, the new coefficients and can be written as finite polynomials of the old coefficients and . The explicit expressions are given in appendix B. It is crucial that we perform the HF-mapping before we switch to the exponential form of our correlators.

These additional redefinitions allow us to cast Eq. (32) into the form

(38)

Note that the site index restriction disappeared. Eq. (33) can be rewritten as

(39)

Now, we are in the position to apply the linked cluster theorem (LCT), as described, e.g., in [Fetter and Walecka, 2003]. We find

(40)

where the summation is performed over all subsets of the lattice . The first few diagrams that are needed for the evaluation of the potential energy are shown in Fig. 3.

Some of the polynomials of in Eq. (38) vanish due to the nil-potency of the Graßmann variables. In contrast to the usual application of the LCT in many-body lattice theories, we can apply our expansion for a finite lattice as well. The nil-potency property allows us to add virtually as many nodes as we need to regroup our diagrams in all orders. Note that after the application of the LCT the nodes are contracted in such a way that all nodes have to be connected to the external nodes . This invalidates the nil-potency of the Graßmann variables inside the curly brackets. Therefore, several nodes can be located on the same site as long as these nodes are only connected indirectly. In a third step we will eliminate all internal nodes with two lines as described in the next section.

Figure 3: The first few connected diagrams that contribute to . The blue square represents the external node. The red circles represent the internal nodes. Black lines stand for the single-particle density matrix. The second and the forth diagram cancel out after the introduction of a gauge in the variational parameters as shown in subsection II.4.

ii.4 Limit of infinite dimensions

A scaling analysis of the ‘kinetic energy operator’ showsMetzner and Vollhardt (1989); Gebhard (1990, 1990) that the lines of the density matrix scale with the lattice dimension as

(41)

where we dropped the spin-band index for notational clarity, and gives the ‘one-norm’ (or ‘Manhattan metric’) of the displacement vector . One can show that all diagrams vanish if two nodes are connected to each other by at least three independent paths. This means, that there are at least three distinctive paths from node to node , so that none of the subsegments coincide. A trivial example is the diagram in which two internal nodes are connected by at least three independent lines. It will even be possible to eliminate all nontrivial diagrams if the nodes that have a single outgoing and incoming line are eliminated, as shown in Fig. 4. For this reason, a gauge in the variational parameters is introduced which sets the weight of these nodes to zero

(42)

This constraint must be incorporated in the optimization of the variational parameters . However, we can show that this constraint will not reduce the variational freedom in our model.zu Münster (2015)

Figure 4: The internal nodes with only two lines will be eliminated by a gauge in the variational parameters.

Then, the scaling of the hopping parameters

(43)

shows that all contributions with an internal node or two external nodes that are connected by three or more lines scale at least as . In the limit , the only remaining terms are given by

(44)
(45)

as shown in Fig.5.

Figure 5: Diagrammatic analysis in infinite dimensions.

In the rest of our work we will refer to the terms in Eq. (44), already derived in [Bünemann et al., 1998], as the ‘infinite- limit’. From these results we conclude that the constraints (42) ensure that the leading order terms of our diagrammatic expansion correspond to the exact Gutzwiller ground-state energy in infinite dimensions. As demonstrated for the single-band model in one dimension,Bünemann et al. (2012) the constraints are also essential for a rapid convergence of our expansion beyond the infinite- limit.

ii.5 Optimization of

In this work, we use the optimization algorithm which was introduced in [Bünemann et al., 2012]. The energy (15) depends on the variational parameters and the state , where the latter enters the functional only through the single-particle density matrix (21),

(46)

As shown, e.g., in appendix A of Ref. [Schickling et al., 2014], the minimization of with respect to leads to the following effective single-particle Hamiltonian

(47)
(48)

which has as a ground state. Hence, for the minimization of (46) we need to solve

(49)

and minimize with respect to ,

(50)

In order to solve these equations self-consistently, we usually start with the ground state of the free system, i.e., we set in (48),(49). Then we compute the density matrix, solve (50), and determine new parameters (48). The optimization terminates if the change of the effective hopping parameters between two cycles drops below some threshold. In order to test the stability of the algorithm, we can start from a different initial state. This initial state may be constructed from a perturbed kinetic energy operator . Usually the optimization algorithm remains stable against these perturbations but in some cases the fix-point of this map does not need to be unique, as shown in [Bünemann et al., 2012] where a symmetry breaking of the Fermi surface (Pomeranchuk phase) has been investigated.

Iii Results

iii.1 Magnetism

The occurrence of a ferromagnetic phase is favored by two conditions. The local Hamiltonian favors the formation of local magnetic moments for positive values of the Hund’s-rule coupling . Then, the two-particle eigenstates of the on-site energy with maximal local spin are lowest in energy, in accordance with Hund’s first rule. Therefore, for large values of , the ground state of the lattice system may show global ferromagnetism if the pre-formed local moments align. In contrast to that, the Stoner picture gives a different explanation for the origin of ferromagnetism. In this picture, a splitting between majority and minority bands reduces their mutual Coulomb repulsion due to the Pauli principle. This effect becomes significant when the density of states at the Fermi energy is large. In the Gutzwiller variational approach, the number of energetically costly multiple occupancies is reduced by an adjustment of the variational parameters. Therefore, we can expect that the Gutzwiller wave function predicts ferromagnetism at much larger interaction strengths than the uncorrelated SPPS. The Gutzwiller variational description leaves room both for the Stoner band splitting and the local moment formation as a source for itinerant ferromagnetism.

Figure 6: The density of states as a function of the density . The hopping amplitudes are .

Throughout this section, we focus on the kinetic energy operator with some generic amplitudes

(51)

where the coefficients and give the hopping amplitudes for a transition process between two orbitals to its horizontal and vertical neighbors, respectively. The coefficients and give the hopping amplitude for the transition process to the next-nearest neighbors with and without an inter-orbital transition. A detailed analysis of the lattice symmetries and all remaining hopping coefficients can be found in appendix A.

The density of states in Fig. 6 shows three peaks at , and . Below we investigate the ferromagnetic transition at , , , . The densities are located near the second peak in the density of states. The corresponding kinetic energies of the free system are , , and , respectively. Therefore, the average kinetic energy is very similar for all cases under investigation.

We define the quantity as

(52)

with , so that , and the total density remains constant. The total magnetization will be given by when both bands are considered. In order to obtain the optimal magnetization we will perform a scan in the - plane while we keep the ratio of and fixed to . Then, we optimize the band symmetric Gutzwiller variational parameters for each magnetization. Our diagrammatic expansion includes all lines with .

As seen in Fig. 7, in the Hartree-Fock approximation the magnetization jumps to a finite value at . The magnetization then increases monotonically until the ground-state is fully polarized at . A detailed analysis of the ground state energy shows that the nature of the jumps can be understood as a first-order phase transition. The Gutzwiller approach reveals a different picture. In second order, the magnetization shows a finite magnetization for and becomes fully magnetized for . This shows that the magnetization is shifted to much larger interaction strengths in the Gutzwiller wave function. The infinite- approximation becomes magnetized at and becomes fully magnetized for . Therefore, the second order diagrams in our diagrammatic expansion do not change the results on ferromagnetism significantly.

Figure 7: Magnetization of the Gutzwiller wave function. The black, red and blue crosses give the HF, the infinite- and the second order approximation, respectively.

Next, we analyze the density that lies very close to the second peak in the density of states, see Fig. 6. In the Stoner picture, the large density of states causes a finite magnetization at much smaller interaction strength. In our second-order Gutzwiller approach, the ground-state becomes already magnetized at although a precise evaluation of the threshold is hindered by numerical difficulties. The magnetization in the second-order approximation jumps to the fully magnetized state at . The infinite- approximation lies almost on top of the second-order expansion except at the transition to the fully magnetized state which occurs at . The HF-result shows the same qualitative behavior but the onset of ferromagnetism is at . Moreover, the magnetization increases more rapidly as function of the interaction strengths and saturates already at .

For the second-order magnetization result jumps to a finite value at and becomes fully spin polarized at . The transition points of the infinite- (HF) approximation lie at () and (), respectively. The magnetization curve shows the same qualitative behavior in all three approximations: The magnetization jumps to a small but finite value. Then the magnetization increases gradually as a function of whereby the slope is much steeper in HF than in Gutzwiller theory. Lastly, the magnetization jumps to full saturation at . In general, the critical values are much larger in Gutzwiller theory than in the Hartree-Fock approach. Note that the second-order terms to the result in lead to fairly small quantitative corrections. The magnetization onset requires a larger interaction strength for because the density of states is lower for than for . Furthermore, we can see that the transitions to the fully magnetized state occur at almost the same interaction strength as for in all approximations. This shows, that the transition to the fully polarized state depends on the density but not on the density of states.

For we still recover qualitatively the same behavior as for but the region between the onset of ferromagnetism and the transition to the fully polarized phase becomes smaller. Again, the critical values in Gutzwiller theory are about a factor of two larger than in Hartree-Fock theory.

In summary we can state that a large density of states at the Fermi energy promotes ferromagnetism. The Gutzwiller approach shows, however, that ferromagnetism, in general, requires large Coulomb interactions. Moreover, the Gutzwiller approach leaves room in parameter space for non-saturated ferromagnetism. For the system parameters, considered in this work, phases with long-range magnetic order are already well described within the GA. This supports the use of this approximation in many earlier works, see, e.g., Refs. [Schickling et al., 2014,Schickling et al., 2011,Schickling et al., 2012].

iii.2 Fermi surface deformations

In this subsection we show that the optimization of the SPPS can lead to a deformed Fermi surface. In some cases, the Fermi surface even changes its topology. Note that Fermi surface deformations within the Gutzwiller variational approach have already been studied in [Bünemann et al., 2012] for the single-band Hubbard model.

For our degenerate Hubbard model in the infinite- limit, neither for the Gutzwiller wavefunction nor within a more sophisticated DMFT calculation, we would find any correlation-induced changes of the Fermi-surface. Hence, all results in this section can be understood as an effect of the finite-dimensional evaluation provided by our higher-order expansion.

We examine the Fermi-surface deformations for the following parameter set,

(53)

The amplitudes are chosen in such a way that the topology of the Fermi surface changes near half-filling where the effect of the Gutzwiller correlator is strongest. The density of states is shown in Fig. 8. The peak near is caused by the change in the topology of the Fermi surface.

Figure 8: The density of states as a function of the density . The hopping amplitudes are .

Our diagrammatic expansion includes all lines with . In some cases, the optimization algorithm alternates between two fix points which are energetically very close. However, the Fermi surface of these fix points may differ significantly. In these cases it is useful to introduce some damping for the effective hopping parameters (48) in the self-consistency cycle. In our calculations, we usually find a fix-point after iteration steps.

In Fig. 9, the Fermi edges of the initial and optimized SPPS are shown for the densities , , , , and . The Hubbard/Hund parameters are set to and . Although the SPPS can have a small but finite magnetization in this parameter regime, we restrict ourselves to a paramagnetic wave function. The deformation of the inner Fermi surface between and start for densities . The outer Fermi edge between and is more robust. For densities , the optimized inner Fermi edge still has a closed topology while the initial Fermi surface topology is open. The optimization becomes difficult for densities larger than . Alternatively, a particle-hole symmetry can be used to determine the optimal Gutzwiller wave function.zu Münster (2015) In this way, we can show that the deviations in the Fermi surface are small for where the topology of the optimized Fermi surface becomes open.

The dependence of the Fermi surface deformations on the Hund’s-rule coupling is shown in Fig. 10. The Fermi edge for remains open for vanishing . An increase of the Hund’s-rule interaction strength to leads to the appearance of small islands in which both bands are filled. These islands collapse when we further increase the interaction strength to so that the Fermi surface becomes closed. The left panel of Fig. 12 shows that the energy gain increases linearly in and becomes vanishingly small for . From an energetic point of view, the transition from an open to a closed inner Fermi surface is gradual as a function of . The existence of intermediate islands also shows that the hopping matrix elements change gradually.

Figure 9: Fermi surface deformations for densities , , , , and . The local interaction strengths are set to and . The dashed lines give the initial Fermi edge and the solid lines give the optimized Fermi surface. Both bands are occupied in the region between the origin and the solid (dashed) black line. In the region between the black and the red lines only the lower band is occupied. For densities , the optimized inner Fermi edge has a closed topology while the initial Fermi surface topology is open.
Figure 10: Fermi surface deformations for different interaction strength . The density and the interaction strength are set to and , respectively. The deformations increase for larger values of . For , the Fermi surface topology (for ) is still open.
Figure 11: Fermi surface deformation for different interaction strength . The density and the Hund’s-rule coupling are set to and respectively. Small islands appear for in which .

The change in the Fermi-surface topology as a function of for and is shown in Fig. 11. For an interaction strength , the Fermi surface starts to deform from an open to a closed topology and small islands appear. The islands at the border of the Brillouin zone vanish for again. The energy gain increases linearly in as shown in the right panel of Fig. 12.

In this section, we showed that the interaction-induced Fermi surface deformations are clearly visible and, therefore, do not match the assumptions made in Fermi liquid theory that the Fermi surface of the non-interacting electrons is identical to the quasi-particle Fermi surface. Moreover, we showed that even the topology of the Fermi surface may change as a function of the Coulomb interaction-strength. The contributions beyond our second-order approximation still affect the Fermi surface and the density matrix. However, a higher order expansion of the single-band Hubbard model Bünemann et al. (2012); Kaczmarczyk (2015b) showed that the Fermi surface deformations are true features of the Gutzwiller wave function. Therefore, we can assume that the qualitative findings are valid in all orders of the approximation.

Figure 12: Left: The energy gain increases linearly in . Right: Energy gain as a function of .

Iv Conclusion and Outlook

In this work, we have given a comprehensive derivation of a diagrammatic method that allows us to evaluate multi-band Gutzwiller wave functions in finite dimensions. Our approach constitutes a systematic improvement of the widely used Gutzwiller approximation which corresponds to the zeroth order of our expansion.

As our first application, we studied the ferromagnetic phase transition in a two-band Hubbard model on a square lattice as a function of the model parameters for various band fillings. In general, a large density of states and a strong Hund’s-rule exchange favor ferromagnetism. In the Gutzwiller wave function, the ferromagnetic order is strongly suppressed so that much larger interaction strength are needed than predicted by the Hartree-Fock solution. Moreover, the regions in parameter space where non-saturated ferromagnetism occurs are much broader in Gutzwiller theory. As shown in earlier studies, this gives room for the experimental observations of non-saturated ferromagnetism, e.g., in transition metals such as nickel and iron. It turned out that long-range ferromagnetic order in our model is already well described within the Gutzwiller approximation.

As a second application, we investigated the interaction-induced deformation of the Fermi surface. These effects occur for large interaction strength, when the potential energy of the system is twice as large as the kinetic energy. For weaker interactions and small densities, the deformations of the Fermi surface can be neglected. Close to half band-filling and for special choices of the electron transfer parameters, the interactions can induce a change in the Fermi-surface topology from open to closed constant-energy contours. These effects are a result of the finite-order diagrams and cannot be seen in the Gutzwiller approximation.

It will be an interesting question for future work whether the Fermi surface deformation can lead to symmetry broken phases (Pomeranchuk phase) as in the single-band Hubbard model. In our two-band model, such a broken rotational spatial symmetry leads to different orbital densities which are energetically unfavorable. Hence, it is an open question if the ground state of our two-band model can have an asymmetric Fermi surface. Another open questions concerns the appearance of superconductivity in our model, as seen in earlier work on two-orbital Hubbard models.Zegrodnik and Spałek (2012); Zegrodnik et al. (2013, 2014)

Acknowledgements.
We like to thank Florian Gebhard for many valuable discussions in all stages of this project. Furthermore we thank Jan Kaczmarczyk for a valuable discussion of his higher-oder study of a single-band Hubbard model.

Appendix A Lattice symmetries

Fig. 13 shows the hopping processes to nearest and next-nearest neighbors. The amplitude for the transition from the () orbitals on site to the site equals the amplitude for the transition from the orbitals () on site to the site . The same holds for the amplitudes and for the transitions between the orbitals on and and the orbitals on sites and respectively. The amplitudes for the hopping processes from to between the orbitals is the same as between the orbitals. The symmetry of the orbitals does not allow any - transition to nearest neighbors. For transitions between next-nearest neighbors, the sign of the amplitudes will change after a rotation of so that . The symmetry of the inter-orbital hopping processes leads to a diagonal local density matrix . Furthermore, the rotational symmetry of the hopping processes within the same orbitals guarantees that all diagonal entries of the local density matrix are the same, .

Figure 13: Hopping amplitudes

The calculation of the on-site Coulomb interaction (3) requires the evaluation of two-particle expectation values of the Coulomb energy

(54)

with . These coefficients can be simplified after a decomposition of the - orbitals in terms of Laguerre and Legendre polynomials, respectively. The explicit derivation of the coefficients can be found in many text books (e.g. in [Sugano et al., 1970]) so that we simply state the result. The matrix representation of the two-particle sector of is given by

(55)

where the standard ordering , , , , , of the two-particle states has been used. Thus we can exclude any terms where only one electron switches from a to a orbital or that violates spin-conversation.

The variational coefficients in the Gutzwiller correlator (5) have the same structure (non-vanishing elements) as the matrix elements of the on-site Coulomb interaction in Eq. (55). In the paramagnetic case, the variational coefficients obey a spin-band symmetry. In the ferromagnetic case, all parameters are symmetric under an interchange of the band index.

Appendix B HF-operators and diagrammatic weights

In this section, we give explicit results for the mapping of an arbitrary operator product to its corresponding HF-operator. We use this mapping to define the coefficients , , , and that give the weights of the internal and external nodes in our diagrammatic expansion. A detailed derivation of this mappings, its inversion, and a derivation of all coefficients can be found in [zu Münster, 2015].

Consider a product of fermionic creation/annihilation operators . For the corresponding HF-operator , the evaluation of

(56)

shall, by definition, not include contractions between any pairs of operators . We use the notation

(57)

where denotes the number of internal contractions, e.g., for ,

(58)

Each internal contraction reduces the number of operators by two. With the abbreviation (57) we can give the following closed expression for the HF-operator of an arbitrary operator

(59)

where denotes the next smallest even number. This result agrees with the expressions for the definition of the HF-operators in [Bünemann, 2009].

In our diagrammatic expansion, we must express all operators in terms of HF-operators,

(60)

Let us contract both sides with an arbitrary operator . In order to determine the value of the coefficient , we have to evaluate the term where the operators and with , form the external contractions with . We need to shift the operators that are reserved for the external contractions to the front of the operator , and contract all remaining operators internally. The operator can be chosen without any restrictions and just indicates which operators are reserved for an external contraction. The contraction with can be carried out symbolically by a replacement of with the HF-operators . Thus, we get