The time-dependent Gutzwiller theory for multi-band Hubbard models

We formulate a multi-band generalisation of the time-dependent Gutzwiller theory. This approach allows for the calculation of general two-particle response functions, which are crucial for an understanding of various experiments in solid-state physics. As a first application, we study the momentum- and frequency-resolved magnetic susceptibility in a two-band Hubbard model. Like in the underlying ground-state approaches we find significant differences between the results of our method and those from a time-dependent Hartree-Fock approximation.


Introduction
The study of materials with medium to strong Coulomb-interaction effects has been a central subject for experimental and theoretical solid-state physics over many years. Despite enormous efforts and significant progress in some fields, however, our theoretical toolbox is still far from satisfactory for such systems. For quite some time, theoreticians in many-particle physics have focused on relatively simple model systems, such as the Heisenberg or the single-band Hubbard models. Only in the past ten years, attention shifted towards the study of more realistic models, e.g., multi-band Hubbard models. A very important impulse in that direction came from the limit of infinite spatial dimensions (D → ∞). The exact solution of Hubbard models in this limit leads to the Dynamical Mean Field Theory (DMFT), in which the original lattice model is mapped onto an effective single-impurity system that has to be solved numerically [1,2,3,4,5]. Although significant progress has been made in recent years in developing numerical techniques for the solution of the DMFT equations, it is still quite challenging and can be carried out only with limited accuracy. It is particularly difficult for the DMFT to study multi-orbital Hubbard models when the full (local) Coulomb and exchange interaction is included.
The theoretical interpretation of a number of experiments requires the study of two-particle response functions. For example, in magnetic neutron scattering the frequency-and momentum-resolved magnetic susceptibility is measured. The textbook method for the calculation of such response functions is the random-phase approximation (RPA), which can be interpreted as a time-dependent generalisation of the Hartree-Fock (HF) theory in the small amplitude limit, i.e., where the perturbation is considered to be sufficiently small. For electronic systems with medium or strong correlation effects, however, the ground-state description of a HF theory is well known to be often inaccurate. Therefore the RPA, as the time-dependent generalisation of the HF theory, is also a questionable approach for such systems.
A time-dependent Gutzwiller theory for the calculation of two-particle response functions was developed for single-band Hubbard models by Seibold et al. [31,32]. In recent years this approach has been applied with astonishing success to quite a number of such models and response functions [33,34,35,36,37,38,39,40,41,42,43]. It is the main purpose of the present work to generalise the time-dependent Gutzwiller theory for the investigation of multi-band models. A brief introduction into our method has already been given in Ref. [44]. All technical details, however, will be first presented here.
Our presentation is organised as follows: In chapters 2 and 3 we summarise the main results of the Gutzwiller variational theory for multi-band Hubbard models. In chapter 4 the reader will be reminded of the derivation which introduces the RPA as a timedependent generalisation of the HF theory. In a very similar way the time-dependent Gutzwiller theory ('Gutzwiller RPA') is introduced in chapter 5. The general Gutzwiller RPA equations are used in chapter 6 for the calculation of response functions for Hubbard-type lattice models. As a first application we study the magnetic susceptibility in a two-band model in chapter 7. A summary and conclusions close our presentation in chapter 8. The more technical parts of our derivation are referred to four appendices.

Multi-band Hubbard models and Gutzwiller wave functions
We study the general class of multi-band Hubbard modelŝ (2) contains all local terms, i.e., the two-particle Coulomb interactions (∼ U i ) and the orbital onsite-energies (∼ ǫ i ). We further introduce the eigenstates |Γ i ofĤ loc,i and the corresponding energies E loc Γ,i , i.e., H loc,i |Γ i = E loc Γ,i |Γ i .

Variational energy
As shown in Refs. [7,46], the expectation value of the Hamiltonian (1) with respect to the variational wave-function (4) can be evaluated in the limit of infinite spatial dimensions. We consider the expectation values of the local HamiltonianĤ loc , and the one-particle HamiltonianĤ 1 separately in sections 3.1 and 3.2. The additional constraints, which arise through the derivation in infinite dimensions are discussed in section 3.3. In section 3.4, we recall how the standard Gutzwiller energy functional for a single-band model is recovered from our general multi-band results.

Local energy
The expectation value of the local Hamiltonian (3) in infinite dimensions reads [48] where and To further evaluate the expectation value (11), we introduce the basis of Fock states (i.e., 'Slater determinants') in which certain spin-orbit states σ ∈ I are occupied. Mathematically, the indices I = (σ 1 , σ 2 , . . . , σ n ) are considered as ordered sets of spin-orbit states σ. Therefore we can use all standard set operations, such as I ∪ σ or I\σ. In addition, we define the number of orbitals in I as |I|. The states |I provide a basis of the local (atomic) Hilbert space. Hence, we can use them for an expansion of the eigenstates |Γ , and write the expectation value (11) as Finally, the uncorrelated expectation values of the transfer operators |I I ′ |, can be written as the determinant Here, Ω I,I ′ are the matrices in which the entries are the elements of the uncorrelated local density matrix that belong to the configurations I = (σ 1 , . . . , σ |I| ) and I ′ = (σ ′ 1 , . . . , σ ′ |I ′ | ). The matrix Ω J,J in (16) is defined as with σ i ∈ J ≡ (1, . . . , N)\(I ∪ I ′ ). In applications of the Gutzwiller theory to multi-band systems it would be quite cumbersome to evaluate the determinants (16) if the local density matrix (18) is nondiagonal. Fortunately, we are free to chose the local orbital basis in a way that suits us best. Therefore, we introduce an orbital basis, defined by (local) operatorsĥ With such a basis the expectation value (16) has the simple form Note that, for simplicity, we always use the same variable I for configuration states of the form (12) irrespective of the underlying orbital basis (e.g.,ĉ (21)). As will be shown in chapter 5, the time-dependent Gutzwiller theory requires to calculate the first and second derivatives of the energy with respect to all elements of the local density matrix, including the non-diagonal terms. Since the local density matrix enters the energy functional solely through matrices of the form (16) we only need to expand these matrices with respect to small perturbations up to second order in δC 0 γ,γ ′ around the diagonal ground-state matrix (20). This expansion is explicitly carried out in Appendix A.
In our derivation of the ground-state energy all local onsite energies were considered as part of the local Hamiltonian (2). For later use, however, we also need an expression for the expectation value of a general local one-particle Hamiltonian This expectation value is given as is the 'correlated' local density matrix.

