Linear-Response Dynamics from the Time-Dependent Gutzwiller Approximation

Within a Lagrangian formalism we derive the time-dependent Gutzwiller approximation for general multi-band Hubbard models. Our approach explicitly incorporates the coupling between time-dependent variational parameters and a time-dependent density matrix from which we obtain dynamical correlation functions in the linear response regime. Our results are illustrated for the one-band model where we show that the interacting system can be mapped to an effective problem of fermionic quasiparticles coupled to"doublon"(double occupancy) bosonic fluctuations. The latter have an energy on the scale of the on-site Hubbard repulsion $U$ in the dilute limit but becomes soft at the Brinkman-Rice transition which is shown to be related to an emerging conservation law of doublon charge and the associated gauge invariance. Coupling with the boson mode produces structure in the charge response and we find that a similar structure appears in dynamical mean-field theory.


Introduction
Recent advances in ultra-fast spectroscopy allow us to monitor the dynamics of electrons on a femtosecond scale. This is especially interesting for strongly correlated materials, such as high-temperature superconductors [1,2,3], since in their case the spectroscopic probe is able to investigate the intra-electronic redistribution of excitation energies before the relaxation via the lattice starts. From the theoretical point of view, this is obviously a challenging problem since it requires a method capable of treating the relaxation dynamics of a strongly correlated system out of equilibrium. In this regard, a state-of-the art approach is the dynamical mean-field theory (DMFT) which has recently been applied [4] to the single-band Hubbard model in order to study the doubleoccupancy relaxation after laser excitation. However, for the application to real systems this method is rather demanding from a numerical point of view since it requires the self-consistent solution of complex single-impurity models driven out of equilibrium.
In this regard, the time-dependent Gutzwiller approximation (TDGA) is a promising alternative since it joins the numerical simplicity of standard random phase approximation (RPA) with the ability to capture important many-particle effects, as the Mott-Hubbard transition. In a series of papers [5,6,7,8], we have developed the TDGA which is based on a variational Ansatz for the Hubbard model [9,10] evaluated in the limit of infinite spatial dimensions [11]. This approach has recently been generalised for the study of multi-band Hubbard models [12,13] and is based on the expansion of the Gutzwiller energy functionals which depend on the density matrix and variational parameters related to the atomic eigenstates. In previous work [5,6] the latter have been eliminated by assuming that they instantaneously adjust to the density fluctuations ('antiadiabaticity approximation'). As a result, one obtains an energy functional which only depends on the density matrix and therefore the RPA for density-dependent forces [14,15] can be applied in order to evaluate response functions. This (approximate) version of the TDGA has been applied successfully to the evaluation of dynamical correlation functions in cuprate superconductors [16] including the optical conductivity [17] and the magnetic susceptibility [18,19]. It has also been related to Auger spectroscopy by calculating pair excitations in one- [8] and three-band [20] Hubbard models.
More recently the TDGA was extended by Schiró and Fabrizio [21,22,23] towards the inclusion of time-dependent variational parameters. Concerning the evaluation of response functions this approach can supersede the 'antiadiabaticity assumption', mentioned above, since the double-occupancy dynamics follows from a time-dependent variational principle. However, in Refs. [21,22] the authors focused on the study of quantum quenches for systems with a homogeneous ground state. In this case, the time dependence is captured by the double-occupancy dynamics whereas the single-particle density matrix is time-independent. Recent developments consider simultaneously the dynamics of the double occupancy and of the density matrix [24,25].
In this paper we will re-derive the TDGA for a time-dependent Gutzwiller variational wave-function applied to multi-band Hubbard models. Our resulting equations of motion will explicitly capture the coupling between the time-dependent variational parameters and the density matrix. We will analyse these equations in the small-amplitude (i.e., linear response) limit and apply the theory to the evaluation of dynamical charge correlations in the single-band case. It turns out that the previous formulation of the TDGA [5,6] is recovered in the low-frequency limit. However, the incorporation of fluctuations in the time-dependent density of double occupied states ("doublons") leads to additional spectral weight above the band-like excitations in very good agreement with exact diagonalisation and DMFT. The paper is organised as follows: In Sec. 2.1 we introduce the time-dependent variational principle which is underlying the present work. Since the corresponding expectation values are evaluated with multi-band Gutzwiller wave-functions the latter are presented in Sec. 2.2.
The evaluation of time-dependent matrix elements is performed in Sec. 2.3 which allows for the derivation of the Lagrangian and corresponding equations of motion in Sec. 2.4. Our investigations are specified for the single-band Hubbard model in Sec. 3, where we also discuss a two-site example which can be treated analytically. Finally, the small amplitude limit of the TDGA is derived in Sec. 4 and discussed in the context of response functions in Sec. 5. Numerical results for the dynamical charge susceptibility are presented in Sec. 6 and compared with dynamical mean-field theory (DMFT) and exact diagonalisation. We finally conclude our investigations in Sec. 7.

