Quantum fluctuations beyond the Gutzwiller approximation

TL;DR

This paper introduces a straightforward scheme to incorporate quantum fluctuation corrections into the Gutzwiller approximation for lattice Hamiltonians, improving the accuracy of response function calculations.

Contribution

The method is applicable to generic multi-band models without assumptions on the dynamics of variational parameters, extending the Gutzwiller approximation's capabilities.

Findings

Quantum fluctuations reproduce the magnetic susceptibility behavior seen in dynamical mean field theory.

The method accurately captures the smooth behavior of susceptibility across the Mott transition.

It recovers known results and provides new insights into quantum fluctuation effects.

Abstract

We present a simple scheme to evaluate linear response functions including quantum fluctuation corrections on top of the Gutzwiller approximation. The method is derived for a generic multi-band lattice Hamiltonian without any assumption about the dynamics of the variational correlation parameters that define the Gutzwiller wavefunction, and which thus behave as genuine dynamical degrees of freedom that add on those of the variational uncorrelated Slater determinant. We apply the method to the standard half-filled single-band Hubbard model. We are able to recover known results, but, as by-product, we also obtain few novel ones. In particular, we show that quantum fluctuations can reproduce almost quantitatively the behaviour of the uniform magnetic susceptibility uncovered by dynamical mean field theory, which, though enhanced by correlations, is found to be smooth across the paramagnetic Mott transition. By contrast, the simple Gutzwiller approximation predicts that susceptibility to diverge at the transition.

Equations

$$\mid \Psi ⟩=\prod _i\mathcal{P}(i)\mid {\Psi }_0⟩,$$

$$\mid \Psi (t)⟩=\prod _i\mathcal{P}(i,t)\mid {\Psi }_0(t)⟩,$$

$$\mathcal{W}(i,t)=\exp\bigg{(}i\sum_{\alpha\beta}\,K_{\alpha\beta}(i,t)\,c^{\dagger}_{i\alpha}\,c^{\phantom{\dagger}}_{i\beta}\bigg{)}\,,$$

$$d_{i\alpha }^{}=\mathcal{W}(i,t{)}^{†}{c}_{i\alpha }^{}\mathcal{W}(i,t)=\sum _{\beta }{U}_{\alpha \beta }(i,t)c_{i\beta }^{}.$$

$$\mathcal{P}(i,t)$$

$$\mid {\Psi }_0(t)⟩$$

$$\mathcal{P}(i,t)$$

$$P_{0\,\bar{n}\bar{m}}(i,t)=\big{\langle}\Psi_{0}(t)\Big{|}\,\mid\!\bar{m};i\rangle\langle\bar{n};i\!\mid\,\Big{|}\Psi_{0}(t)\big{\rangle},$$

$$\hat{\Phi }(i,t)\equiv \hat{\lambda }(i,t)\sqrt{{\hat{P}}_0(i,t)},$$

$$\displaystyle\mathrm{Tr}\Big{(}\hat{\Phi}(i,t)^{\dagger}\,\hat{\Phi}(i,t)\Big{)}$$

$$\displaystyle\mathrm{Tr}\Big{(}\hat{\Phi}(i,t)^{\dagger}\,\hat{\Phi}(i,t)\,\hat{c}^{\dagger}_{i\,\alpha}\hat{c}^{\phantom{\dagger}}_{i\,\beta}\Big{)}$$

$$=$$

$$⟨{\Psi }_0(t)\mid c_{i\gamma }^{†}\mathcal{P}(i,t{)}^{†}{c}_{i\alpha }^{}\mathcal{P}(i,t)\mid {\Psi }_0(t)⟩$$

$$=\sum _{\beta }{n}_{\gamma \beta }(i,t)R_{\alpha \beta }(i,t{)}^{},$$

$$\hat{R}(i,t)=\hat{Q}(i,t)\hat{S}(i,t),$$

$$Q_{\alpha\beta}(i,t)=\mathrm{Tr}\Big{(}\hat{\Phi}(i,t)^{\dagger}\,\hat{c}^{\phantom{\dagger}}_{i\alpha}\,\hat{\Phi}(i,t)\,\hat{c}^{\dagger}_{i\beta}\Big{)}\,,$$

$$4\hat{S}(i,t{)}^{-2}=1-\hat{\Delta }(i,t{)}^2$$

$$\Delta_{\alpha\beta}(i,t)=\mathrm{Tr}\bigg{(}\hat{\Phi}(i,t)^{\dagger}\,\hat{\Phi}(i,t)\,\Big{[}\,\hat{c}_{i\alpha}^{\phantom{\dagger}}\,,\,\hat{c}_{i\beta}^{\dagger}\,\Big{]}\bigg{)}\,.$$

$$\mathcal{P}(i,t{)}^{†}{\mathbf{c}}_i^{}\mathcal{P}(i,t)\to \hat{R}(i,t){\mathbf{c}}_i$$

$$\hat{R}(i,t)\to \hat{R}(i,t{)}^W=\hat{R}(i,t)\hat{U}(i,t{)}^{†},$$

$$\mathcal{W}(i,t)\mathcal{P}(i,t{)}^{†}{\mathbf{c}}_i^{}\mathcal{P}(i,t)\mathcal{W}(i,t{)}^{†}\to \hat{R}(i,t{)}^W{\mathbf{d}}_i.$$

$$\mathrm{Tr}\Big{(}\hat{\Phi}(i,t)^{\dagger}\,\hat{\Phi}(i,t)\,\hat{d}^{\dagger}_{i\,\alpha}\hat{d}^{\phantom{\dagger}}_{i\,\beta}\Big{)}=\delta_{\alpha\beta}\,n_{\alpha}(i,t)\,,$$

$$R_{\alpha \beta }(i,t{)}^W$$

$$\mathcal{H}=\sum _{i\ne j}{\mathbf{c}}_i^{†}{\hat{t}}_{ij}{\mathbf{c}}_j^{†}+\sum _i{\mathcal{H}}_i,$$

$$E(t)$$

$${\mathcal{H}}_{\ast }(t)$$

$$\mathcal{S}$$

$$i\frac{\partial \hat{\Phi }(i,t)}{\partial t}$$

$$i\mid {\dot{\Psi }}_0(t)⟩$$

$$\frac{\partial E_{\ast }(t)}{\partial \hat{\Phi }(i,t{)}^{†}}$$

$$\hat{\Phi }(i,t)+\delta \hat{\Phi }(i,t)$$