Kinetic energy
The expectation value of a hopping term inĤ 0 , Eq. (1), is given as where we have introduced the (local) renormalisation matrix [48] The matrix H σ ′ I 1 ,I 4 contains three different contributions depending on whether the index σ ′ is an element of I 1 ∩I 4 , I 4 \(I 1 ∩I 4 ), or J = (1, . . . , N)\(I 1 ∪I 4 ). With the abbreviation The expectation value m 0;σ ′ I 1 \σ ′ ,I 4 in (28) has the same form as the one in (16), except that the index J has to be replaced by J\σ ′ . In case of a diagonal local density matrix one finds in agreement with results derived earlier [7]. Note that, in general, the renormalisation matrix is not Hermitian, i.e., it is Using (26), the expectation value of the one particle HamiltonianĤ 0 can be written as where we introduced the tensor

Physical constraints
As it turns out through the evaluation of expectation values in infinite dimensions, the variational parameters λ Γ,Γ ′ need to obey certain local constraints. These are Note that moving the operatorP †P relative toĉ † σ orĉ σ ′ in (34) would not alter the whole set of constraints. With the explicit form (5) of the correlation operatorP , the constraints read

Recovery of the 'standard' single-band energy functional
In case of a single-band model, the atomic eigenstates |Γ coincide with the configuration states |I . If we assume a translationally invariant ground state and the most general form of a local density matrix [48] where for those of the expectation values (15) which are finite. Here we used the notation As a consequence, the expectation value of the local Coulomb interaction in (6) reads For the single-band model with the correlation operator (7) and the local density matrix (37) the elements of the renormalisation matrix have the form Finally, the constraints in this case are given as As mentioned before, it is possible to overcome the complications that arise from a non-diagonal local density matrix by a simple transformation to a new orbital basis for which the local density matrix is diagonal by definition, In this new basis the constraints have the rather simple form where the non-diagonal constraints are automatically fulfilled by working with a diagonal correlation operator (λ 1,2 = 0). Equations (51)-(52) can be readily solved by introducing the expectation values which leaves us with only one variational parameter, the expectation valuem d for a double occupancy. The resulting renormalisation matrix is then diagonal and its elements have the well known form [49,50,51] Hence, the single-particle energy (31) is given as where we used the orthonormality relation In order to formally show the equivalence of both approaches we write the parameters λ Γ,Γ ′ in (7) as With these relations it is easy to show that the constraints (46)-(48) are indeed fulfilled. For example, the first constraint (46) can be written as If we use and the orthonormality relation (58) we readily find that (62) is indeed solved by the parameters (59)- (61). In the same way one can show the equivalence of the ground-state energy functionals.

The Time-Dependent Hartree-Fock Approximation
The approximation most frequently applied to two-particle Green's functions is the 'random-phase approximation' (RPA). This approach can be derived in various ways, e.g., by an equation of motion technique or in diagrammatic perturbation theory [52]. In this section, we use a different derivation which introduces the RPA as a time-dependent generalisation of the Hartree-Fock theory; see, e.g., Refs. [53,54]. If derived in this way, the approach can be generalised quite naturally in order to formulate a time-dependent Gutzwiller theory. This will be the subject of chapter 5.

The Hartree-Fock approximation
In the Hartree-Fock approximation a single-particle product wave function |Ψ 0 is used in order to investigate the ground-state properties of a many-particle system. Note that such wave functions are included in the Gutzwiller variational space by setting λ i;Γ,Γ ′ = δ Γ,Γ ′ . The expectation value of a many-particle Hamiltonian with respect to a Hartree-Fock wave function is a function of the single-particle density matrix. For example, for the Hamiltonian (1) it reads are the elements of the single-particle density matrixρ and is the expectation value of the Coulomb interaction in the Hamiltonian (2). Note that it will be more convenient both in the time-dependent Hartree-Fock and Gutzwiller theory, to use a different order of subscripts in the definition (65) of density matrices than, e.g., in section 3 or in previous work on the Gutzwiller theory.
To keep notations simple, we use the abbreviations υ ≡ (i, σ) for local single-particle states and Y = (υ, υ ′ ) for pairs of these indices. For example, the elements ofρ can then be written as With these new notations, the Hartree-Fock energy (64) reads where and for indices υ k = (i, σ k ) that belong to the same lattice site i. Further, we introduced the 'inverse' indexȲ ≡ (υ ′ , υ) for Y = (υ, υ ′ ). Note the symmetries which will be employed in the following section. The energy functional (68) has to be minimised with respect to all density matrices which belong to a single-particle product state. Such matrices are idempotent, i.e., they obey the matrix equatioñ If one imposes this constraint via a Lagrange parameter matrixη with elements η υ,υ ′ , the following equation has to be solved This condition leads tõ where we introduced the matrixh(ρ) with the elements Equation (74) is solved ifρ satisfies both (72) and [h(ρ),ρ] = 0 .
Starting with a certain density matrixρ we can introduce the 'Hartree-Fock' basis of states which diagonalise the Hamilton matrixh(ρ), i.e., Equation (76) is then solved by setting where the Fermi energy E F is determined by the total number of particles The density matrix (79) has to be reinserted into (68),(75) until self-consistency is reached. We denote the solution of these equations asρ 0 and introduce the corresponding Hamilton matrix