Variational Principle
The time-dependent Schrödinger equation for a general time-dependent Hamiltonian H(t) ( = 1), can be obtained by requesting that the action is stationary with respect to variations of the wave function. It is in general convenient to perform this variation based on a real Lagrangian [26] which leads to equations of motion that are independent of the phase and the norm of the wave-functions. If one restricts the wave-function |Ψ(z i (t)) to a certain trial form, depending on a set of (in general complex) 'variational parameters' z i (t), the Euler-Lagrange equations for z i (t) provide an approximation for the time evolution of the system. For example, restricting the wave-function to Slater determinants and using the amplitude of the single particle orbitals as variational parameters yields the time-dependent Hartree-Fock approximation. In this work, we will consider variational wave functions of the Gutzwiller form [9,10,11] for general multi-band Hubbard models leading to the timedependent Gutzwiller approximation.

Gutzwiller energy functionals for multi-band Hubbard models
We first recall some results of the conventional Gutzwiller approximation for the ground state properties of multi-band systems. We aim to study the physics of the following family of models, where t σ,σ ′ i,j denotes the 'hopping parameters' and the operatorsĉ i,σ annihilate (create) an electron with spin-orbital index σ on a lattice site i. The local Hamiltonian is determined by the orbital-dependent on-site energies ε i;σ 1 ,σ 2 and by the two-particle Coulomb interaction U σ 1 ,σ 2 ,σ 3 ,σ 4 i . Upon introducing the eigenstates |Γ i and eigenvalues E i;Γ of (4) (which can be readily calculated by means of standard numerical techniques) the local Hamiltonian can be written aŝ Multi-band Gutzwiller wave-functions have the form where |Ψ S is a normalised Slater determinant and the local Gutzwiller correlator is defined asP Here we introduced the variational-parameter matrix λ i;Γ,Γ ′ which allows us to optimise the occupation and the form of the eigenstates ofP i . The evaluation of expectations values with respect to the wave function (6) is a difficult many-particle problem, which cannot be solved in general. As shown in Refs. [27,28], one can derive analytical expressions for the variational ground-state energy in the limit of infinite spatial dimensions (D → ∞). Using this energy functional for the study of finite-dimensional systems is usually denoted as the 'Gutzwiller approximation'. This approach is the basis of most applications of Gutzwiller wave functions in studies of real materials and it will also be used in this work. One should keep in mind, however, that the Gutzwiller approximation has its limitations and the study of some phenomena requires an evaluation of expectation values in finite dimensions [29,30].
For the evaluation of Gutzwiller wave functions in infinite dimensions it is most convenient to impose the following (local) constraints [27,28] With these constraints, the expectation value of the local Hamiltonian (5) reads where can be calculated by means of Wick's theorem.
The expectation value of a hopping operator in infinite dimensions has the form where the (local) renormalisation matrix q σ ′ σ is an analytic function of the variational parameters and of the non-interacting local density matrix The explicit form of the renormalisation matrix is given, e.g., in Ref. [31] but will not be needed for our further considerations in this work. In the following, we assume that the correlated and the non-correlated local density matrix are equal, This is the case when all non-degenerate orbitals on a lattice site belong to different representations of its point symmetry group. In summary, the expectation value of the Hamiltonian (3), is a function of all variational parametersλ ( * ) ≡ {λ ( * ) i;Γ,Γ ′ } and of the non-interacting density matrixρ with the elements The same holds for the constraints (8), (9), for which we will use the abbreviation where n c is the maximum number of independent constraints. In Sec. 3 we apply these results to the special case of a single-band Hubbard model.

Evaluation of time-dependent matrix elements
In this section, we will apply the concept, introduced in Sec. 2.1, to the general class of Gutzwiller wave functions where the single-particle product states |Ψ S (t) and the local correlation operatorŝ are now time-dependent quantities. The state |Ψ S (t) can be written as Here, n γ ∈ (0, 1) determines which of the single particle states |γ(t) , described by the operatorsĥ † are occupied and u υ,γ (t) is a (time-dependent) unitary transformation, The functions u υ,γ constitute a second set of (time-dependent) variational parameters (in addition to λ i;Γ,Γ ′ (t)) and determine the single-particle wave function |Ψ S . Note that the time dependence of the operators (21) implies that the non-interacting density matrix is also time dependent. We start with a consideration of the time derivative in Eq. (2) which requires the evaluation of and its complex conjugate. With Eqs. (20) and (21), we find This equation allows us to evaluate the contribution of the first term on the r.h.s. of Eq. (24) as In the last line, we have used that, in all relevant applications,P G and |Ψ S have the same symmetry and, therefore, all contributions with γ = γ ′ vanish in (27). We now proceed with a consideration of the second term on the r.h.s. of Eq. (24). With the definition of the correlation operatorP G we find The r.h.s. of (29) can be evaluated by the standard diagrammatic techniques in infinite dimensions [27]. This leads to where m i;Γ 2 ,Γ 3 = m i;Γ 2 ,Γ 3 (t), as defined in Eq. (11), depends on the local elements of the density matrix (23).