$$\mid {\Psi }_0(t)⟩+\mid \delta {\Psi }_0(t)⟩$$

$$\mathcal{K}(i,t)=\sum _{\alpha \beta }{K}_{\alpha \beta }(i,t)c_{i\alpha }^{†}{c}_{i\beta }^{},$$

Full text

We apply the method to the standard half-filled single-band Hubbard model. We are able to recover known results, but, as by-product, we also obtain few novel ones. In particular, we show that quantum fluctuations can reproduce almost quantitatively the behaviour of the uniform magnetic susceptibility uncovered by dynamical mean field theory, which, though enhanced by correlations, is found to be smooth across the paramagnetic Mott transition. By contrast, the simple Gutzwiller approximation predicts that susceptibility to diverge at the transition.

I Introduction

The Gutzwiller approximationGutzwiller (1963, 1965) is likely the simplest tool to deal with strong correlations in lattice models of interacting electrons. It consists in a recipe for approximate analytical expressions of expectation values in a class of wavefunctions, named Gutzwiller wavefunctions, of the form

[TABLE]

where $|\Psi_{0}\rangle$ is a variational Slater determinant, and $\mathcal{P}(i)$ a linear operator that acts on the local Hilbert space at site $i$ and depends on a set of variational parameters.

Curiously, the Gutzwiller approximation often provides physically more sound results than a direct evaluation of expectation values in wavefunctions like Eq. (1). For instance, the numerical optimisation on a finite-dimensional lattice of a variational Gutzwiller wavefunction for a single-band half-filled Hubbard model never stabilises a genuine Mott insulating phaseYokoyama and Shiba (1987); Capello et al. (2005), i.e. an insulator that does not break any symmetry, which intuitively is to be expected beyond a critical strength of the on-site repulsion. By contrast, the Gutzwiller approximation is instead able to describe such a genuine Mott transitionBrinkman and Rice (1970). The explanation of this strange outcome relies on the following observations. The first is that, in order to describe a genuine Mott insulator, one needs to add to the Gutzwiller wavefunction, Eq. (1), long range density-density Jastrow factorsCapello et al. (2005). However, the effect of such Jastrow factors disappears in lattices with coordination number $z\to\infty$, therefore, only in that limit, wavefunctions like Eq. (1) can faithfully describe Mott insulators. Moreover, right in that limit of $z\to\infty$, the Gutzwiller approximation provides the exact expression of expectation valuesMetzner and Vollhardt (1989); Bünemann et al. (1998). Therefore the Gutzwiller approximation should better be regarded as a recipe to evaluate approximate expectation values in Gutzwiller-Jastrow wavefunctions, which becomes exact when the coordination number tends to infinity, rather than in Gutzwiller-only wavefunctions. In other words, the Gutzwiller approximation applied on a lattice with finite $z$ is just the variational counterpart of dynamical mean field theory (DMFT)Georges et al. (1996) applied on that same lattice.

Recently, several attempts to include the Gutzwiller approximation inside DFT electronic structure codes have been performed with quite encouraging outcomesDeng et al. (2008); Ho et al. (2008); Deng et al. (2009); Wang et al. (2010); Tian et al. (2011); Yao et al. (2012); Schickling et al. (2012); Lanatà et al. (2013, 2014); Borghi et al. (2014); Tian et al. (2015); Lanatà et al. (2015). In this perspective, it might be useful to have at disposal a simple and flexible method to calculate linear response functions within the Gutzwiller approximation, in view of an extension of the so-called linear response TDDFTRunge and Gross (1984); Gross and Dobson (1996) to the case when DFT is combined with the Gutzwiller approximation.

There are already several works dealing with linear response in the Gutzwiller approximation, most of which limited to the single-band Hubbard modelVollhardt (1984); Seibold and Lorenzana (2001); Seibold et al. (2003, 2004); Schiró and Fabrizio (2011); B nemann et al. (2013). Extensions to multi-band models have been attempted Oelsen et al. (2011); von Oelsen et al. (2011), though under an assumption about the dynamics of the variational parameters that determine the linear operators $\mathcal{P}(i)$ in Eq. (1).

Here we shall instead present a very simple and general method to evaluate linear response functions within the Gutzwiller approximation without any preliminary assumption. The method is essentially an extension of the time-dependent Gutzwiller approximation of Ref. Schiró and Fabrizio, 2010 to a generic multi-band Hamiltonian, where the dynamics of the linear operators $\mathcal{P}(i)$ and of the Slater determinant $|\Psi_{0}\rangle$, see Eq. (1), are treated on equal footing. Linearisation of the equations of motion around the stationary solution, which is the equilibrium state, thus allows calculating linear response functions.

We note that the results of the Gutzwiller approximation at equilibrium coincide with the saddle point solution of the slave-boson theory in the path-integral formulationKotliar and Ruckenstein (1986), which, in multi-band models, corresponds to the so-called rotationally invariant slave boson formalism (RISB)Lechermann et al. (2007). Our present results in the linear response regime can therefore be considered equivalent to the quantum fluctuations corrections above the RISB saddle-point solution. We preferred here to derive such corrections to the action directly from the time-dependent Gutzwiller approximation rather than from the RISB theory, since the former is at least a well controlled variational scheme in lattices with infinite coordination number. However, both the notations as well as the language we shall use are actually closely related to RISB theory.

The paper is organised as follows. In Sec. II we briefly present the time-dependent Gutzwiller approximation, with some additional technical details postponed to the Appendix. In Sec. III we linearise the equations of motion around the stationary solution and derive an effective action for the fluctuations in the harmonic approximation. In Sec. IV we apply the method to the single-band half-filled Hubbard model, which allows a comparison with already existing results. Section V is devoted to concluding remarks.

II The Gutzwiller approximation in brief

Besides the original worksGutzwiller (1963, 1965) where M. Gutzwiller introduced a novel class of variational wavefunctions as well as an approximate scheme to compute expectation values, after him called Gutzwiller wavefunctions and approximation, and the subsequent demonstration that such an approximation becomes exact in the limit of infinite-coordination latticesMetzner and Vollhardt (1989); Bünemann et al. (1998), there are by now many articles where the Gutzwiller approximation is described in detail. Here we shall follow Ref. Fabrizio, 2013 and use its same notations.

The time-dependent Gutzwiller wavefunction is defined through Seibold and Lorenzana (2001); Schiró and Fabrizio (2010); Fabrizio (2013)

[TABLE]