Equation of Motion for the Density Matrix
We consider two-particle Green's functions of the form where |Φ 0 is the exact ground state of our multi-band Hubbard Hamiltonian (1), and c ( †) υ (t) is the Heisenberg representation of the operatorsĉ ( †) υ with respect toĤ. As shown in most textbooks on many-particle physics, the Green's functions (82) naturally arise in 'linear-response theory' because they describe the time-dependent changes of the density matrixρ in the presence of a small time-dependent perturbation added toĤ [55,56,57]. After a Fourier transformation and using again the abbreviation and f Y (ω) and δρ Y (ω) defined accordingly. Ideally, we would like to calculate the time dependence of the density matrix where |Ψ(t) is the exact solution of the time-dependent Schrödinger equation for the HamiltonianĤ The expectation value (87) obeys the Heisenberg equation which contains the commutator In the time-dependent Hartree-Fock approximation, it is assumed that the solution |Ψ(t) of the Schrödinger equation at any time t is approximately given by a singleparticle product wave function. In this case, the expectation value of the commutator (90) can be evaluated by means of Wick's theorem. This leads to the equation of motion forρ(t), where the matrixh(ρ) has been introduced in (75). Equations (75) and (91) will be crucial also for our formulation of a time-dependent Gutzwiller theory in chapter 5.

Expansion for Weak Perturbations
We are only interested in cases wherê is a weak perturbation to the time-independent HamiltonianĤ. In this case, the density matrixρ(t) and the Hamilton matrixh(t) are given as where δρ(t) describes a 'small' time-dependent perturbation around the ground-state density matrixρ 0 , and With the expansion (93)-(94), the equation of motion (91) becomes These equations have to be solved for density matricesρ(t) that obey the matrix equation (72). After applying the expansion (93), Eq. (72) reads (to leading order in δρ(t))ρ Note that Eqs. (97),(99) just recover the time-independent Hartree-Fock equations derived in section 4.1.

Random -phase approximation (RPA) equations
Mathematically, the density matrix is a projector onto 'hole'-states,ρ h ≡ρ 0 . In addition, we define the projector onto 'particle'-states as With these two operators, we can decompose all matrices into their four components where v, w ∈ {p, h}. Note thath 0;vw has the elements An evaluation of the condition (100) for the components δρ vw (t) yields and δρ ww (t) = 0 (107) where v = w. Hence the components δρ pp (t) and δρ hh (t) can be neglected in the following compared to the leading fluctuations δρ hp (t) and δρ ph (t). We express the time-dependent quantities δρ vw (t) and δf vw (t) by their respective Fourier transforms δρ vw (ω) and δf vw (ω). The equation of motion (98) then leads to where the plus and minus signs correspond to vw = ph and vw = hp, respectively. With the abbreviation A = (α 1 , α 2 ) for pairs of indices α we find Here, the elements of the matrixŨ are given as The coefficients u υ,α in (110) have been introduced in Eq. (77) and determine the solutions |α of the Hartree-Fock equations. Equations (108) and (109) then yield with a matrixẼ defined as By comparing Eqs. (111) and (85) we find for the inverse of the two-particle Green's function Here we have added an increment iδ with δ = 0 + in order to ensure the correct boundary conditions of a retarded Green's function. ForŨ = 0, the inverse Green's function (113) which leads to Note thatΓ is not the exact Green's function for the single-particle HamiltonianĤ 0 since we just setŨ = 0 in (113), but kept finite the 'Hartree-Fock self-energy' contributions where, in the second line, we expanded [1 +ŨΓ(ω)] −1 into a power series with respect tõ UΓ. Both Eqs. (118),(119) are familiar expressions for the two-particle Green's function in the random-phase approximation.

Time-Dependent Gutzwiller Theory
The time-dependent Gutzwiller approximation has been first introduced for single-band Hubbard models by Seibold et al. [31,32]. In this section, we generalise this approach for the investigation of multi-band models. To this end, we set up an effective energy functional of the density matrix in section 5.1, which is used in sections 5.2-5.4 to derive the Gutzwiller RPA equations.