Lagrangian and equations of motion
From Eqs. (27), (30), together with the expectation value of the Gutzwiller energy derived in Sec. 2.2, we are now in the position to derive the Lagrangian Eq. (2). However, we also need to include two sets of constraints, (i) the unitarity of u υ,γ and (ii) the Gutzwiller constraints (17). Therefore we finally obtain the following Lagrangian where Λ n (t) and Ω γ,γ ′ (t) are (real) Lagrange-parameters. As will be exemplified below in the single-band case, the original Hamiltonian can be time dependent which will reflect in a time dependence of E GA and allows for a coupling with arbitrary external fields.
From Eq. (31) the Euler-Lagrange equations can now be derived in the standard way. The equation for the variational parameters u * υ,γ reads which in terms of the density matrix (23) can be rewritten as Here we have introduced the 'Gutzwiller Hamiltonian' and a potential V λ which depends on the (time-dependent) phases of λ Γ,Γ ′ Note that the same equation of motion forρ is obtained in the previous formulation of the TDGA [5,6,7,8]. The new ingredient in the present formulation is the implementation of the explicit time dependence of the variational parameters λ Equations (33) and (36) forρ(t) and λ i;Γ,Γ ′ (t) constitute the time-dependent Gutzwiller theory for multi-band Hubbard models.

Evaluation of the time-dependent GA energy
In order to make the results, derived in the previous section, more transparent we consider the case of the single band Hubbard model, where t i,j denotes the 'hopping parameters' (t i,i ≡ 0) and the operatorsĉ ( †) i,σ annihilate (create) an electron with spin index σ on a lattice site i. We further introduced andn i,σ ≡ĉ † i,σĉ i,σ . All parameters in the Hamiltonian can be time and spatial-dependent which allows us to study the response to scalar fields v i (t), vector potential fields through the Peierls substitution, and modulations of the on-site interaction. Here, we introduced q = −|e| for electrons. The local Hamiltonian can be diagonalised by the states |Γ i = |d i , |σ i , |∅ i in which the site i is double occupied (i.e. has a doublon), single occupied by an electron with spin σ, or empty. Restricting for simplicity to paramagnetic states, we can work with a diagonal matrix of variational parameters, λ iΓΓ ′ = λ iΓ δ ΓΓ ′ . Thus the local 'Gutzwiller projection operator' reads, where the variational parameters λ i;Γ are related to the probability p i,Γ of finding a configuration Γ at site i according to We have four configuration probabilities per site which we denote as p i,Γ ≡ E i , S i,σ , D i for empty, single, and doublon occupied sites. In the present case, the constraints (8),(9) read, where ... i indicates that the index i is implicit everywhere in the expression. The first constraint is the statement Γ p Γ = 1, as it should be, while the second constraint implies that local charges are the same in the correlated and the uncorrelated state. Obviously this also guarantees that the total charge per spin in the system is the same in both states. We will show below that these constrains lead to equations of motions that nicely respect charge conservation.
According to Eq. (12) the expectation value of the one-band hopping operator in infinite dimensions can be written as where the (local) renormalisation factors are given by and we used the notation The q σ factors renormalise the probability amplitude for the annihilation of an electron on site j and the creation of an electron on site i. Each one of these processes has two possible channels. For example the creation of an electron at i with spin up can be seen as a transition from an empty state to an |↑ state [leading to the first term in Eq. (45)] or a transition from |↓ to a doublon state [leading to second term in Eq. (45)].
Since the variational parameters are now complex, these two channels can interfere either constructively or destructively affecting the total hopping amplitude. This issue is discussed further in Appendix A where the physical origin of this renormalisation is exemplified in a simple two-site case.
It is convenient to write the parameters λ Γ in terms of a real positive amplitude and a phase which are used in the following as the dynamical variables. With these definitions the hopping renormalisation factors read, with the definitions for site i, The phases ϕ i,σ , ϕ d,i have been eliminated in favor of χ i,σ and η i,σ . Note that ϕ ∅,i does not appear anywhere in the functional and therefore can be disregarded as a dynamical variable. In addition, we have used the constraints Eqs. (42),(43) to eliminate E and S σ in favor of D. Therefore our dynamical variables are the single-particle amplitudes u v,γ , the double occupancy D i and the phases χ i,σ and η i . In summary, the expectation value of the time-dependent single-band Hubbard model, is a function of the variational parameters and of the non-interacting density matrixρ,