which is the analogous of Eq. (1) where now $\mid\!\!\Psi_{0}(t)\rangle$ is a time-dependent variational Slater determinant, and $\mathcal{P}(i,t)$ linear operators on the local Hilbert space that depend on time-dependent variational parameters. For sake of simplicity, we shall not include in our analysis BCS wavefunctions nor operators $\mathcal{P}(i,t)$ that are charge non-conserving. The extension to those cases is simple, though notations get more involved.

Suppose that the Hamiltonian is written in terms of fermionic operators $c^{\phantom{\dagger}}_{i\alpha}$ and $c^{\dagger}_{i\alpha}$, $\alpha=1,\dots,2M$, that correspond to annihilating or creating a fermion at site $i$ in a chosen basis of Wannier functions $\phi_{i\,\alpha}(\mathbf{x},t)$, where $\alpha$ indicates both spin and orbital indices. Let us imagine a $U(2M)$ unitary transformation

[TABLE]

with $K_{\alpha\beta}(i,t)=K_{\beta\alpha}(i,t)^{*}$, which maps $c_{i\,\alpha}$ into a new basis set $d_{i\,\alpha}$ of single particle operators

[TABLE]

Evidently, if we consider the gauge transformation

[TABLE]

the Gutzwiller wavefunction $|\Psi(t)\rangle$ in (2) stays invariant and the transformed $|\Psi_{0}(t)\rangle$ remains a Slater determinant. Such gauge invariance, analogous to that of the RISB theoryLechermann et al. (2007), repeatedly appears in the calculations that follow.

The most general $\mathcal{P}(i,t)$ can be writtenFabrizio (2007, 2013) as

[TABLE]

where $n$ and $\bar{m}$ can be chosen to belong to the local basis of Fock states built with the operators $c^{\phantom{\dagger}}_{i\,\alpha}$. Alternatively, one can use a mixed-basis representation where $n$ labels Fock states in the original basis $c^{\phantom{\dagger}}_{i\,\alpha}$, and $\bar{m}$ Fock states in a different basisLanatà et al. (2008), e.g. the basis of the operators $d^{\phantom{\dagger}}_{i\alpha}$ in Eq. (4), which is also used to built the Slater determinant $|\Psi_{0}(t)\rangle$. We define the uncorrelated local probability distribution $\hat{P}_{0}(i,t)$, which is positive definite, by its matrix elements

[TABLE]

as well as the Gutzwiller variational matrix

[TABLE]

with matrix elements $\Phi_{n\bar{m}}(i,t)$. Expectation values of local and non-local operators in the Gutzwiller wavefunction (2) can be calculated explicitly in infinite coordination lattices if one imposes the following two constraints at any timeBünemann et al. (1998); Fabrizio (2013):

[TABLE]

where the fermionic operators within the spur must be regarded as their matrix representation in the local Fock space. The second constraint Eq. (11) plays the role of a gauge-fixing condition, exactly as in the RISB modelLechermann et al. (2007).

Another important ingredient is the wavefunction renormalisation matrix $\hat{R}(i,t)$ with elements $R_{\alpha\beta}(i,t)$, defined by solving the set of equations

[TABLE]

where the left hand side can be straightforwardly evaluated by the Wick’s theorem. As shown in the Appendix A, the solution of the above equation reads

[TABLE]

where $\hat{Q}(i,t)$ has matrix elements

[TABLE]

and the hermitian matrix $\hat{S}(i,t)$ is defined through

[TABLE]

where the matrix elements of $\hat{\Delta}(i,t)$ are

[TABLE]

The meaning of $\hat{R}(i,t)$ is that the action of the annihilation operator $c^{\phantom{\dagger}}_{i\alpha}$ on the Gutzwiller wavefunction is equivalent to the action of the operator

[TABLE]

on the Slater determinant $|\Psi_{0}(t)\rangle$, where $\mathbf{c}_{i}$ is a spinor with components $c_{i\alpha}$. One can readily show that under the gauge transformation Eq. (5),

[TABLE]

where $\hat{U}(i,t)$ has the matrix elements $U_{\alpha\beta}(i,t)$ of Eq. (4), so that Eq. (17) transforms into

[TABLE]

Since we have complete freedom in choosing $\mathcal{W}(i,t)$, a convenient choice is the unitary transformation that diagonalises the local single-particle density matrix, in which case the operators $d^{\phantom{\dagger}}_{i\alpha}$ are associated to the natural orbitals and satisfy

[TABLE]

while the matrix elements of $\hat{R}(i,t)^{W}$ acquire the simple expression

[TABLE]

The matrix $\hat{\Phi}(i,t)$ is in this case conveniently defined in the mixed-basis representation, where $n$ in $\Phi_{n\bar{m}}(i,t)$ refers to a Fock state in the original basis, and $\bar{m}$ to a Fock state in the natural one. Such a mixed-basis representation is useful since,throughout all calculations, one does not actually need to know what the natural basis is in terms of the original oneLanatà et al. (2008). Such a nice property is linked to the gauge-invariance, equations (5) and (6), of the theoryLechermann et al. (2007).

II.1 The model

We shall assume the generic Hamiltonian

[TABLE]

where $\mathcal{H}_{i}$ includes all on-site terms. If the constraints Eq. (10) and Eq. (11) are satisfied at any time $t$, then, in infinite coordination lattices, it holds thatBünemann et al. (1998); Fabrizio (2013)

[TABLE]

where

[TABLE]

may be interpreted as the Hamiltonian of the quasiparticles. Evidently, all expectation values can be straightforwardly evaluated since the uncorrelated wavefunction $|\Psi_{0}(t)\rangle$ allows using Wick’s theorem.

II.2 The action

In the time-domain the variational principle corresponds to searching for the saddle point of the actionSchiró and Fabrizio (2010)

[TABLE]

where the equivalence holds on provision that the constraints (10) and (11) are fulfilled at any time. The saddle point equations are readily obtained:

[TABLE]

where

[TABLE]

$\hat{T}(i,t)$ is a tensor with components $T_{nm;n^{\prime}m^{\prime}}(i,t)$, which is still functional of the matrices $\hat{\Phi}$ and $\hat{\Phi}^{\dagger}$ at site $i$ as well as at all sites connected to $i$ by the hopping. One can show that this tensor is hermitean, $\hat{T}(i,t)=\hat{T}(i,t)^{\dagger}$, which implies that the normalisation Eq. (10) is conserved by the time evolution.

II.3 Fate of the constraint