Effective Energy Functional
As summarised in chapter 3, the expectation value of the multi-band Hamiltonian (1) in the Gutzwiller theory is a function of the variational parameters λ Γ,Γ ′ and of the one-particle wave function |Ψ 0 . Like in the Hartree-Fock theory, the single-particle wave function |Ψ 0 enters the energy functional solely through the elements (65) of the non-interacting density matrixρ. It is therefore possible to consider the energy as a function of the density matrixρ and of the 'vector' of n p variational parameters λ Γ,Γ ′ (and λ * Γ,Γ ′ for Γ = Γ ′ ). The density matrix in the energy functional (120) must be derived from a single-particle wave function and, therefore, it has to obey the condition (72). Note that, in the following considerations, the density matrix will either be considered as a matrix (with respect to its two indices (i, σ) and (j, σ ′ )) or as a vector (with respect to its single index Y ). To distinguish both cases, we will denote the density matrixρ in some equations as ρ in order to indicate its vector interpretation.
The constraints (35)- (36) are also functions of λ and ρ and will be denoted as Here, n c is the (maximum) number of independent constraints, which, due to symmetries, is usually smaller than its maximum value N 2 so + 1, where N so is the number of spin-orbital states per lattice site. We assume that the functions (122) are real, i.e., in case of complex equations (33)-(34) their real and imaginary parts are treated separately.
By solving Eqs. (122) we can, at least in principle, express n c of the variational parameters (≡ λ d X ) through the density matrix ρ Y and the remaining 'independent' In this way, we obtain an energy functional which has to be minimised without constraints apart from Eq. (72) and the condition that the total particle number is conserved. For a fixed density matrixρ, the minimisation of (124) with respect to the determines these parameters as a function of ρ. This allows us to define the 'effective' energy functional which, for a fixed density matrix ρ, is given as the minimum of E GA with respect to λ i . With this effective functional we will formulate the time-dependent Gutzwiller theory in the following section. Using a Lagrange-parameter matrixη as in chapter 4.1, we find Here we introduced the matrixh(ρ) with the elements and used again the notationȲ ≡ (j, σ ′ ; i, σ) for Y = (i, σ; j, σ ′ ). The self-consistent solution of Eqs. (130)-(131) then yields the ground-state density matrixρ 0 , the matrix h 0 ≡h(ρ 0 ), and the corresponding single-particle 'Gutzwiller-Hamiltonian'

Gutzwiller RPA Equations
The derivation of RPA-type equations within the time-dependent Gutzwiller theory goes along the same lines as discussed in chapter 4 for the time-dependent Hartree-Fock theory. We add a small time-dependent field to our multi-band Hamiltonian (1). With the particular time dependence the expectation value of δV (t) reads where The renormalisation matrixq and the (correlated) local density matrixC c are defined in equations (27) and (25), respectively. With Eq. (127) they can both be considered as functions ofρ. The time-dependent field induces small fluctuations of the density matrix, Our main assumption is now that δρ Y (t) obeys the same equation of motion, as the density matrix in the time-dependent Hartree-Fock theory; see Eq. (98). Here, however, the Hamilton matrix is not derived from the Hartree-Fock functional (68), but from the effective energy functional (128), where the matrixK is given as The diagonalisation ofh 0 (or equivalently of the Gutzwiller Hamiltonianĥ 0 ) yields a basis |α with and a ground-state density matrix that is given as With the projectorsρ h ≡ρ 0 andρ p ≡ 1−ρ 0 , we define the particle and hole components of all matrices, as we did in Eqs. (102)-(104). The components δρ vw (t) of the densitymatrix fluctuations obey Eqs. (106)-(107), i.e., to leading order we can neglect δρ hh (t) and δρ pp (t). Hence, after a Fourier transformation we end up with the same form of RPA equations, for the two-particle Green's function matrix within the time-dependent Gutzwiller approximation.
One should keep in mind that the external 'fields' δf (ω) in (144) are 'renormalised', i.e., they are not the bare fields as they appear in (85), see Eq. (136). On the other hand, on the l.h.s. of (85) appears the 'correlated' expectation value of the density matrix, while in (144) we work with the fluctuations of the uncorrelated density matrix. Therefore, the 'true' Green's function seen in experiments may, in fact, be given as with certain frequency independent factors c Y,Y ′ . These factors, however, are of minor importance since they only affect the overall spectral weight and not the frequency dependence of the Green's function matrixG(ω). We can calculate them with the assumption that correlated and uncorrelated density-matrix fluctuations are related through the same renormalisation factors as the corresponding ground-state density matrices. For the one-band model it has been checked that this prescription is in fact the correct procedure for which the correlation functions fulfil the standard sum rules [34,37,43].

Second-Order Expansion of the Energy Functional
For an evaluation of the Gutzwiller RPA equations (144), we need to determine the matrixK which is given by the second derivatives (141) of the effective energy functional (128). To this end, we expand E GA up to second order around the ground state values ρ 0 and λ i;0 ≡ λ i ( ρ 0 ), Here, we introduced the matricesM ρρ ,M λρ ,M ρλ ,M λλ with the elements where the second derivatives on the r.h.s. are evaluated forρ =ρ 0 and λ i = λ i;0 . Note that there is no linear term ∼ λ i Z in (147) because of the minimisation condition (126). For our further evaluation, it is useful to write the second order terms in Eq. (147) in a more compact form by means of matrix-vector products, Here we used the symmetrỹ In the effective energy functional (128) the parameters λ i are determined by the minimisation condition (126). Applied to our second-order expansion (151) this condition yields which gives us the multiplet-amplitudes as a linear function of the densities δ ρ. This result leads to the quadratic expansion of the effective energy as a function of the density fluctuations δ ρ. In earlier work on the time-dependent Gutzwiller theory, Eqs. (153) and (154) have been denoted as the 'antiadiabaticity assumption'. In fact, these equations have the physical meaning that the local multiplet dynamics, described by fluctuations δλ i Z (t), are fast compared to those of the density-matrix fluctuations δρ Y (t). We will use the phrase 'antiadiabaticity assumption' in this work too although, strictly speaking, in our derivation it does not constitute an additional approximation.
With the functional (155), we could now proceed with our evaluation of the Gutzwiller RPA Eqs. (144). For practical applications, however, it is more convenient to determine the 'interaction kernel' (156) in a way that avoids the explicit solution of the constraint equations (122). This alternative procedure is the subject of the following section.