Lagrangian and equations of motion
We are now in the position to evaluate the Lagrangian Eq. (31) which can be written as Note that, since we have implemented the constraints (17) explicitly, the corresponding Lagrange-parameter terms are not needed. The Lagrangian is invariant with respect to a gauge transformation of the form, . Notice that the hopping amplitude and the site energy transform in a way that generalises to the lattice the usual gauge transformation in the continuum (qA → qA − ∇χ, v → v +χ). Hence, χ i,σ plays the role of a gauge phase and implements charge conservation. Indeed, the Euler-Lagrange equations for χ i,σ yield the usual charge conservation law, where the current in a bond is given by The Euler-Lagrange equations for the variational parameters u * υ,γ yield again the equation of motion for the density matrix with the 'Gutzwiller Hamiltonian' matrix and the last term ensuring gauge invariance. In addition to (59), one obtains equations of motion for the double-occupancy parameters and the phases. For η i we obtain the Euler-Lagrange equatioṅ From Eq. (55) we see that η plays for D a similar role as the gauge phase χ for the charge. A time dependent η is equivalent to a change in the Hubbard U. However, there are important differences. In the case of a uniform system for n > 1 and U → ∞ the probability to find empty sites E → 0 which leads to q ∅,j,σ → 0 [c.f., Eq. (49)]. Then η becomes a gauge phase and D i is conserved as can be easily checked from the charge constraints. Using E instead of D as variational parameter one arrives to the conclusion that E is conserved for n < 1 and U → ∞. In general, however, η is not a gauge phase and D is of course not conserved. This reflects the fact that, when an electron jumps from a doubly-occupied site i to site j, one may have the process |d i |σ j → |σ i |d j which conserves D but one can also have the process |d i |∅ j → |σ i |σ j which does not conserve D. Therefore, in general, η i (t) is not arbitrary and has observable physical consequences as will be explained for the 2-site example (cf. the following section). On the other hand since χ i (t) is a gauge phase we can work in a gauge in which χ i (t) = 0. Finally form requiring stationarity with respect to D we geṫ

Non interacting limit
As a check of the consistency of the equations of motion it is instructive to see how the non interacting limit is recovered when U i → 0. In this case the double occupancy should factorise as D i = ρ i↑ ρ i↓ . Using this as an ansatz together with η i = 0 it is easy to check that Eq. (62) is satisfied. This follows from the fact that q i,∅,σ +q d,i,σ is the hopping renormalisation of the static theory and attains its maximum value q ∅,i,σ + q d,i,σ = 1 as a function of D i precisely when D i = ρ i↑ ρ i↓ . Thus its derivative as a function of D i evaluated at D i = ρ i↑ ρ i↓ vanishes (as can be checked from an explicit computation) and Eq. (62) is satisfied.
Using the same ansatz we notice that in this limit Eq. (61) can be written as, which completes the consistency check. On the last passage, we used Eqs. (50) and (57).
For small U, one would recover the time dependent Hartree-Fock approximation which in the small amplitude limit corresponds to the usual RPA.

Two-site example
In order to clarify the meaning of the TDGA equations it is interesting to consider the following two-site, two-electron example whose exact time-dependent evolution can be found analytically. The Hamiltonian is assumed to be time-independent with parameters v iσ = 0, t 1,2 = −t 0 ; the interaction on site 2 is infinite, U 2 = ∞, while U 1 is a variable. Albeit simple, the problem is in the strong coupling limit and provides a non-trivial test for the performance of the theory. Even more, been a zero dimensional problem we expect it to be the more demanding test bed for the TDGA which is based on infinite dimension results.

Exact solution:
The two-site Hamiltonian defined above is diagonalised by the states, and eigenvalues E ± given by, The time-dependent wave function can be expanded as such that the double occupancy is given by where Independently of the initial condition, as long as there is a finite overlap with both eigenstates, the double occupancy has a fluctuating part going like ∼ cos(ω 0 t).
As an example we choose the starting state at t = 0 as The probability of finding the system in the state |d∅ is then given by where the last approximate equality is valid for U >> t 0 . Clearly the probability to find the system in the state |s is 1 − D 1 from which one finds, with n i = n i↑ + n i↓ .