Concerning the second constraint, Eq. (11), we now prove that, if it is satisfied at the initial time, it will remain so at the saddle point solutions of Eq. (25) and Eq. (26). Suppose we have indeed found the saddle point $\hat{\Phi}(i,t)$ and $|\Psi_{0}(t)\rangle$. By definition, any small variation with respect to that solution must lead to a vanishing variation of the action. Let us consider the infinitesimal gauge transformation

[TABLE]

where the operator

[TABLE]

has infinitesimal matrix elements $K_{\alpha\beta}(i,t)=K_{\beta\alpha}(i,t)^{*}$, and $\hat{K}(i,t)$ is its matrix representation in the Fock space. We already mentioned that the energy $E(t)$ is gauge invariant so that the variation of the action, $\delta\mathcal{S}=\mathcal{S}^{\mathcal{W}}-\mathcal{S}$, simply reads

[TABLE]

Since $\hat{\Phi}(i,t)$ and $|\Psi_{0}(t)\rangle$ are solutions of the saddle point equations, it follows that $\delta S$ must strictly vanish for any choice of the infinitesimally small matrix elements $K_{\alpha\beta}(t)$, which implies

[TABLE]

thus just the desired result. It actually means that the term in parenthesis is conserved in the evolution. Therefore, if it is initially vanishing, it will remain so at any time, which thus implies that the constraint Eq. (11) is fulfilled during the whole time evolution.

II.4 Stationary problem

At equilibrium one needs to find the minimum of the energy with the two constraints Eqs. (10) and (11), which can be enforced e.g. by Lagrange multipliers, leading to the set of equations

[TABLE]

where $\Lambda(i)$ enforces Eq. (10), and the hermitean matrix $\hat{\mu}(i)$ with components $\mu_{\alpha\beta}(i)$ enforces Eq. (11). In whatever follows we shall assume to work in a mixed-basis representation where the operators $d^{\phantom{\dagger}}_{i\alpha}$ are associated to the natural orbitals, so that we must also ensure that

[TABLE]

The quasiparticle Hamiltonian in the natural basis, including explicitly the Lagrange multipliers, is therefore

[TABLE]

with $\hat{R}$ defined in Eq. (20). Working in the mixed-basis representation with the natural orbitals considerably simplifies all calculations.

Recalling that $\hat{T}(i)$ is still functional of $\hat{\Phi}$, Eq. (29) looks like a stationary non-linear Schrœdinger equation Lanatà et al. (2012, 2015). One can for instance solve it as in any Hartree-Fock calculation. Namely, one can find the eigenstates and eigenvalues of Eq. (29) assuming $\hat{T}(i)$ fixed, and impose that, when $\hat{T}(i)$ is calculated substituting the actual expression of the lowest energy solution $\hat{\Phi}_{0}(i)$, the two values coincide. The Lagrange multiplier $\hat{\mu}$ is fixed by imposing Eq. (11) and Eq. (19). In this way one finally gets the self-consistent $\hat{T}(i)$, which we shall hereafter denote as

[TABLE]

Once the latter is known, as well as the value of $\hat{\mu}$, one can also solve (29) for all eigenvectors, $\hat{\Phi}_{n}(i)$ and corresponding eigenvalues $E_{n}(i)$, with $E_{0}(i)=\Lambda(i)$. We shall denote $\mathcal{H}_{*}$, $\hat{R}$, $\hat{Q}$, $\hat{n}$ and $\hat{S}$ calculated with $\hat{\Phi}_{0}$ as $\mathcal{H}^{(0)}_{*}$, $\hat{R}^{(0)}$, $\hat{Q}^{(0)}$, $\hat{n}^{(0)}$ and $\hat{S}^{(0)}$, respectively, with the latter two matrices diagonal in the natural basis,

[TABLE]

We conclude by noting that the saddle point Hamiltonian Eq. (31) with the inclusion of the Lagrange multipliers is not anymore invariant under the most general $U(2M)$ gauge transformation, but only under a subgroup $G$ with generators $\hat{T}^{a}$ that commute with $\hat{\mu}$. This is common in theories where the gauge invariance implements constraints about physical states. In the natural basis representation, $\mu_{i,\alpha\beta}=\delta_{\alpha\beta}\,\mu_{i\alpha}$ is diagonal, so that the matrix elements of $\hat{T}^{a}$ must satisfy

[TABLE]

whose solution is straightforward. For any non-degenerate $\alpha$, i.e. such that $\mu_{i\alpha}\not=\mu_{i\beta}$, $\forall\,\beta\not=\alpha$, we associate the generators $T^{\alpha}_{i,\gamma\beta}=\delta_{\alpha\beta}\,\delta_{\gamma\beta}$ of $U(1)$ abelian groups. On the contrary, for any set of $\alpha_{i}$, $i=1,\dots,k$, such that $\mu_{i\alpha_{i}}=\mu_{i\alpha_{j}}\not=\mu_{i\beta}$, $\forall\,\beta\not=\alpha_{1},\dots,\alpha_{k}$, we can associate generators of a $U(k)$ Lie algebra.

III Fluctuations above the saddle point solution

Our goal is to determine the action of the fluctuations beyond the saddle point within the harmonic approximation. To that purpose we assume that

[TABLE]

where $\phi_{n}(i,t)$ for $n>0$ is regarded as a first order fluctuation, while, to enforce normalisation,

[TABLE]

In addition, the Slater determinant is defined through

[TABLE]

where $|\Psi_{0}(t)\rangle$ is properly normalised and includes the zeroth order $|\Psi^{(0)}_{0}\rangle$, solution of the saddle point, as well as a fluctuation correction $|\delta\Psi_{0}(t)\rangle$. The unitary operator

[TABLE]

where $\hat{\mu}(i)$ is the equilibrium Lagrange multiplier, and $\hat{W}(i,t)$ is the matrix representation of $\mathcal{W}(i,t)$.

Through the above definitions, the action becomes

[TABLE]

where $E_{*}(t)=\langle\Psi_{0}(t)\mid\mathcal{H}_{*}(t)\mid\Psi_{0}(t)\rangle$, being now

[TABLE]

and

[TABLE]

We expand $\mathcal{H}_{*}(t)$ up to second order in the fluctuations. The zeroth order is just $\mathcal{H}_{*}^{(0)}$. Since the stationary solution is the saddle point of the action, the expectation value of the first order expansion $\mathcal{H}^{(1)}(t)$ over the saddle point Slater determinant $|\Psi^{(0)}_{0}\rangle$ cancels with the first order expansion of the local energy $\sum_{i}\,\sum_{nm}\phi_{n}(i,t)^{*}\,V_{nm}(i)\,\phi_{m}(i,t)$. Therefore $\mathcal{H}^{(1)}(t)$ contributes to $E_{*}(t)$ with a second order term that, by linear response theory, reads