Lagrange-functional expansion
In the second-order expansion, described in section 5.3, we implemented the constraints (122) by explicitly eliminating a certain set of n c variational parameters. Although such a procedure can, at least in principle, always be applied, for the numerical implementation it is more convenient to impose the constraints by means of Lagrange parameters. To this end, we define the 'Lagrange functional' which depends on all variational parameters λ, the density matrixρ(= ρ) and the n c Lagrange parameters Λ n . The optimum variational parameters λ 0 Z , density-matrix elements ρ 0 Y , and Lagrange parameters Λ 0 n are then determined by the equations which have to be solved simultaneously.
We expand the Lagrange functional to leading order with respect to parameter (δλ Z , δΛ n ) and density fluctuations (δρ Y ). The second-order contribution has the form with matricesL ρρ ,L λρ ,L λλ defined as in Eqs. (148)-(150) only with E GA replaced by L.

The antiadiabaticity conditions
yield the n c equations and the n p equations Together these equations allow us to express the n p + n c parameter fluctuations δΛ n , δλ Z in terms of the density fluctuations δρ Y . These can be reinserted into (159) to obtain the desired quadratic functional solely of the density fluctuations, In Appendix B.1, we prove that the interaction matrixK Y,Y ′ in (164) is, in fact, identical to K Y,Y ′ in Eqs. (155)-(156).

Two particle response functions for lattice models
In the previous chapter we have developed the general formalism of the time-dependent Gutzwiller theory for the calculation of two-particle Green's functions. We will be more specific in this section and explain in detail how the response functions which are of interest in solid-state physics can be calculated within our approach.

Two-particle response functions
In solid-state physics one is usually not interested in the full two-particle Greens-functionG as it has been defined in (82). The properties, relevant for experiments, are certain linear combinations of elements ofG. For our translationally invariant model Hamiltonians (1) these are in particular the two-particle response functions or, more importantly, their Fourier transforms Here, we introduced the fermionic operatorŝ and the usual notation for the Fourier transform of a Green's function with arbitrary operatorsÔ,Ô ′ . With the abbreviation v = (σ, σ ′ ) for spin-orbit indices and the operatorŝ we can write (166) as The Green's functions (166) are still quite general since they include all possible channels of local coupling σ 1 ↔ σ 2 , σ 3 ↔ σ 4 . In experiments one usually measures response functions which are certain linear combinations, of some of the Green's functions (166), defined by the matrix κ v = κ σ,σ ′ . For example, the transversal spin-susceptibility χ( q, ω) is given as are the usual spin-flip operators and b is an index for the orbitals at each lattice site i. The spin susceptibility of a two-band Hubbard model will be investigated in chapter 7.

Response functions in the time-dependent Gutzwiller approximation
In order to apply the time-dependent Gutzwiller approximation, as developed in chapter 5, we have to expand the Lagrange functional (157) up to second order with respect to density-matrix (δρ) and variational-parameter fluctuations (δλ Γ,Γ ′ ). This means that we need an expansion of the constraints (35)-(36), of the local energies (9)- (11) and (24), and of the kinetic energy (31)- (32). The second-order expansion of the kinetic energy is more involved than that of the local energies and of the constraints.
In the latter there are only contributions from fluctuations at same lattice sites while in the kinetic energy local and non-local fluctuations (such as δ ĉ † i,σĉj,σ ′ Ψ 0 ) couple. Nevertheless, the calculation of the second-order Lagrange functional is tedious but otherwise straightforward. We therefore refer to Appendix C where the details of this derivation are presented. As shown in that Appendix, it is useful to introduce the operatorsB and to define the auxiliary Green's function matrixΠ( q, ω) with the elements We are actually interested only in the first 'element' of this matrix, i.e., the Green's functions (170) since they allow us to determine any response function of the form (171). As shown in Appendix D, however, the time-dependent Gutzwiller approximation leads to the following equation for the entire matrix (178) from which (170) can be extracted, Here,Ṽ q is the effective second-order interaction matrix, introduced in (C.40), and Π 0 ( q, ω) is the Green's function matrix (178) evaluated for the single-particle Gutzwiller Hamiltonian (132). As shown in Refs. [10,46], this Gutzwiller Hamiltonianĥ 0 ≡Ĥ eff 0 for our lattice Hamiltonian (1) has the form where the Lagrange parameters η σ 1 ,σ 2 are determined by the minimisation of the variational ground-state energy andǭ σ 1 ,σ 2 k is defined as The creation and annihilation operatorsĥ ( †) k,α of the effective single-particle Hamiltonian (180) can be written aŝ where the coefficients u σ,α are determined by a diagonalisation of (180). With these eigenstates the calculation ofΠ 0 ( q, ω) is now a simple task. For example, the first is the ground-state distribution function (143). In the same way, we can calculate all other elements ofΠ 0 ( q, ω) . The result is always the same as in (184) only with additional factors ∼ ǫ σ,σ ′ k or ∼ ǫ σ,σ ′ k+q due to the definition of the operators (176)-(177). For example, the second element in (178) leads to To summarise, with Eqs. (179), (184), (186), and the interaction matrix (C.40) we are now in the position to investigate any two-particle response function for our general class of multi-band models (1). As a first example, we study the magnetic susceptibility for a two-band model in the following section.