Time-dependent Gutzwiller approximation:
To solve the TDGA equations it is useful to note that U 2 = ∞ leads to D 2 (t) → 0. If the constraints where exactly satisfied in the GA this would also imply E 1 (t) = 0. However, the numerical solution of the static GA shows that E 1 is nonzero but very small for any U 1 (E 1 < 0.08) and vanishes for U 1 → −∞. In the following we assume for simplicity E 1 = 0. Then the constraint Eq. (42) implies D 1 = n 1 − 1 and one can evaluate the TDGA equations analytically. The hopping renormalisation factors simplify to, q 1σ = e −i(χ 1σ +η 1 ) q d,1σ , q 2σ = e −iχ 2σ q ∅,2σ , and the energy becomes The bonding single-particle state is defined bŷ Without loss of generality we set u 1 = n 1 2 e iφ and u 2 = n 2 2 which yields The Lagrangian reads, where we chose a gauge in which Furthermore, since φ and η 1 play the same role, we can set φ = 0. Thus we get a Lagrangian with two dynamical variables n 1 and η 1 which are conjugate. Their equations of motion have the form, The 'potential' v(n 1 ) is shown in Fig. 1. The static solution is given by η 1 = 0 and v ′ (n 1 ) = −U 1 /(4t 0 ). Hence, for U 1 = 0 the static charge is n 0 1 = 1 2 1 + √ 5 ≈ 1.62 and it decreases towards n 1 = 1 when U 1 → U c 1 = 4t 0 where a spurious Brinkman-Rice transition occurs. For negative U 1 instead the charge tends asymptotically to 2 when U 1 → −∞.
Equations (67),(68) can be readily solved for small oscillations around the static solution. We find that the oscillatory frequency is ω 0 = 4t 0 v ′′ (n 0 1 )(−v(n 0 1 )). For negative U 1 and until U 1 = 0 is approached we find that this frequency is in excellent agreement with the exact oscillation frequency of the two-site problem (c.f., Fig. 2). For positive U 1 as the spurious Brinkman-Rice point is approached, not surprisingly the approximations fail: the exact frequency increases monotonously as a function of U 1 while the approximate one vanishes at the Brinkman-Rice point. We attribute this to the failure of the GA in this finite-site system. Indeed we will show below that for high-dimensions a similar softening approaching a Mott state is real. The softening can be traced back to the fact that in the Brinkman-Rice phase the doublon charge gets frozen at zero therefore it is conserved and η 1 becomes a gauge phase which leads to an energy independent of η 1 . We will see that a similar phenomena is found in the usual Brinkman-Rice transition.

Linear response: The small amplitude limit
In the previous two-site model, the expansion of the effective potential v(n 1 ), Eq. (66), around the stationary value n 0 1 yields the dynamics in the vicinity of the GA saddlepoint. Such an expansion is also important for the evaluation of response functions, where a (weak) external perturbation drives the system out of equilibrium.
Based on our general scheme derived in Sec. 2.4, the small amplitude limit is obtained by expanding the equations of motion, Eqs. (33), (36), around the ground state valuesρ 0 and λ 0 Γ,Γ ′ , For example, the first term on the r.h.s of Eq. (36) becomes After the linearisation of Eqs. (33) and (36) and with the harmonic Ansatz we finally end up with a linear set of equations for λ i;Γ,Γ ′ (ω) and δρ υ,υ ′ (ω) which can be solved numerically. Note that at zero frequency, the l.h.s. of Eq. (36) vanishes. This equation then recovers exactly the 'antiadiabaticity assumption' which was used in the previous TDGA formulation [5,6,7,8]. Within the single band Hubbard model we will demonstrate in the following that the inclusion of the full time-dependence of the variational parameters λ i;Γ,Γ ′ (t) generates additional features in the dynamical charge correlations which are absent in the 'antiadiabatic approximation'.

Response functions in the single-band Hubbard model
In the single-band Hubbard model the three most relevant response channels are related to the coupling to time-dependent magnetic, charge and pair fields. This requires the computation of the (transversal) 'magnetic susceptibility' the 'charge susceptibility' and the 'pairing susceptibility' or their respective Fourier transforms where we have introducedn i = σn i,σ .

The magnetic and the pairing susceptibility
If the state |Ψ S is paramagnetic or is restricted to longitudinal magnetic order, the charge and (transverse) magnetic susceptibilities are decoupled. Moreover, if the ground state does not contain superconducting correlations, also the pairing fluctuations are decoupled from the magnetic and charge correlations. In this case, the (mixed) second derivative of the energy with respect to ĉ † i,↑ĉi,↓ and λ i;d vanishes. In a similar way, one can show that the second derivative with respect to pairing fluctuations ĉ † i,↑ĉ † i,↓ and δλ i;d vanishes. Therefore, in both cases, the linearised differential equations (33) and (36) are decoupled and the time dependence of δλ i;d (t) does not enter the computation of (transverse) magnetic and pairing correlation functions. The susceptibilities are then solely determined by the solution of Eq. (33) for the single-particle density matrix and the present approach agrees with the previous formulation of the TDGA involving the 'antiadiabaticity approximation' [5,6,7,8]. Therefore we concentrate in the following on the investigation of the dynamical charge-charge correlation function where the present approach is able to capture high-energy excitations on the scale of the Hubbard repulsion due to the explicit time dependence of the variational parameters.