[TABLE]

where, hereafter, $\langle\dots\rangle_{0}$ will denote average over $|\Psi^{(0)}_{0}\rangle$, and the operators in Eq. (43) have an additional time dependence since are evolved with the saddle point Hamiltonian $\mathcal{H}_{*}^{(0)}$. The explicit expression of $\mathcal{H}^{(1)}(t)$ is

[TABLE]

where $\hat{R}^{(0)}(i)$ is the stationary value, while the explicit expression of the first order Taylor expansion $\hat{R}^{(1)}(i,t)$ is given in Appendix A.1, see Eq. (128).

There are several second order terms upon expanding $\mathcal{H}_{*}(t)$, which we shall consider separately. The first is simply

[TABLE]

whose expectation value over $|\Psi^{(0)}_{0}\rangle$ is an additional second order contribution

[TABLE]

which, together with $\delta_{1}E_{*}(t)$ in Eq. (43), endow the action with spatial correlations among the $\phi_{n}(i,t)$’s at different sites.

The next second order corrections to $\mathcal{H}_{*}(t)$ derive from the second order expansion of $\hat{R}(i,t)$

[TABLE]

where we distinguish two different contributions, see equations (130) and (131) in Appendix A.1. The reason of this distinction is that

[TABLE]

reproduces the bare excitation energy of the fluctuations. The last contribution to the energy of the fluctuations is therefore

[TABLE]

If we define new variables

[TABLE]

and the quadratic potential

[TABLE]

which has a retarded component $\delta_{1}E_{*}(t)$, see Eq. (43), the action of the fluctuations reads, upon defining $\omega_{n}=E_{n}-E_{0}$,

[TABLE]

which is just the action of coupled harmonic oscillators.

$\delta S$ in Eq. (III) can be for instance used to evaluate the fluctuation corrections to linear response functions of local operators. For any local observable $\hat{O}(i)$, let us define the matrix element

[TABLE]

Suppose we add a perturbation that couples to the local density matrix

[TABLE]

where the matrix $\hat{V}(i,t)$ with elements $V_{\alpha\beta}(i,t)$ represents the external field. Without loss of generality we can assume that the expectation value of $\delta\mathcal{H}(t)$ in Eq. (55) vanishes at the stationary solution. Since by assumption the external field is first order, the perturbation adds a second order correction to the action (III) that is

[TABLE]

In the presence of $V(t)$ the action transforms into that of forced harmonic oscillators, whose solution allows calculating the expectation value of any local operator $\hat{O}(i)$, see Eq. (54),

[TABLE]

at linear order in the external field.

III.1 Residual gauge invariance and would-be Goldstone modes

As we mentioned, the action Eq. (40), with the time dependent quasiparticle Hamiltonian defined in Eq. (41), is invariant under a subgroup $G$ of the initial $U(2M)$ gauge symmetry. This implies the existence of massless modes with singular propagators that diverge as $1/\omega^{2}$ at low frequency, which are the would-be Goldstone modes related to the fact that the saddle-point $\hat{\Phi}_{0}(i)$ is not invariant under $G$. Let us consider for instance a $U(1)$ subgroup of $G$ related to the non-degenerate state $\alpha$ in the natural basis. The associated adjoint charge is

[TABLE]

and its conjugate variable is readily found to be

[TABLE]

The role of $\varphi_{\alpha}(i,t)$ is just to enforce the constraint Eq. (11), i.e.

[TABLE]

Indeed we can always perform a gauge transformation on the fermions

[TABLE]

which makes $\varphi_{\alpha}(i,t)$ to disappear from the energy leaving just the time derivative term in the action,

[TABLE]

The condition of vanishing derivative with respect to $\varphi_{\alpha}(i,t)$ is therefore just the condition that the constraint is conserved.

It follows that we can always drop from the action all terms that contain the variables conjugate to the adjoint charges associated with the gauge symmetry $G$, on provision that, wherever $n_{\alpha}(i,t)$ appears, we replace it with $\langle\;c^{\dagger}_{i\alpha}c^{\phantom{\dagger}}_{i\alpha}\;\rangle_{t}$.

However, the above procedure does not involve all the coefficients $\phi_{n}(i,t)$; some of their linear combinations are untouched by gauge-fixing and remain genuine independent dynamical degrees of freedomJolicoeur and Le Guillou (1991). This fact, rather than being a limitation, it endows the theory with a richer dynamics.

IV Application to the half-filled Hubbard model

We now apply the above formalism to the simple case of a single band Hubbard model at half-filling, where all calculations can be worked out analytically and which also allows for a direct comparison with previous worksVollhardt (1984); Jolicoeur and Le Guillou (1991); Raimondi and Castellani (1993); Seibold and Lorenzana (2001); Seibold et al. (2003, 2004); Schiró and Fabrizio (2011); Sandri et al. (2012); B nemann et al. (2013). We will show that we can indeed recover known results, but also find few novel ones.

The Hamiltonian is in this case

[TABLE]

where $<\!ij\!>$ means nearest neighbour bonds on a $d$-dimensional hyper cubic lattice, and $z=2d$ is the lattice coordination number that must be sent to $+\infty$ for the calculation to be really variational.

The local basis comprises four states which we choose to be, in order, the empty configuration, $\mid\!\!0\rangle$, the doubly occupied one, $\mid\!\!2\rangle$, the singly occupied by a spin up electron, $\mid\!\,\uparrow\rangle$, and that occupied by a spin down one, $\mid\!\,\downarrow\rangle$. The most general charge-conserving $\hat{\Phi}$ has the following form, dropping for the meanwhile the site index,

[TABLE]

where the charge component, i.e. the matrix elements in the subspace \big{(}\mid\!0\rangle,\mid\!2\rangle\big{)}, is

[TABLE]

with $\sigma_{0}$ the $2\times 2$ identity matrix, and $\sigma_{i}$, $i=1,\dots,3$ the Pauli matrices, whereas the spin component, namely the matrix elements in the subspace \big{(}\mid\!\,\uparrow\rangle,\mid\!\,\downarrow\rangle\big{)}, is instead

[TABLE]