Magnetic susceptibility of a two-band Hubbard model
In this chapter we investigate the magnetic susceptibility of a two-band Hubbard model in three spatial dimensions. The model Hamiltonian and the Gutzwiller wave functions which we use for its investigation are introduced in section 7.1. In section 7.2 we discuss the Green's function matrices that we need to study in order to calculate the magnetic susceptibilities within the RPA and the Gutzwiller-RPA schemes. The numerical results for the two-band model are presented in section 7.3. Table 1. Two-particle eigenstates with symmetry specifications and energies.

Model and variational ground state
We investigate a Hubbard model with two degenerate e g orbitals per site on a cubic lattice. The local Hamiltonian (2) Here, e = 1, 2 labels the e g orbitals, s =↑, ↓ is the spin index and we use the convention ↑ ≡↓,↓ ≡↑. Due to the cubic symmetry the Coulomb parameters U, U ′ and the exchange parameter J are related to each other through Hence, only two of these three parameters can be chosen independently.
There are four spin-orbital states σ = (e, s) per atom, leading to a 2 4 = 16dimensional atomic Hilbert space. All eigenstates |Γ ofĤ 2b I with particle numbers N = 2 are simple Slater determinants of spin-orbital states |σ and their energies are The two-particle eigenstates are slightly more complicated because some of them are linear combinations of Slater determinants. We introduce the basis of two-particle states, which are used to set up the eigenstates ofĤ loc;i , see   For the variational ground state we can work with a wave function (4) that contains only diagonal parameters λ Γ,Γ . Non-diagonal parameters could only arise if we break the cubic symmetry or want to study states with magnetic orders not collinear to the chosen spin-quantisation axis. Note, however, that for the study of spin excitations we must allow for non-diagonal variational parameters, see below.
In our numerical analysis of this two-band model we will consider a tight-binding HamiltonianĤ 0 with generic hopping parameters which were already used in previous works and lead to the density of states at the Fermi energy shown in Fig. 1 (left). Due to the maximum in the density of states at approx. n σ = 0.29, in that range of band fillings there is the strongest tendency for a ferromagnetic state to be lower in energy than the paramagnet. This has already been demonstrated in Ref. [7]. Another important finding in that work is the huge importance of the exchange interaction J for the appearance of ferromagnetic order. This can be seen from the Gutzwiller phase diagram for our model in Fig. 1 (right). In contrast, the HF phase diagrams shows almost no dependence on the size of J; see also Ref. [44] where similar results have been reported for a two-band model in infinite dimensions.

The magnetic susceptibility
For the calculation of the spin susceptibility (172), we need to determine a Green's function matrix of the form (178) in which the operatorsÂ q v ,B q w ,B q w are given aŝ The matrix (178), which results from these operators is 4 + 16 + 16 = 36 dimensional. Due to symmetries, this dimension can be reduced to 20 for a general wave-vector q. Along symmetry lines the symmetry reduction could even go further. In our numerical calculations, however, we did not exploit such symmetry considerations since the numerical efforts for a two-band model are still moderate, even in three dimensions. Note that there is a difference between Hartree-Fock and Gutzwiller RPA calculations concerning the elements ofΠ 0 ( q, ω) which have to be taken into account in our calculation of the susceptibility In Hartree-Fock RPA, due to the locality of the interaction terms in Hubbard models and the symmetries of e g orbitals, the only elements ofΠ 0 ( q, ω) which contribute are , those which are diagonal with respect to the local orbital indices. This is different from the Gutzwiller RPA equations in which all Green's functions defined by the operators (193)-(195) have to be taken into account. In particular, Green's functions The reason for this difference is the non-locality of the interaction matrixṼ q in the time-dependent Gutzwiller theory.

Results
We prepare a ferromagnetic ground state in both HF and Gutzwiller approximation at band filling n σ ≈ 0.2987 in order to be close to the maximum of the DOS. In general both schemes will give different magnetisations for the same set of interaction parameters. Therefore, one could either perform the comparison for fixed parameters or fixed magnetisation (cf. also Ref. [59]). To avoid this inconsistency, we present results for interaction parameters which lead to a fully polarised ferromagnetic ground state in both approximations, i.e. m = 2n σ . Note that due to numerical reasons we have to stay slightly below this value in case of the Gutzwiller approximation. The corresponding interaction parameters are specified in the captions to Figs. 2, 3 which display the magnetic excitations obtained within both approximations.
These spectra are composed of a low-energy magnon part due to the breaking of spin-rotational invariance and a high energy Stoner continuum which reflects the particle-hole spin-flip excitations of the 'bare' system, i.eΠ 0 ( q, ω), cf. Eq. (179). For both methods we show the excitations along the (100) and (111) directions. The difference in these directions mainly arises due to the orientation dependence of the particle-hole dispersion which is significantly stronger along the diagonals.
One first important difference between HF and Gutzwiller approximation concerns the difference in the magnetic band splitting λ − = E ↑ F − E ↓ F . In HF theory this value is just given by λ − (HF ) = (U + J)m and thus for strongly correlated systems produces a large gap O(U) between the low energy magnon and a high energy Stoner continuum. On the other hand, we find λ − (GA) is significantly reduced with regard to its HF The magnon dispersion is fitted by ε magnon ( q)/t = D| q| 2 (1 + β| q| 2 ) with D 100 ≈ D 111 = 100 × 10 −3 . Scaling: t ≡ |t (1) ddσ | (c.f., Ref. [58]) and q x,y,z ∈ (−π, π).
counterpart. For the present system λ − (GA) ≈ 1/4λ − (HF ). Given the broadening of the Stoner continuum with increasing transferred momentum, the low energy magnon thus rapidly merges with the continuum in the time-dependent GA as can be seen from Fig. 3. As a consequence the excitation at ω ∼ J, corresponding to a respective spin-flip in the two orbitals, is only visible in HF+RPA along the (100) direction, whereas in the time-dependent Gutzwiller approach it is already within the continuum. Note that the overestimation of the Stoner excitation energy within HF+RPA is a longstanding problem in solid state theory as discussed in Ref. [60].
At q = 0 all the weight is contained in the zero frequency Goldstone mode. The existence of this excitation provides an important consistency check of the Gutzwiller+RPA approach similar to the analogous finding in HF+RPA. The positive dispersion of the magnon further demonstrates that the underlying Gutzwiller solution is a stable energy minimum which is not destroyed by the fluctuations. The spin-wave stiffness, i.e. the quadratric coefficient of the magnon dispersion, is significantly larger in HF+RPA than in time-dependent Gutzwiller theory. Note, however, that to a certain extend this huge difference is caused by an instability of the ferromagnetic ground state with respect to an incommensurate phase which is found for interaction parameters not much smaller than those used in Fig. 3.