Dynamical charge susceptibility
As derived in Sec. 3, the time dependence of the system is governed by small deviations of the density matrix δρ, the double-occupancy parameters δD i and the phase δη i from their saddle-point values (indicated by a '0' superscript). Note that we consider a GA ground state with η 0 i = 0 such that δη i = η i . The corresponding equations of motion have to be solved for small deviations from the GA saddle-point. For this purpose we expand the Gutzwiller energy functional Eq. (54) up to second order in the fluctuations around the saddle-point value E GA,0 . The first order contribution yields and the derivatives have to be taken at the saddle-point. Here, the first term describes particle-hole excitations within the 'bare' Gutzwiller Hamiltonian. The second term ∼ δD i vanishes due to the saddle-point condition, whereas the last term can be expressed in the small amplitude limit as where we have made use of Eq. (80). Thus in the small amplitude limit it is convenient to introduce a 'displaced double occupancy' such that the dynamics of η i and δD i can be expressed via the second order contribution of E GA as The following evaluation of E GA,2 is exemplified for a homogeneous paramagnet but can be straightforwardly generalised to arbitrary ground states. In momentum space the second order expansion involves fluctuations of η q and δD q which are coupled to fluctuations of the local δρ q = k,σ δρ σ k+q,k and transitive δT + q = k,σ ε 0 k+q + ε 0 k δρ σ k+q,k charge densities: where at the saddle point the renormalisation factors become site and spin independent (i.e. q 0 ∅,i,σ = q 0 ∅ , q 0 d,i,σ = q 0 d , q 0 iσ = q 0 ) and the primed letters denote derivatives which are specified in Appendix B. We have also defined the non-renormalised single-particle dispersion ε 0 k , whereas the GA quasi-particle dispersion will be denoted as ε k = (q i,σ ) 2 ε 0 k . Then from Eqs. (85),(86) the phase and double-occupancy dynamics is given by which upon Fourier transforming in the time domain yields for the double occupancỹ Here we have defined the renormalised couplings and the 'double occupancy' Green's function which has poles at Ω 2 q = K q /M. In case of a half-filled system one obtains where ε 0 is the energy of the non-interacting system, u ≡ U/|8ε 0 |, and κ q = 1 D D i=1 cos(q i ) . Interestingly the Brinkman-Rice transition u = 1 can therefore be associated with a 'soft mode' where Ω q=0 → 0. This softening is clearly associated to the extra gauge invariance condition that appears at the Brinkman-Rice point which makes the mass M to diverge.
Equations (87)-(96) show that in the present approach the interacting electron problem can be mapped to an effective electron-boson problem. Electron quasiparticles are coupled to a boson representing fluctuations of the double occupancy ("doublon fluctuations") or of its conjugate variable. Notice that the mass M of the fluctuation field diverges in the two situations discussed after Eq. (55) either E → 0 of D → 0 (c.f Eqs. (49),(61) and (93). This follows from the fact that in those cases η becomes a gauge phase and therefore gauge invariance requires that the last term of Eq. (87) vanishes.
In analogy with electron-phonon problems the dressed fluctuations are most conveniently evaluated by defining a (non-interacting) susceptibility matrix and an effective interaction matrix which is composed of the bare electron-electron interaction and the second order 'bosonic contribution' The susceptibility for the interacting system is then obtained from the following RPA series Note that in the static limit ω → 0 the matrix V ee q (0) is exactly the effective interaction obtained within the antiadiabaticity condition in Refs. [5,6].

Results
In this section, we present results for the local dynamical charge correlation function where |0 and |q denote ground and excited states of the single-band Hubbard model. We first study the low-density regime where we compare the TDGA spectra with exact results for the case of two particles. For higher densities we study the performance of the TDGA by comparing with DMFT and exact diagonalisation results.

Low-density regime
In the low-density limit n → 0 the energy of the non-interacting system is determined by ε 0 = −Bn where 2B would be the total bandwidth in case of a rectangular density of states. The saddle point double occupancy in this limit reads as which allows for the computation of the frequency of double-occupancy oscillations as Ω q = 2B + U, i.e., it is independent of momentum. Since the coupling to these fluctuations Eqs. (98),(99) vanishes with n we expect the local TDGA charge correlations to be composed of a renormalised low-energy part for 0 < ω < 2B and a high-energy excitation at ω = 2B + U with spectral weight proportional to the density. In case of two particles one can determine the eigenenergies in Eq. (104) from [32] where −4t ≤ ǫ k ≤ 4t denotes the single-particle dispersion with bandwidth 2B = 8t. The ground state is obtained for q = 0 and both particles at the bottom of the band, i.e., E 0 ≈ −2B. A particular solution in Eq. (106) is obtained for Q = q = (π, π) at E Q = U so that the excited state E Q − E 0 = 2B + U corresponds to our TDGA result discussed above. In addition, the exact solution involves two-particle excitations which are not present in the TDGA. The maximum excitation energy is obtained for q = 0 which can be estimated for a rectangular density of states as ω = 4B/(e 1/U − 1). The weight of these excitations in Eq. (106) vanishes for zero momentum transfer but clearly the exact solution displays high energy features in addition to the TDGA for small q due to the coupling between particle-hole and particle-particle excitations. Figure 3 displays the low-density local charge susceptibility Eq. (106), evaluated with TDGA and exact diagonalisation for U/t = 5 and U/t = 10, respectively. Results have been obtained for 2 particles on a 8 × 8 square lattice lattice (only nearest neighbor hopping −t).
As anticipated above, the spectra consist of the (dominant) low-energy bandlike particle-hole excitations in the range 0 < ω < 2B (cf. main panels) and a high-energy part at ω ≈ 2B + U which is resolved in the upper insets. Similar to the previous investigations, which were based on the 'antiadiabaticity assumption' [5,6], the TDGA gives a very good account of the low-energy part with respect to both, energy and spectral weight of the excitations. The new feature, which was previously missing [5,6] in the TDGA, is the high-energy part due to the explicit consideration of the doubleoccupancy time dependence. In order to estimate the associated spectral weight, we show in the lower insets to Fig. 3 the first moment which fulfills the sum rule Here T denotes the kinetic energy which in case of the TDGA has to be evaluated on the basis of the GA. From Fig. 3 it turns out that the onset of the high-energy excitations is accurately captured by the TDGA, however, it overestimates the associated spectral weight. On the other hand, this 'additional' weight is partially compensated for in the bandlike excitations such that the kinetic energy of the Gutzwiller approximation (i.e. M 1 (∞)) again gives a good account of the exact result.

Density dependence
We proceed by studying the doping dependence of χ loc (ω) as a function of doping which is shown in Fig. 4 for U/t = 8 for a square lattice. The spectra again separate into lowenergy band-like excitations and a high-energy part due to the double-occupancy time dependence. Starting form the low density limit the overall weight of the high-energy  excitations increases approaching half filling. In addition the high-energy feature shifts to lower frequencies upon doping as shown in the lower inset to Fig. (4).
We show also in the upper inset the imaginary part of the local double-occupancy propagator, i.e. ImD 0 loc (ω) = Im and D(q, ω) has been defined in Eq. (100). The double-occupancy excitations evolve from Ω = 2B + U (= 16t for the present parameters) in the limit n → 0 to the frequencies Ω q defined in Eq. (101) for the case of half-filling. In order to understand the doping dependence of the high-energy feature one has also to take into account the couplings Eqs.  Fig. (4) corresponds to the 'antiadiabaticity' result derived in Refs. [5,6]. With decreasing doping the increasing coupling between local density and double fluctuations increases and induces the shift of the high-energy feature to large frequencies.
In order to check the quality of the TDGA at larger doping vs. other approaches we compare our results with dynamical mean-field theory (DMFT) [33]. Despite freezing the spatial fluctuations beyond mean-field, DMFT takes fully into account the local quantum dynamics and it is in particular reliable to describe the evolution of the spectral weight as a function of temperature and other control parameters like doping. DMFT maps the lattice model onto an Anderson Impurity model, which is solved with an 'impurity solver', which in the present work is exact diagonalisation [34]. Within this method the bath of the AIM is discretised into N b levels, which here is taken to be 9. A Matsubara grid defined by an effective temperature β = 80/t is used and the stability of the results as a function of the two control parameters has been checked. Within DMFT, the local dynamical correlation functions can be obtained without further approximation. As a result of the discrete bath, the spectra appear more spiky than in the actual solution, but it has been shown that this approach describes accurately the evolution of the spectral weight (for instance of the optical conductivity) as a function of the various control parameters [35,36].   ) shows the local charge susceptibility obtained within DMFT for different fillings and U/t = 8. Despite the peaky structure (which hampers a direct comparison with the TDGA results) it is obvious that the main features correspond to those obtained within the TDGA, i.e., bandlike excitations on the scale of 8t and additional higher energy excitations which soften and gain spectral weight upon increasing density and approaching the insulating phase.
In order to test the performance of the TDGA we show in Fig. 6 (panels b,d) a comparison of the first moment M 1 (ω), Eq. (107), evaluated within TDGA (black solid lines) and DMFT (green circles) for U/t = 8 and densities n = 0.6 (panel b) and n = 1.0 (panel d). As anticipated above the TDGA gives a rather good account of the spectral weight evolution at lower densities where it is in good agreement with the DMFT data ( Fig. 6, panel b) given the uncertainties due to the finite number of bath states. However, due to the vanishing coupling between electrons and double-occupancy fluctuations at half-filling, all the TDGA spectral weight is contained in the band-like excitations in this limit so that the 'antiadiabaticity' result of Refs. [5,6] is recovered. Therefore the corresponding first moment increases much faster than DMFT but nevertheless both moments approach for ω → ∞ due to the agreement in the kinetic energies. One should note that at half-filling the GA ground state actually corresponds to a spin-density wave and also for such symmetry-broken states we find that at half-filling the antiadiabaticity result of Refs. [5,6] is valid. In fact, the TDGA charge excitations on top of such a ground state are in reasonable agreement with exact diagonalisation as demonstrated in Fig. 5 of Ref. [6].
For comparison, we also show in Fig. 6 the result of HF+RPA theory together with the spectra of our previous 'approximate' TDGA [5,6] where the double-occupancy fluctuations have been antiadiabatically eliminated. As discussed above, at half-filling the antiadiabatic approximation agrees with the 'exact' evaluation of the TDGA (cf. panels c,d in Fig. 6). The difference becomes pronounced at lower doping where the highenergy feature gets significant spectral weight and due to repulsion induces a softening of the bandlike excitations. From panel a of Fig. 6 it is clear that both effects are essential for reproducing the very good agreement of the first moment with the DMFT result at n = 0.6 whereas the approximate TDGA interpolates the spectral weight between the bandlike and high-energy excitations. Note, however, that for small frequencies ω → 0 the approximate result agrees with the 'exact' TDGA as has already been discussed in Sec. 4. In addition, the first moment of both approaches agrees in the limit ω → ∞ since it is set by the static GA kinetic energy.
In case of the conventional HF+RPA approach, where the local charge susceptibility is given by and χ 0 q (ω) has the same structure than the (11)-element of Eq. (102) but with nonrenormalised single-particle dispersions ε 0 k . Since HF+RPA obviously does not capture the double-occupancy dynamics, the local charge fluctuations originate from the bandtype excitations which are renormalised due to the RPA denominator and high frequency excitations are absent. Moreover, HF overestimates the kinetic energy so that the first moment of χ HF+RPA loc overshoots both the TDGA and DMFT result. As can be seen from Fig. 6 this discrepancy is most apparent close to half-filling where the renormalisation of the kinetic energy due to correlation effects is more pronounced. Upon reducing the band filling the low-energy part of the HF+RPA spectra approaches the TDGA result where, however, the latter approach additionally captures the high frequency excitations at 2B + U with spectral weight ∼ n.