which allows a full spin-$SU(2)$ invariant analysisLi et al. (1989); Seibold et al. (2004). Normalisation implies that

[TABLE]

One can readily verify that the matrix $\hat{Q}$ with components

[TABLE]

can be written as

[TABLE]

where

[TABLE]

with $i=1,\dots,3$. Seemingly,

[TABLE]

IV.1 Stationary solution

As common when discussing the Mott transition in the single band Hubbard model, we shall be interested in the stationary solution within the paramagnetic sector, i.e. neglecting spontaneous breakdown of spin $SU(2)$ symmetry. Such solution at half-filling is characterised by a site independent

[TABLE]

with

[TABLE]

Under this assumption

[TABLE]

so that the quasiparticle Hamiltonian is just a tight-binding model with renormalised hopping, i.e.

[TABLE]

and natural and original orbitals coincide. It follows that the stationary Slater determinant is the non-interacting Fermi sea. We define

[TABLE]

where $\langle\dots\rangle_{0}$ is the average over the Fermi sea. Therefore $-T_{0}$ is the hopping energy per site, and $-2T_{0}/z$ the hopping energy per bond of the Fermi sea.

The saddle point equations for $\hat{\Phi}_{0}$ can be readily found

[TABLE]

The lowest energy eigenvalue is

[TABLE]

and is characterised by

[TABLE]

with $\tan\theta=8T_{0}\,R^{(0)}/U$. Since through Eq. (66) $R^{(0)}=\sin\theta$, the self-consistency condition implies

[TABLE]

namely

[TABLE]

$U_{c}$ is the well known value of the Brinkman-RiceBrinkman and Rice (1970) metal-insulator transition within the Gutzwiller approximation.

In conclusion, the lowest energy eigenstate is

[TABLE]

where $\cos\theta=\text{min}\left(1,U/U_{c}\right)$, and has eigenvalue

[TABLE]

We can now find all other eigenvalues and eigenvectors. The highest energy one is

[TABLE]

with eigenvalue

[TABLE]

This eigenstate actually corresponds to the high energy Hubbard bands.

The lowest excited eigenstate is threefold degenerate ($i=1,2,3$)

[TABLE]

with eigenvalue

[TABLE]

and describes spin fluctuations. We note that above the Brinkmann-Rice transition, $U>U_{c}$, this magnetic state becomes degenerate with the ground state. In what follows we shall anyway expand always around $\hat{\Phi}_{0}$, and, to avoid problems, we will mostly consider the metal phase at $U\leq U_{c}$.

Finally, the last eigenstate is

[TABLE]

with eigenvalue

[TABLE]

and describes instead charge fluctuations. This mode becomes degenerate with $\hat{\Phi}_{3}$ above the transition.

IV.2 Action of the fluctuations

Following section III we write

[TABLE]

with $\phi_{0}(i,t)$ fixed by normalisation. Through equations (62), (63) and (64) we find that

[TABLE]

where we have introduced the conjugate variables associated with $\phi_{n}$ and $\phi_{n}^{*}$. Eq. (44) reads explicitly

[TABLE]

where $\boldsymbol{\nabla}$ is the lattice divergence, $\mathbf{J}_{s}(i)$ and $\mathbf{J}_{c}(i)$ the spin and charge currents, respectively, defined through the continuity equations

[TABLE]

and finally $h_{*}(i)$ the Hamiltonian density

[TABLE]

Therefore $\delta_{1}E_{*}(t)$ defined in Eq. (43) becomes, due to particle-hole and spin $SU(2)$ symmetry

[TABLE]

where $\chi_{{}_{\boldsymbol{\nabla J\nabla J}}}$ is the linear response function of $\boldsymbol{\nabla J}$ with the Hamiltonian $\mathcal{H}^{(0)}_{*}$, which is actually the same for charge and spin currents, and $\chi_{h_{*}h_{*}}$ the response function of $h_{*}$. We observe that, because of charge and spin continuity equations, in Fourier space the following equivalence holds

[TABLE]

where $\chi(\mathbf{q},\omega)$ is the density-density response function, which is the same both in the charge and spin channels, and by definition

[TABLE]

Without going into further details, we find that the following expressions for the remaining contributions $\delta_{2}E_{*}(t)$ in Eq. (46), and $\delta_{4}E_{*}(t)$ in Eq. (49):

[TABLE]

We have now all ingredients required to evaluate linear response functions of local operators within the harmonic approximation for the fluctuations.

IV.3 Hubbard-band dispersion mode

As we mentioned, the Hubbard bands may be associated with the excited state $\hat{\Phi}_{3}$, hence with the operators $x_{3}$ and $p_{3}$. Their equations of motion in Fourier space are

[TABLE]

Within the metal phase, $U<U_{c}$, $\omega_{3}=E_{3}-E_{0}=4T_{0}$, so that, upon defining $\cos\theta=U/U_{c}\equiv u$, and noting that, for small $|\mathbf{q}|$, $\chi_{hh}(\mathbf{q},\omega)=O(q^{4})$, the eigenmode energy is solution of the equation

[TABLE]

thus describes an optical mode that softens at the metal insulator transition, $\omega_{3\mathbf{0}}=4T_{0}\,\sqrt{1-u^{2}\,}\,\to 0$ when $u\to 1$. We observe that the continuum of quasiparticle-quasihole excitations extends up to an energy of order T_{0}\big{(}1-u^{2}\big{)}, so that, upon approaching the transition, $\omega_{3\mathbf{q}}$ must detach from the continuum and become a genuine coherent excitation.

This coherent mode actually corresponds to the spin-wave excitations of the Ising field within the $Z_{2}$ slave-spin representation of the Hubbard modelHuber and Rüegg (2009); Rüegg et al. (2010); Schiró and Fabrizio (2011). This is not surprising since, as shown in Ref. Schiró and Fabrizio, 2011, the Gutzwiller wavefunction is just the mean-field variational state of the $Z_{2}$ slave-spin theory. At the mean-field level, the Mott transition in this representation translates into the order-disorder transition of a quantum Ising model. Therefore the mode $x_{3}$ seems to be the real fingerprint of the Mott transition.

IV.4 Dynamical charge susceptibility

We assume to perturb the system in the metal phase, $u\leq 1$, by an external potential that couples to the charge deviation from half-filling, namely

[TABLE]

Since \omega_{2}=E_{2}-E_{0}=2T_{0}\big{(}1+u\big{)} and by means of Eq. (86), we find in the presence of the field the following equations of motion for the conjugate variables $x_{2}$ and $p_{2}$