Summary
In this paper we have given a detailed derivation of the time-dependent Gutzwiller approximation for multi-band Hubbard models. The basic assumptions which underlie the method can be summarised as follows. First, it is assumed that the dynamics of the Slater-determinant (upon which the Gutzwiller projector acts in the starting Ansatz) is determined by the so-called Gutzwiller Hamiltonian, Eq. (140), which leads to an equation of motion similar to standard RPA, Eq. (138). Second, the dynamics of the variational parameters is determined from the assumption that at each instant of time the energy is minimised. This leads to a linear relation between variational parameter and density fluctuations, Eq. (154). Third, as in the standard HF+RPA approach it is assumed that the external perturbation and thus the density fluctuations are small, Eqs. (106), (107). We have seen that also in the multi-band case these assumptions lead to a consistent theory in the sense, that an instability which is signalled within the Gutzwiller+RPA corresponds to a (second-order) phase transition which one would obtain from the bare variational Gutzwiller approximation. We have further demonstrated that for ferromagnetic ground states the Gutzwiller+RPA leads to the appearance of the Goldstone mode as expected for systems which break continuous spin symmetry.
The formalism as developed in its present form can now be straightforwardly applied to the investigation of correlation functions in strongly correlated multi-band systems as e.g. pnictides, manganites, cobaltates etc. On the other hand, a natural application of the theory would also comprise the investigation of e.g. orbital quenches for which the small amplitude assumption for the density matrices has to be abandoned. For single-band Hubbard models such a fully time-dependent formulation of the Gutzwiller approximation has been recently presented by Schiró and Fabrizio [61,62] where also the second assumption above has been replaced by separate equations of motion for the variational parameters. Future work should thus address the question whether their approach reduces to the present theory in the small amplitude limit and how it eventually can be extended to the multi-band case.

Appendix A. Second order expansion of determinants
According to Eq. (16), section 3.1, expectation values m 0 I,I ′ can be written as determinants of certain matrices A with elements that are linear functions of the local density matrix C 0 γ,γ ′ . In the variational ground state, C 0 γ,γ ′ is diagonal, and, if we chose a proper order of the orbitals γ, to 0-th order the matrix A is diagonal too, up to second order with respect to the matrix elements δa i,j . For this expansion one readily finds Note that for A 0 i,i = 0 the right-hand side is defined by the corresponding limit A 0 i,i → 0.