Conclusions
In this paper we have developed the time-dependent Gutzwiller approximation for multiband Hubbard models. This approach is based on a time-dependent variational principle where expectation values are evaluated with the Gutzwiller variational wave-function in the limit of infinite dimensions. In contrast to the standard Gutzwiller approximation [9,10,11] both, variational parameters and the underlying Slater determinant, acquire a time dependence. In this regard our calculations generalise earlier investigations by Schiró and Fabrizio [21,22] who have studied quantum quenches in homogeneous systems where the time dependence of the density matrix does not couple to that of the variational parameters. On the other hand, momentum (or space) dependent out-of equilibrium displacements of the system require such a coupling as evident from our generalised equations of motion Eqs. (33), (36).
We have applied this theory in the small-amplitude, i.e., linear-response limit and exemplified for the case of dynamical charge correlations in the single-band Hubbard model. In an earlier formulation of the TDGA the so-called 'antiadiabaticity approximation' [5,6] has been applied, where the dynamics of the double-occupancy parameters was slaved by that of the density matrix. In contrast, the present approach explicitly incorporates the time dependence of the double-occupancy variational parameters which agrees with the previous formulation in the static limit. In addition it improves the theory in Refs. [5,6] by incorporating the high-energy features which are on the scale of the Hubbard repulsion for small densities and which position is in good agreement with that of exact diagonalisation. On the other hand, the spectral weight of the high-energy excitations is overestimated within the TDGA although it significantly improves the standard HF+RPA approach in this regard. Further refinement of the theory could be achieved by including the coupling between particle-hole and particleparticle excitations which have been studied in Ref. [8] in the framework of the GA.
It is interesting that in the present approach the Brinkman-Rice transition appears signaled by a collective mode whose frequency goes to zero. This is not due to the doublon fluctuation stiffness becoming soft but because the mass of the fluctuations diverges. We have shown that this divergence appears each time the double occupancy becomes a conserved quantity which is the case in the Brinkman-Rice case where D = 0. It remains to be seen which of these feature remain in an exact description although the similarity of the TDGA results with dynamical mean field theory (DMFT) suggest that at least in an approximate way this physics survives in real Mott transitions.
Within the DMFT it is quite diffult to study systems in which the momentum dependence of collective excitations is important as, for example, spin waves in insulators [37]. In such cases, the TDGA provides us with an important additional tool which complements the DMFT.