[TABLE]

from which it follows that the dynamical charge susceptibility is

[TABLE]

where it is evident the analogy with conventional RPA, though with a renormalised coupling constant

[TABLE]

We note that

[TABLE]

where

[TABLE]

is the quasiparticle density of states (DOS) at the chemical potential, as opposed to the bare DOS $\mathcal{N}_{0}$, and diverges approaching the Mott transition. Therefore, through Eq. (94), the charge compressibility is readily obtained

[TABLE]

and defines the Landau $F_{0}^{S}$ parameter

[TABLE]

Since approaching the transition, $u\to 1$, $F^{S}_{0}\sim(1-u)^{-2}$ diverges faster than $\mathcal{N}_{*}\sim(1-u)^{-1}$, we find that the charge compressibility correctly vanishes at the MIT. The expression of $F^{S}_{0}$ coincides with that originally obtained by VollhardtVollhardt (1984).

In the opposite limit of small $|\mathbf{q}|$ with respect to frequency,

[TABLE]

which, inserted into Eq. (94), allows calculating the poles of the dynamical charge susceptibility, which are

[TABLE]

This acoustic mode is above the quasiparticle-quasihole continuum and actually corresponds to the Landau’s zero sound. Once again this result is compatible with Vollhardt’s description of the correlated metal within the Gutzwiller approximation in the framework of Landau-Fermi liquid theoryVollhardt (1984). Indeed the zero sound velocity has the expected Landau’s expression, once one realises that in a lattice with infinite coordination $F_{1}^{S}=0$ and it is unrelated to the enhancement of the effective mass.

We conclude highlighting that the velocity of the zero sound stays constant approaching the Mott transition. In particular, for \omega^{2}\gg T_{0}\,\big{(}1-u^{2}\big{)}\,\big{(}\gamma_{\mathbf{0}}-\gamma_{\mathbf{q}}\big{)}, the dynamical charge susceptibility can be written as

[TABLE]

hence the pole at the zero sound has vanishing weight as the transition $u\to 1$ is approached, in agreement with the expectation that spectral weight is transferred at high energy.

We conclude by observing that the propagator $\Pi_{2}(\mathbf{q},\omega)$ of $p_{2}(\mathbf{q},\omega)$

[TABLE]

is singular at $\omega=0$, although this singularity does not appear in the physical response function, which is proportional to the propagator of the conjugate variable $x_{2}(\mathbf{q},\omega)$. Indeed, $p_{2}(\mathbf{q},\omega)$ is one of the would-be Goldstone modes that we mentioned in section III.1. The action of the single-band Hubbard model is $U(2)=U(1)\times SU(2)$ gauge invariant, and $p_{2}(\mathbf{q},\omega)$ is just the would-be Goldstone mode associated with the abelian $U(1)$, whereas we shall see that $\mathbf{p}_{1}(\mathbf{q},\omega)$ are instead those associated with $SU(2)$. In fact, the RPA form of the charge susceptibility could be very easily obtained by the gauge-fixing prescription of section III.1. If we drop all terms that contain $p_{2}(i,t)$ and replace

[TABLE]

we get an effective Hamiltonian of the quasiparticles, neglecting for convenience all other variables but $x_{2}(i,t)$,

[TABLE]

where

[TABLE]

which readily leads to Eq. (94).

IV.5 Dynamical spin susceptibility

In order to study the spin response, we imagine to add an external field that couples to the spin density, e.g. to its $z$ component, namely

[TABLE]

In the metal phase \omega_{1}=E_{1}-E_{0}=2T_{0}\big{(}1-u\big{)}, and repeating all calculations done for the charge susceptibility, we finally obtain the dynamical spin susceptibility

[TABLE]

where

[TABLE]

The above expression reproduces the small $u$ Stoner’s enhancement of the magnetic susceptibility. In addition it satisfies the relationship $\Gamma_{s}(U)=\Gamma_{c}(-U)$ valid at particle-hole symmetryVollhardt (1984). Since $\Gamma_{s}\sim(1-u)$ vanishes linearly approaching the transition, the Landau’s parameter

[TABLE]

is constant for $u\to 1$, which implies that the uniform static spin susceptibility diverges at the MIT. This result agrees with previous onesVollhardt (1984); Seibold et al. (2004) also obtained within the Gutzwiller approximation, but contrasts DMFT, which instead finds a finite uniform spin susceptibility at the transition.

Such negative outcome critically depends from the fact that the effective interaction $\Gamma_{s}$, Eq. (103), vanishes at the transition. We are going to show that beyond the harmonic approximation this cancellation does not occur anymore.

We note that $p_{1a}(i,t)$, $a=1,\dots,3$, are now the Goldstone modes associated with $SU(2)$ gauge invariance, and their propagators

[TABLE]

diverge at $\omega=0$. We can, as in section IV.4, drop $\mathbf{p}_{1}(i,t)$ from the action and replace

[TABLE]

whose effect could be absorbed into an effective magnetic field

[TABLE]

that straightforwardly leads to Eq. (102).

IV.6 Beyond RPA in the $x_{3}$ mode

We observe that all the above results in the metal phase correspond to expanding the action at second order in the fluctuations but treating the linear coupling between the latter and the fermions just within RPA, i.e. not accounting for exchange processes. While this procedure is somehow forced by gauge invariance for what it concerns charge and spin modes, see the ending parts of sections IV.4 and IV.5, it is not really compulsory for the $x_{3}(i,t)$ mode that describes the Hubbard bands. We can therefore take a first step forward when dealing with $x_{3}(i,t)$ in the direction of the so called RPA+Exchange. According to Eq. (81), promoting $x_{3}$ and $p_{3}$ to quantum conjugate variables, after defining $t_{*}=t\,\sin^{2}\theta$ and

[TABLE]

the Hamiltonian reads

[TABLE]

where the effective fields are those in Eqs. (100) and (105). The last term in Eq. (106), linear in $x_{3}$, derives from Eq. (III) and cancels the linear term of the hopping when the latter is averaged over the Fermi sea, which is just the saddle point condition for $x_{3}$.

Near the Mott transition from the metal side, $u\lesssim 1$, since $t_{*}$ is small with respect to $\omega_{3}$, we can integrate out $x_{3}$ and neglect the frequency dependence of its propagator $D_{3}(\mathbf{q},\omega)$, which, through Eqs. (90) and (91), implies that

[TABLE]

where we have furthermore neglected the momentum dependence.