Appendix B. Invariance of Second-Order Expansions
Appendix B.1. Equivalence of the Lagrange-functional expansion In this section, we show that the interaction kernelK ρ,ρ Y,Y ′ in (164) is identical to K ρ,ρ Y,Y ′ in Eqs. (155)-(156). To this end, we choose again some arbitrary independent and dependent variational parameters λ i Z and λ d X , c.f. Eq. (123). By construction, the constraints (122) are automatically fulfilled as a function of λ i and ρ, i.e., we have Consequently, all first or higher-order derivatives of (B.1) with respect to λ i Z and ρ Y vanish. For example, the first-order derivatives lead to Using the matrices With the classification of dependent and independent variables we are in the position to evaluate the antiadiabaticity conditions (162)-(163). First, Eq. (162) leads to which, together with Eqs.
Here we introduced the six matrices Inserting this expression into the first row of Eqs. (B.10) and using we eventually find Equations (B.14), (B.12), and (B.9) now enable us to write all fluctuations δ λ i , δ λ d , and δ Λ as functions of the density fluctuations δ ρ. These relations can be inserted into the second-order expansion of the Lagrange functional, in order to calculateK ρρ Y,Y ′ in Eq. (164). However, to prove just the identity ofK ρρ (155) it is sufficient to apply only Eq. (B.9) to the expansion (B.15). This leads to As we will show below, the matrices (148)-(150) which determine the second order expansion (151) are the same as the corresponding matrices in (B.16). Hence, we have δE (2) = δL (2) . (B.17) Since the antiadiabaticity condition for δE (2) reproduces Eq. (B.14), the identity ofK ρρ Y,Y ′ and K ρρ Y,Y ′ is then finally demonstrated.
It remains to be shown that the matrices (148)-(150) agree with those in (B.16). To this end, we use the explicit form (124) of the energy functional (120) that appears in the definition of the matrices (148)-(150). As an example, we consider the matrix M ρρ and show that it is identical to the matrix in the first line of (B.16). With similar derivations one can prove the same for the other matrices (149),(150) and their counterparts in (B.16).
Using (124) and (148) we find Here, the matrices andg ρρ n , . . . ,g dd n are defined as in (B.11) only with L replaced by E or g n respectively. Obviously, the matrix in the first line of (B. 16 To prove (B.21), we use the fact that the second (total) derivatives of (B.1) with respect to the densities ρ Y vanish This equation, however, holds trivially, since (158) leads to for all parameters λ Z and in particular for λ Z = λ d X as it appears in (B.24).

Appendix B.2. Linear transformations of the density matrix
In investigations of our translationally invariant lattice systems (1) it turns out to be more convenient to work with fluctuations δ µ which are linearly related to the densitymatrix fluctuations, The antiadiabaticity condition for δ µ then reads Inserted into (B.28) this equation yields as claimed above.

Appendix C. Explicit form of the second-order expansion
We calculate the second-order expansion of the Lagrange functional with respect to the variational parameters λ i;Γ,Γ ′ and the density matrix (65). For the general consideration in section 5 and Appendix B it was convenient to subsume the parameters λ Γ,Γ ′ and their conjugates λ * Γ,Γ ′ in a set of n p parameters λ Z , c.f., Eq. (121). Here in this appendix, where we aim to resolve the explicit structure of the second-order expansion, it is better to take the difference between λ Γ,Γ ′ and λ * Γ,Γ ′ into account. The constraints (35)- (36), the local energy (24), and the renormalisation matrix (27) are all functions only of λ * i;Γ,Γ ′ , λ i;Γ,Γ ′ and of the local density matrix C 0 i;σ,σ ′ . For simplicity we use the joint variables A i v , (A i v ) * for all these local variables, i.e., it is either With respect to the parameters λ * i;Γ,Γ ′ , λ i;Γ,Γ ′ the second derivatives of (35)- (36), (24), and (27) are quadratic functions of the form ∼ (A i v ) * A i v ′ . Due to the Hermiticity of the density matrix the same can be achieved with respect to the local density matrix. Then the only finite second derivatives of the Lagrange functional whereas The second-order expansion of the constraints and the local energy is straightforward since only local fluctuations δA i v couple, δL where we introduced and the Fourier transforms of the local fluctuations All derivatives in this section (e.g., (C.7)) have to be evaluated for the ground-state values of the variational parameters λ i;Γ,Γ ′ , the density matrixρ, and the Lagrange parameters Λ i,n . Note that the density-matrix fluctuations δA q σ 2 ,σ 1 can be written as where the operatorÂ q v has been defined in (169). In addition to (C.6), we need to take into account the mixed terms ∼ δA i v δΛ i,n . In real space, their contribution is given as If we introduce the Fourier transforms δΛ q n of the fluctuations δΛ i,n , we can write (C.10) as Here, we used that the constraints g i,n are assumed to be real and lattice-site independent such that More involved than the calculation of (C.6) is the expansion of the kinetic energy. Here we find and The fact that the complex conjugates give the terms not explicitly shown in Eqs. (C.14)-(C.15) follows from the relations For our translationally invariant ground state it is more convenient to write (C.14)-(C.15) in momentum space: With the Fourier transforms of the local fluctuations the term (C.14) reads Here we assumed that the renormalisation matrix is lattice-site independent and introduced the tensor Note that for q = 0 the tensor (C.21), 23) has already been defined in (32). For the evaluation of the second ('transitive') term (C.15) we write the non-local density-matrix fluctuations as (C.24) Together with (C.8) this yields In principle, Eqs. (C.25)-(C.26) allow us to calculate all second-order couplings of density-matrix and parameter fluctuations that arise from δT (2) t . For numerical calculations, however, these equations are not very useful due to the explicit k dependence of (C.26). It is much easier to introduce the two auxiliary fluctuations is an abbreviation for quadruples of indices σ. With these definitions we can write (C.25) as Note that we introduced the two different fluctuations (C.27),(C.28) only because they allow us to write the second-order expansion in a relatively simple form. In fact, these fluctuations are not independent but related through Altogether we end up with the following second-order expansion of the Lagrange functional As described in section 5.4, the antiadiabaticity condition leads to an effective secondorder functional only of the density matrix. This condition can be evaluated directly for the second-order expansion (C.33) since the fluctuations δ A q , δ B q , δ B q are some linear functions of the density-matrix fluctuations δ ĉ † k,σ 1ĉ k+q,σ 2 , c.f., Appendix B.2. To this end, we distinguish the fluctuations of the local density matrix δ A q ρ and of the variational parameters δ A q λ as well as the corresponding blocks in the matrix (C.34), Note thatṼ q (unlikeK q ) includes finite couplings also between the fluctuations δ B q , δ B q . The calculation ofṼ q (for fixed q) only involves the handling of finite-dimensional matrices. In contrast, the evaluation of the functional (C.25) (instead of (C.29)) would have lead to significantly more complicated equations. 1 √ L s (n 0 k,α 2 − n 0 k+q,α 1 )δf (k+q,α 1 ),(k,α 2 ) .