In this approximation the mode $x_{3}$ simply induces a non-retarded electron-electron interaction, which, within RPA+Exchange, leads to a change of the charge and spin susceptibilities,

[TABLE]

where

[TABLE]

which also implies that the Landau parameters change into

[TABLE]

The charge $F_{0}^{S}>0$ keeps its singularity $(1-u)^{-2}$, so that the charge compressibility still vanishes. On the contrary,

[TABLE]

so that the uniform spin susceptibility

[TABLE]

is now finite. Remarkably, this expression agrees with that obtained by DMFTGeorges et al. (1996), although the numerical value of $U_{c}$ in DMFT is smaller than in the Gutzwiller approximation.

The quantum Hamiltonian (106) also allows calculating the optical conductivity. In the presence of a small transverse vector potential $A_{i\to j}(t)=-A_{j\to i}(t)$ the Hamiltonian acquires an additional term

[TABLE]

The calculation of the optical conductivity is straightforward, and follows exactly that obtained within slave-bosons in Ref. Raimondi and Castellani, 1993. Besides the Drude peak that is obtained taking $X(i)=1$, and vanishes like $\sin^{2}\theta=1-u^{2}$ at the transition, the optical conductivity gets high-frequency contributions from the absorption spectrum of the mode $x_{3}$Raimondi and Castellani (1993).

V Conclusions

In this paper we have presented a quite simple method to calculate linear response functions within the Gutzwiller approximation, including in a consistent way quantum fluctuations in the harmonic approximation. The calculation is straightforward and just requires a little more effort than the equilibrium one. In fact, besides the variational matrix $\hat{\Phi}_{0}$ that minimises the energy at equilibrium, and which can be regarded as the lowest energy eigenstate of a local HamiltonianLanatà et al. (2012, 2015), see Eq. (29), one also needs all excited eigenstates and eigenvalues. In a model that involves $M$ correlated orbitals in each unit cell, this local Hamiltonian is defined in a Hilbert space of dimension $\binom{4M}{2M}$, and can be conveniently recast into the problem of an impurity with $M$ orbitals hybridised to a single bath site with the same number of orbitals, the coupled system being at half-fillingLanatà et al. (2015).

As a check we have applied the method to the single-band Hubbard model at half-filling and recovered all known resultsVollhardt (1984); Jolicoeur and Le Guillou (1991); Raimondi and Castellani (1993); Seibold and Lorenzana (2001); Seibold et al. (2003, 2004); Schiró and Fabrizio (2011); Sandri et al. (2012). As a by-product, we also showed how to cure one flaw of the Gutzwiller approximation, i.e. the divergence of the uniform magnetic susceptibility approaching the Mott transition from the metal side.

Acknowledgments

This work has been supported by the European Union under H2020 Framework Programs, ERC Advanced Grant No. 692670 “FIRSTORM”.

Appendix A The wavefunction renormalisation matrix $\hat{R}(i)$

At equilibrium and in the natural basis, the constraint Eq. (11) reads

[TABLE]

where $\hat{P}_{0}(i)$ is the local probability distribution of the Slater determinant. Hereafter we shall drop for simplicity the site index $i$.

We can always write $\hat{P}^{(0)}_{0}$ as the Boltzmann distribution of a non-interacting Hamiltonian

[TABLE]

where f\big{(}\epsilon_{\alpha}\big{)}=n^{(0)}_{\alpha} is the Fermi distribution function. If $\hat{\Phi}$ is varied, also the probability distribution must vary in such a way as to preserve the constraint. This change will generally correspond to

[TABLE]

Since $H$ must still be a one body Hamiltonian it follows that

[TABLE]

where $\hat{H}$ is the matrix representation of $H$ in the single-particle basis, so that $d^{\phantom{\dagger}}_{\alpha}(\tau)$ remains a combination of creation operators. Since $\hat{U}(\tau_{1})\;\hat{U}(\tau_{2})=\hat{U}(\tau_{1}+\tau_{2})$, it trivially holds that $\hat{U}(\tau)\,\hat{U}(-\tau)=1$ and

[TABLE]

The local probability distribution

[TABLE]

so that

[TABLE]

namely

[TABLE]

which relates $\hat{U}(\beta)$ to $\hat{n}$. It also follows that

[TABLE]

The renormalisation coefficients $R$ is obtained by solving for any $\alpha$ and $\gamma$

[TABLE]

where

[TABLE]

Therefore, once we define

[TABLE]

then Eq. (115) is equivalent to

[TABLE]

or, in matrix form, and observing that $\hat{n}=\hat{U}(\beta)-\hat{n}\;\hat{U}(\beta)$,

[TABLE]

so that, multiplying both sides on the right by $\hat{U}(-\beta)$ we finally get

[TABLE]

We denote as

[TABLE]

since $\hat{S}=\hat{S}^{\dagger}$, so that

[TABLE]

namely the desired result

[TABLE]

One can rewrite

[TABLE]

where the matrix elements of $\hat{\Delta}$ are

[TABLE]

At equilibrium

[TABLE]

are diagonal, which allow an explicit evaluation of matrix derivatives. It follows that the equilibrium renormalisation matrix has elements

[TABLE]

A.1 Derivatives of $\hat{R}$

We write

[TABLE]

where $\hat{\Phi}_{n}$ is a basis set,

[TABLE]

with $\hat{\Phi}_{0}$ the equilibrium solution. By inspection we realise that

[TABLE]

where the tensor \hat{\Gamma}_{\alpha\beta}\Big{[}\hat{\Phi},\hat{\Phi}^{\dagger}\Big{]} is still functional of $\hat{\Phi}$ and $\hat{\Phi}^{\dagger}$. Therefore

[TABLE]

The equilibrium value is obtained by setting $\phi_{n}=\delta_{n0}$.

In particular, exploiting the fact that $\hat{S}$ is diagonal at equilibrium, the first order derivatives evaluated at equilibrium read explicitly

[TABLE]

while the second derivative, still calculated at equilibrium, is

[TABLE]

where

[TABLE]

The terms that appear in the above equations are

[TABLE]

and, lastly,

[TABLE]

In addition

[TABLE]

where the right hand sides are obtained straightforwardly through Eq. (125). The above derivatives calculated at the equilibrium solution allow calculating the Taylor expansion of $\hat{R}$. In particular, through equations (121) and (122), the first order expansion is

[TABLE]

while the second order expansion mentioned in Eq. (47), is

[TABLE]

where, explicitly,

[TABLE]

and

[TABLE]