Undamped Higgs Modes in Strongly Interacting Superconductors

Abstract

In superconductors, gauge $U\left(1\right)$ symmetry is spontaneously broken. According to Goldstone’s theorem, this breaking of a continuous symmetry establishes the existence of the Bogoliubov phase mode while the gauge-invariant response also includes the amplitude fluctuations of the order parameter. The latter, which are also termed ‘Higgs’ modes in analogy with the standard model, appear at the energy of the spectral gap $2\Delta$, when the superconducting ground state is evaluated within the weak-coupling BCS theory, and, therefore, are damped. Previously, we have shown that, within the time-dependent Gutzwiller approximation (TDGA), Higgs modes appear inside the gap with a finite binding energy relative to the quasiparticle continuum. Here, we show that the binding energy of the Higgs mode becomes exponentially small in the weak-coupling limit converging to the BCS solution. On the other hand, well-defined undamped amplitude modes exist in strongly coupled superconductors when the interaction energy becomes of the order of the bandwidth.

1. Introduction

There exists a close analogy between the mechanism at the heart of the mass generation for scalar and vector bosons in the standard model [1,2] and the mechanism responsible for superconductivity [3,4]. In both cases, the spontaneous breaking of a continuous symmetry comes along with the emergence of new elementary excitations. In the context of superconductors, they are identified with the collective fluctuations of the macroscopic (complex) order parameter. These comprise the massless phase fluctuations, which, in a charged superconductor, are pushed to the plasma frequency via the Anderson–Higgs mechanism [5]. In addition, amplitude fluctuations are the analog of the Higgs field in the ‘standard model’ and are, therefore, also referred to as “Higgs” modes in a superconducting system, cf. Figure 1a. In conventional (weakly-coupled or BCS) superconductors, the energy of the Higgs mode coincides with the spectral gap $2\Delta$ for single-particle excitations. The latter influence the dynamics of the Higgs excitation, leading to a non-relativistic strongly overdamped mode. In addition, the effects can hardly be probed experimentally, unless the system is strongly perturbed from its equilibrium [6,7,8,9,10,11,12].

It has been discussed in the literature (see, e.g., Refs. [13,14,15]) that amplitude fluctuations of the superconducting order parameter can give rise to in-gap excitations when the quantum critical point of the superconducting phase transition is approached. In this spirit, the subgap absorption, observed in strongly disordered NbN and InO films [16], has been taken as evidence for Higgs modes, which appear below the spectral gap $2\Delta$. However, this interpretation, which relies on the validity of Lorentz-invariant $\mathcal{O}\left(2\right)$ models, has been questioned in Ref. [17] and, instead, it has been shown that the subgap absorption, observed in Ref. [16] and also Ref. [18], is most likely due to phase modes, which, due to the low superfluid density (and, thus, low energy plasma excitation), can be active at energies below the spectral gap.

In a recent paper [19], we have shown that, in the limit of strong coupling, the Higgs mode is shifted inside the spectral gap, even in the clean system, which leads to undamped Higgs oscillations, cf. Figure 1b. This theory is relevant for systems where the interaction becomes of the same order than that of the bandwidth. Such a situation could be realized in cold atomic Fermion condensates and, in fact, a well-defined collective mode throughout the entire crossover from BCS to Bose–Einstein condensation has been observed recently [20] in a system with weak modulation of the confinement. Our previous work was based on the time-dependent Gutzwiller approximation (TDGA) [21,22,23,24,25,26,27,28,29], which allows for the investigation of dynamical properties of systems with a strong local interaction U (i.e., Hubbard-type models). By evaluating the spectral gap $2\Delta$ and the energy of the Higgs mode, it could be shown that both quantities begin to separate when the magnitude of the local attraction $\left|U\right|$ becomes larger than the bandwidth. However, within this numerical approach, it is difficult to analyze the situation at weak coupling. In particular, it is difficult to distinguish between the situation where a critical coupling must be exceeded to push the Higgs mode into the spectral gap, and a scenario where, at weak coupling, the difference between the energy of the Higgs mode and $2\Delta$ becomes exponentially small.

Here, we analyze the situation in the weak-coupling limit of the TDGA and show that the latter scenario is realized in this case. Section 2 introduces the model and first discusses the appearance of amplitude modes within the standard BCS + RPA approach. The analysis of the Higgs excitation within the TDGA is then performed in Section 4 in the weak-coupling limit and we finally conclude our discussion in Section 5.

2. Model and BCS Approximation

We exemplify the evaluation of amplitude modes within the single-band attractive Hubbard Hamiltonian

(1)$\begin{array}{ccc}\hfill H& =& \sum _{ij}{t}_{ij}{c}_{i,\sigma }^{†}{c}_{j,\sigma }-\mu \sum _i{n}_i-\left|U\right|\sum _i{n}_{i,\uparrow }{n}_{i,\downarrow }\hfill \end{array}$

where, in the following, we denote the Fourier transform of $t_{ij}$ as the dispersion ${\epsilon }_k$.

Decoupling in the pair channel yields the usual BCS Hamiltonian

(2)$H=\sum _{k,\sigma }{\xi }_k{c}_{k,\sigma }^{†}{c}_{k,\sigma }+\sum _k\left({\Delta }^{*}{c}_{-k,\downarrow }{c}_{k,\uparrow }+\Delta c_{k,\uparrow }^{†}{c}_{-k,\downarrow }^{†}\right)+N{\displaystyle \frac{{|\Delta |}^2}{\left|U\right|}}$

with $\Delta =-\left|U\right|⟨c_{-k,\downarrow }{c}_{k,\uparrow }⟩$, ${\xi }_k={\epsilon }_k-\mu$, and N denotes the number of lattice sites.

Equation (2) can be diagonalized with the Bogoljubov transformation,

(3)$\begin{array}{ccc}\hfill c_{k,\uparrow }& =& {\alpha }_k{\gamma }_{k,0}-{\beta }_k^{*}{\gamma }_{k,1}^{†}\hfill \end{array}$

(4)$\begin{array}{ccc}\hfill c_{-k,\downarrow }& =& {\beta }_k{\gamma }_{k,0}^{†}+{\alpha }_k{\gamma }_{k,1}\hfill \end{array}$

and one obtains

(5)$H=\sum _k{E}_K\left({\gamma }_{k,0}^{†}{\gamma }_{k,0}+{\gamma }_{k,1}^{†}{\gamma }_{k,1}\right)+\sum _k({\xi }_k-E_k)+N{\displaystyle \frac{|{\Delta }^2|}{\left|U\right|}}$

where $E_k=\sqrt{{\xi }_k^2+{\Delta }^2}$, ${\alpha }_k^2={\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{{\xi }_k}{E_k}}\right)$, and ${\beta }_k^2={\displaystyle \frac{1}{2}}\left(1-{\displaystyle \frac{{\xi }_k}{E_k}}\right)$. The self-consistency condition (zero temperature) reads

(6)${\displaystyle \frac{1}{\left|U\right|}}={\displaystyle \frac{1}{N}}\sum _k{\displaystyle \frac{1}{2E_k}}.$

3. Collective Modes beyond Weak-Coupling BCS Theory

The mean-field decoupling neglects the following fluctuation contributions in the pairing channel:

(7)$V_{fl}^{pair}=-\left|U\right|{\displaystyle \frac{1}{N}}\sum _q\delta {\Delta }_q^{†}\delta {\Delta }_q$

with

(8)$\begin{array}{ccc}\hfill \delta {\Delta }_q& =& \sum _k\left[c_{-k+q/2,\downarrow }{c}_{k+q/2,\uparrow }-⟨c_{-k+q/2,\downarrow }{c}_{k+q/2,\uparrow }⟩\right]\hfill \\ & =& \sum _k\left[{\alpha }_k^{+}{\alpha }_k^{-}{\gamma }_{k-q/2,1}{\gamma }_{k+q/2,0}-{\beta }_k^{+}{\beta }_k^{-}{\gamma }_{k-q/2,0}^{†}{\gamma }_{k+q/2,1}^{†}\right.\hfill \\ & +& \left.{\alpha }_k^{+}{\beta }_k^{-}{\gamma }_{k-q/2,0}^{†}{\gamma }_{k+q/2,0}+{\beta }_k^{+}{\alpha }_k^{-}{\gamma }_{k+q/2,1}^{†}{\gamma }_{k-q/2,1}\right]\hfill \end{array}$

(9)$\begin{array}{ccc}\hfill \delta {\Delta }_q^{†}& =& \left[c_{k+q/2,\uparrow }^{†}{c}_{-k+q/2,\downarrow }^{†}-⟨c_{k+q/2,\uparrow }^{†}{c}_{-k+q/2,\downarrow }^{†}⟩\right]\hfill \\ & =& \sum _k\left[{\alpha }_k^{+}{\alpha }_k^{-}{\gamma }_{k+q/2,0}^{†}{\gamma }_{k-q/2,1}^{†}-{\beta }_k^{+}{\beta }_k^{-}{\gamma }_{k+q/2,1}{\gamma }_{k-q/2,0}\right.\hfill \\ & +& \left.{\alpha }_k^{+}{\beta }_k^{-}{\gamma }_{k+q/2,0}^{†}{\gamma }_{k-q/2,0}+{\beta }_k^{+}{\alpha }_k^{-}{\gamma }_{k-q/2,1}^{†}{\gamma }_{k+q/2,1}\right],\hfill \end{array}$

and we have defined ${\alpha }_k^{\pm }={\alpha }_{k\pm q/2}$ and ${\beta }_k^{\pm }={\beta }_{k\pm q/2}$.

In the same way, the fluctuation contribution in the charge channel is obtained as

(10)$V_{fl}^{charge}=-{\displaystyle \frac{\left|U\right|}{2}}{\displaystyle \frac{1}{2N}}\sum _q\delta {\rho }_q\delta {\rho }_{-q}$

with

(11)$\begin{array}{ccc}\hfill \delta {\rho }_q& =& \sum _k\left[({\alpha }_k^{+}{\beta }_k^{-}+{\alpha }_k^{-}{\beta }_k^{+})({\gamma }_{k-q/2,1}^{†}{\gamma }_{k+q/2,0}^{†}+{\gamma }_{k-q/2,0}{\gamma }_{k+q/2,1})\right.\hfill \\ & +& \left.({\alpha }_k^{+}{\alpha }_k^{-}-{\beta }_k^{+}{\beta }_k^{-})({\gamma }_{k-q/2,1}^{†}{\gamma }_{k+q/2,1}+{\gamma }_{k+q/2,0}^{†}{\gamma }_{k-q/2,0})\right].\hfill \end{array}$

3.1. Correlation Functions and RPA Resummation

We denote correlation functions by

(12)${\chi }_{nm}(\widehat{A},\widehat{B})=-i\int dt{e}^{i\omega t}⟨\mathcal{T}{\widehat{A}}_n\left(t\right){\widehat{B}}_m\left(0\right)⟩$

where, in the following, $\widehat{A}$ and $\widehat{B}$ correspond to either pair or charge fluctuations in Equations (8), (9), and (11).

It is convenient to define the $3\times 3$ matrices

$\underline{{\chi }_q^0}\left(\omega \right)=\left(\begin{array}{ccc}{\chi }_q^0(\delta {\Delta }_{-q}^{†},\delta {\Delta }_q^{†})& {\chi }_q^0(\delta {\Delta }_{-q}^{†},\delta {\Delta }_{-q})& {\chi }_q^0(\delta {\Delta }_{-q}^{†},\delta {\rho }_{-q})\\ {\chi }_q^0(\delta {\Delta }_q,\delta {\Delta }_q^{†})& {\chi }_q^0(\delta {\Delta }_q,\delta {\Delta }_{-q})& {\chi }_q^0(\delta {\Delta }_q,\delta {\rho }_{-q})\\ {\chi }_q^0(\delta {\rho }_q,\delta {\Delta }_q^{†})& {\chi }_q^0(\delta {\rho }_q,\delta {\Delta }_{-q})& {\chi }_q^0(\delta {\rho }_q,\delta {\rho }_{-q})\end{array}\right)$

where the superscript ‘0’ indicates that the elements are computed with the bare BCS wave function and are given in Appendix A.

The correlations of the interacting system are then defined by an analogous matrix $\underline{{\chi }_q}\left(\omega \right)$ but without superscript. From Equations (7) and (10), the matrix for the local interaction is given by

(13)$\underline{V}=\left(\begin{array}{ccc}0& -\left|U\right|& 0\\ -\left|U\right|& 0& 0\\ 0& 0& -\left|U\right|/2\end{array}\right).$

The RPA resummation can then be written as

$\underline{\chi }=\underline{{\chi }^0}+\underline{{\chi }^0}\underline{V}\underline{\chi }$

which is solved by

(14)$\underline{\chi }={\left[\underline{1}-\underline{{\chi }^0}\underline{V}\right]}^{-1}\underline{{\chi }^0}.$

3.2. Amplitude and Phase Correlations

Introducing amplitude $A_q\equiv (\delta {\Delta }_q+\delta {\Delta }_{-q}^{†})/\sqrt{2}$ and phase operators ${\Phi }_q\equiv (\delta {\Delta }_q-\delta {\Delta }_{-q}^{†})/\sqrt{2}$, the amplitude and phase correlation functions are obtained from

(15)$\begin{array}{ccc}\hfill {\chi }_q^{AA}& \equiv & {\chi }_q(A_q,A_{-q})={\displaystyle \frac{1}{2}}\left[{\chi }_q(\delta {\Delta }_{-q}^{†},\delta {\Delta }_q^{†})+{\chi }_q(\delta {\Delta }_{-q}^{†},\delta {\Delta }_{-q})\right.\hfill \\ & +& \left.{\chi }_q(\delta {\Delta }_q,\delta {\Delta }_q^{†})+{\chi }_q(\delta {\Delta }_q,\delta {\Delta }_{-q})\right]\hfill \end{array}$

(16)$\begin{array}{ccc}\hfill {\chi }_q^{\Phi \Phi }& \equiv & {\chi }_q({\Phi }_q,{\Phi }_{-q})={\displaystyle \frac{1}{2}}\left[{\chi }_q(\delta {\Delta }_{-q}^{†},\delta {\Delta }_q^{†})-{\chi }_q(\delta {\Delta }_{-q}^{†},\delta {\Delta }_{-q})\right.\hfill \\ & -& \left.{\chi }_q(\delta {\Delta }_q,\delta {\Delta }_q^{†})+{\chi }_q(\delta {\Delta }_q,\delta {\Delta }_{-q})\right]\hfill \end{array}$

and analogously for the mixed correlations between amplitude, phase, and charge.

The susceptibility and interaction matrices can now be transformed into the amplitude-phase representation by introducing the matrix

(17)$\underline{\gamma }=\left(\begin{array}{ccc}{\displaystyle \frac{1}{\sqrt{2}}}& {\displaystyle \frac{1}{\sqrt{2}}}& 0\\ {\displaystyle \frac{1}{\sqrt{2}}}& -{\displaystyle \frac{1}{\sqrt{2}}}& 0\\ 0& 0& 1\end{array}\right)$

so that

(18)$\begin{array}{ccc}\hfill \underline{\gamma }{\underline{\chi }}_q\underline{\gamma }& =& \left(\begin{array}{ccc}{\chi }_q^{AA}& {\chi }_q^{A\Phi }& {\chi }_q^{A\rho }\\ {\chi }_q^{\Phi A}& {\chi }_q^{\Phi \Phi }& {\chi }_q^{\Phi \rho }\\ {\chi }_q^{\rho A}& {\chi }_q^{\Phi \rho }& {\chi }_q^{\rho \rho }\end{array}\right)\hfill \end{array}$

(19)$\begin{array}{ccc}\hfill \underline{\gamma }\underline{V}\underline{\gamma }& =& \left(\begin{array}{ccc}-\left|U\right|& 0& 0\\ 0& \left|U\right|& 0\\ 0& 0& -\left|U\right|/2\end{array}\right)\hfill \end{array}$

and the RPA equation, Equation (14), also holds in the transformed representation. Therefore, the interaction in the amplitude channel is attractive, whereas it is repulsive in the phase channel.

3.3. The ph-Symmetric Case in the Limit q = 0

In the case of particle-hole symmetry and for $q=0$ (where ${\alpha }_k^{\pm }={\alpha }_k$, and ${\beta }_k^{\pm }={\beta }_k$), one finds a decoupling of amplitude from phase and charge fluctuations, i.e.,

(20)$\begin{array}{ccc}\hfill {\chi }_q^{\left(0\right),A\Phi }\left(\omega \right)& =& -{\displaystyle \frac{1}{2N}}\sum _k{\displaystyle \frac{{\xi }_k}{E_k}}\left[{\displaystyle \frac{1}{\omega -2E_k}}+{\displaystyle \frac{1}{\omega +2E_k}}\right]=0\hfill \end{array}$

(21)$\begin{array}{ccc}\hfill {\chi }_q^{\left(0\right),A\rho }\left(\omega \right)& =& -{\displaystyle \frac{\Delta }{\sqrt{2}N}}\sum _k{\displaystyle \frac{{\xi }_k}{E_k^2}}\left[{\displaystyle \frac{1}{\omega -2E_k}}-{\displaystyle \frac{1}{\omega +2E_k}}\right]=0.\hfill \end{array}$

As a consequence, the amplitude correlations decouple from the phase-charge sector in the RPA equation, Equation (14), and one finds

(22)${\chi }_{q=0}^{AA}\left(\omega \right)={\displaystyle \frac{{\chi }_{q=0}^{\left(0\right)AA}\left(\omega \right)}{1+\left|U\right|{\chi }_{q=0}^{\left(0\right)AA}\left(\omega \right)}}$

and, from Equation (15) together with the bare correlation functions listed in Appendix A, one obtains

(23)${\chi }_{q=0}^{\left(0\right)AA}\left(\omega \right)={\displaystyle \frac{2}{N}}\sum _k{\displaystyle \frac{{\xi }_k^2}{E_k}}{\displaystyle \frac{1}{{\omega }^2-4E_k^2}}.$

Noting that

(24)${\chi }_{q=0}^{\left(0\right)AA}(\omega =2\Delta )={\displaystyle \frac{2}{N}}\sum _k{\displaystyle \frac{{\xi }_k^2}{E_k}}{\displaystyle \frac{1}{4{\Delta }^2-4E_k^2}}=-{\displaystyle \frac{1}{N}}\sum _k{\displaystyle \frac{1}{2E_k}}$

it becomes apparent that the ‘pole’ of Equation (22) is given by

(25)$1+\left|U\right|{\chi }_{q=0}^{\left(0\right)AA}(\omega =2\Delta )=1-{\displaystyle \frac{\left|U\right|}{N}}\sum _k{\displaystyle \frac{1}{2E_k}}=0$

because this is identical to the self-consistency equation, Equation (6).

Therefore, the energy of the amplitude mode at $q=0$ is given by $\omega =2\Delta$ and coincides with the onset of quasiparticle excitations. Thus, it is not a pole but rather a branch cut in the amplitude correlation function.

It should also be noted that a similar result appears in the case of the repulsive Hubbard model where the amplitude excitations of the SDW order parameter appear at the energy of the SDW gap, as has been shown in Ref. [30].

Figure 2 visualizes the condition to find the Higgs pole in Equation (22) for a two-dimensional system. As a function of $\omega$, the (negative) bare amplitude correlation function Equation (24) is a continuously decreasing function, which, at $\omega =2\Delta$, reaches the value ${\chi }_{q=0}^{\left(0\right)AA}(\omega =2\Delta )=-1/\left|U\right|$, which leads to the Higgs pole at exactly $\omega =2\Delta$. This can change if we hypothesize an interaction between quasiparticles V different from the static interaction U that appears in the BCS equation. In case V would be more negative than U (blue dashed line in Figure 2), the pole of Equation (22) would occur at $\omega <2\Delta$ and, therefore, would correspond to a bound state inside the spectral gap. We will show, in the following section, that this is exactly the situation that is realized within the TDGA because, in this variational scheme, the denominator of Equation (22) is not related to the self-consistency condition for the spectral gap when $\omega =2\Delta$. This can be understood as vertex corrections in the effective interaction between quasiparticles in the same spirit as Ref. [31].

4. TDGA

The TDGA ground state is obtained by optimizing the number of doubly occupied sites in the underlying BCS state $|BCS⟩$ by applying the ‘Gutzwiller projector’ ${\widehat{P}}_G$, i.e.,

$|{\Psi }_G⟩={\widehat{P}}_G|BCS⟩.$

The variational energy $E^{GA}(D,J^{-})=⟨{\Psi }_G|\widehat{H}|{\Psi }_G⟩/⟨{\Psi }_G|{\Psi }_G⟩$ depends on the double occupancy D and the anomalous correlations $J^{-}\equiv ⟨c_{i\downarrow }{c}_{i\uparrow }⟩$, cf., e.g., [32]. A first step in the application of the TDGA is the evaluation of the so-called Gutzwiller Hamiltonian, defined as the derivative of $E^{GA}(D,J^{-})$ with respect to the density matrix [21,22]. One obtains (for the homogeneous system, N lattice sites) [32]

$\begin{array}{ccc}\hfill H^{GA}& =& \sum _k{E}_k\left[{\gamma }_{k,0}^{†}{\gamma }_{k,0}+{\gamma }_{k,1}^{†}{\gamma }_{k,1}\right]+\sum _k({\xi }_k-E_k)-2N\Delta J^{-}-N\left|U\right|D.\hfill \end{array}$

where, now, the quasiparticle energy in ${\xi }_k=z^2{\epsilon }_k-\mu$ is renormalized by

(26)$z={\displaystyle \frac{\sqrt{\frac{1}{2}-D+J^z}\left(\sqrt{D-J^z-J^{-}}+\sqrt{D-J^z+J^{-}}\right)}{\sqrt{\frac{1}{4}-{\left(J^{-}\right)}^2}}}.$

and $J^z$ is related to the density n via $J^z=1/2(n-1)$. In contrast to weak-coupling BCS theory, where the spectral gap parameter is given by $\Delta =-\left|U\right|J^{-}$, in the TDGA, this quantity is obtained from the variational principle as

(27)$\Delta ={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathsf{\partial }{z}^2}{\mathsf{\partial }⟨J^{-}⟩}}{\displaystyle \frac{1}{N}}\sum _k{\epsilon }_k\left[1-{\displaystyle \frac{{\xi }_k}{E_k}}\right].$

On the other hand, the equation for the anomalous correlations still resembles the corresponding BCS result and reads

(28)$⟨J^{-}⟩=-{\displaystyle \frac{1}{N}}\sum _k{\displaystyle \frac{|\Delta |}{2E_k}}.$

4.1. Ground State for the Half-Filled Case

Similar to the BCS case, we focus in the following on the half-filled system where analytical results can be obtained in the weak-coupling limit. The renormalization factor simplifies to

(29)$z^2={\displaystyle \frac{1-2D}{\frac{1}{4}-{\left(J^{-}\right)}^2}}\left[D+\sqrt{D^2-{\left(J^{-}\right)}^2}\right]$

and expanding Equation (29) for small $J^{-}$ (weak coupling) yields

(30)$z^2=8D(1-2D)+32D(1-2D)(1-{\displaystyle \frac{1}{16D^2}}){\left(J^{-}\right)}^2.$

Furthermore, for the half-filled system, we can write $D=1/4+d$ and obtain in the lowest order in $J^{-}$,

(31)$z^2=1-16d^2+32d{\left(J^{-}\right)}^2.$

We proceed by evaluating Equations (27) and (28) for a constant DOS, $\rho \left(\omega \right)=1/\left(2B\right)$ for $-B\le \omega \le B$. For Equation (28), we obtain

(32)$\begin{array}{ccc}\hfill \Delta & =& -32d{J}^{-}{\displaystyle \frac{1}{2B}}{\int }_{-B}^{B}d\omega {\displaystyle \frac{{\omega }^2}{\sqrt{{\omega }^2+{\Delta }^2}}}\hfill \\ & =& -32d{J}^{-}{\displaystyle \frac{1}{2B}}\left[B\sqrt{B^2+{\Delta }^2}-{\displaystyle \frac{{\Delta }^2}{2}}ln{\displaystyle \frac{\sqrt{1+{\Delta }^2/B^2}+1}{\sqrt{1+{\Delta }^2/B^2}-1}}\right]\hfill \\ & \approx & -16d{J}^{-}\left[B-{\displaystyle \frac{{\Delta }^2}{B}}ln{\displaystyle \frac{2B}{\Delta }}\right]\approx -16Bd{J}^{-},\hfill \end{array}$

and, in the last line, we used the limit of weak coupling.

In the same limit, the anomalous correlations Equation (27) are obtained as

(33)$\begin{array}{ccc}\hfill ⟨J^{-}⟩& =& -{\displaystyle \frac{\Delta }{4B}}{\int }_B^B{\displaystyle \frac{d\omega }{\sqrt{{\omega }^2+{\Delta }^2}}}\hfill \\ & =& -{\displaystyle \frac{\Delta }{4B}}ln{\displaystyle \frac{\sqrt{1+{\Delta }^2/B^2}+1}{\sqrt{1+{\Delta }^2/B^2}-1}}\hfill \\ & \approx & -{\displaystyle \frac{\Delta }{2B}}ln{\displaystyle \frac{2B}{\Delta }}.\hfill \end{array}$

Inserting Equation (33) into Equation (32) finally yields

(34)$\begin{array}{ccc}\hfill \Delta & =& 2B{e}^{-{\displaystyle \frac{1}{8d}}}\hfill \end{array}$

(35)$\begin{array}{ccc}\hfill J^{-}& =& -{\displaystyle \frac{1}{8d}}{e}^{-{\displaystyle \frac{1}{8d}}},\hfill \end{array}$

where it should be noted that $B=B_0{z}^2$, with $B_0$ being the bare bandwidth. The dependence of $\Delta$ and $J^{-}$ on U is encoded in the dependence on d with $\Delta \to 0$, $J^{-}\to 0$ when $d\to 0$.

We are now in the position to evaluate the total Gutzwiller approximated energy

$\begin{array}{ccc}\hfill {\displaystyle \frac{E^{GA}}{N}}& =& -{\displaystyle \frac{1}{2B}}{\int }_{-B}^{B}d\omega \sqrt{{\omega }^2+{\Delta }^2}-2\Delta J^{-}+UD\hfill \\ & =& -{\displaystyle \frac{1}{4B}}\left[2B\sqrt{B^2+{\Delta }^2}+{\Delta }^{2}ln{\displaystyle \frac{\sqrt{B^2+{\Delta }^2}+B}{\sqrt{B^2+{\Delta }^2}-B}}\right]-2\Delta J^{-}+UD\hfill \\ & \approx & -B/2-2\Delta J^{-}+UD\hfill \\ & \approx & -{\displaystyle \frac{B_0}{2}}\left(1-16d^2+32d{\left(J^{-}\right)}^2\right)+{\displaystyle \frac{B_0}{2d}}\left(1-16d^2+32d{\left(J^{-}\right)}^2\right)e^{-{\displaystyle \frac{1}{4d}}}-\left|U\right|\left({\displaystyle \frac{1}{4}}+d\right)\hfill \\ & =& -{\displaystyle \frac{B_0}{2}}\left(1-16d^2+{\displaystyle \frac{1}{2d}}{e}^{-{\displaystyle \frac{1}{4d}}}\right)+{\displaystyle \frac{B_0}{2d}}\left(1-16d^2+{\displaystyle \frac{1}{2d}}{e}^{-{\displaystyle \frac{1}{4d}}}\right)e^{-{\displaystyle \frac{1}{4d}}}-\left|U\right|\left({\displaystyle \frac{1}{4}}+d\right)\hfill \end{array}$

and the minimization $\mathsf{\partial }{E}^{GA}/\mathsf{\partial }d=0$ leads to the equation

(36)$\begin{array}{ccc}& & \underset{1}{\underbrace{16d}}\underset{2}{\underbrace{-16e^{-{\displaystyle \frac{1}{4d}}}}}+\underset{3}{\underbrace{{\displaystyle \frac{1}{4d^2}}\left[{\displaystyle \frac{1}{4d}}-1\right]\left[{\displaystyle \frac{1}{d}}{e}^{-{\displaystyle \frac{1}{4d}}}-1\right]e^{-{\displaystyle \frac{1}{4d}}}}}\hfill \\ & +& \underset{4}{\underbrace{{\displaystyle \frac{1}{2d^2}}\left[{\displaystyle \frac{1}{4d}}-1\right]\left[1-16d^2+{\displaystyle \frac{1}{2d}}{e}^{-{\displaystyle \frac{1}{4d}}}\right]e^{-{\displaystyle \frac{1}{4d}}}}}={\displaystyle \frac{\left|U\right|}{B_0}}.\hfill \end{array}$

Contribution ‘1’ determines the double occupancy in the absence of SC, i.e., $d_0={\displaystyle \frac{\left|U\right|}{16B_0}}$.

On the other hand, it can be seen from Figure 3 that, for small d, the contributions ‘2–4’ in Equation (36) are exponentially small; so, by setting $d=d_0+ϵ$, the leading correction for small $U/B_0$ is given by

(37)$ϵ\approx -{\displaystyle \frac{1}{128}}{\displaystyle \frac{e^{-\frac{1}{2d_0}}}{d_0^4}}=-512{\left({\displaystyle \frac{B_0}{\left|U\right|}}\right)}^4{e}^{-{\displaystyle \frac{8B_0}{\left|U\right|}}},$

i.e., the double occupancy decreases with the onset of the SC order.

4.2. TDGA for the SC Half-Filled Case

Within the TDGA and employing the weak-coupling limit, the effective interaction in the amplitude channel can be evaluated by expanding the renormalized kinetic energy (per site)

(38)$e_k=-{\displaystyle \frac{B}{2}}=-{\displaystyle \frac{1}{2}}{B}_0{z}^2=-{\displaystyle \frac{1}{2}}{B}_0\left[1-16d^2+32d{\left(J^{-}\right)}^2\right]$

in terms of the fluctuations $\delta d$ and $\delta J^{-}$.

One obtains

(39)$H^{int}={\displaystyle \frac{1}{2}}\left(\begin{array}{c}\delta d\\ \delta J^{-}\end{array}\right)\left(\begin{array}{cc}16B_0& -32B_0{J}^{-}\\ -32B_0{J}^{-}& -32d{B}_0\end{array}\right)\left(\begin{array}{c}\delta d\\ \delta J^{-}\end{array}\right).$

In the spirit of the anti-adiabaticity principle [21], we eliminate the high-energy double occupancy fluctuations from Equation (39) using the condition $\mathsf{\partial }{H}^{int}/\mathsf{\partial }\delta d=0$, which yields $\delta d=2J^{-}\delta J^{-}$. The resulting interaction in the amplitude channel alone (it should be noted that $\delta J^{-}=\sqrt{2}\delta A_{q=0}$, cf. Section 3.2) is given by $H_{int}^{AA}=1/2V_{eff}\delta A_{q=0}\delta A_{q=0}$ with

(40)$V^{eff}=-16B_0\left[d+2{\left(J^{-}\right)}^2\right]$

or, when we use $d=d_0+\epsilon$,

(41)$V^{eff}=-\left|U\right|-16B_0\left[\epsilon +2{\left(J^{-}\right)}^2\right].$

From Equations (35) and (37), it turns out that, in the considered limit (half-filling, weak coupling), the TDGA provides an exponentially small correction to the bare interaction $-\left|U\right|$.

We proceed by investigating the consequences for the resulting Higgs mode within the TDGA. Similar to the BCS + RPA case, the amplitude mode occurs as a pole in Equation (22), $1+|V_0^{eff}|{\chi }_{q=0}^{\left(0\right)AA}\left(\omega \right)=0$ with

(42)${\chi }_{q=0}^{\left(0\right)AA}\left(\omega \right)={\displaystyle \frac{2}{N}}\sum _k{\displaystyle \frac{{\xi }_k^2}{E_k}}{\displaystyle \frac{1}{{\omega }^2-4E_k^2}}.$

where, now, ${\xi }_k$ contains the renormalization of ${\epsilon }_k$ by the GA factors $z^2$.

The comparison of Equation (25) and Equation (27) reveals that the Higgs mode would occur at $\omega =2\Delta$ when the effective interaction would be given by ${\tilde{V}}_{eff}=\Delta /⟨J^{-}⟩$. From Equation (32), one finds

(43)${\displaystyle \frac{\Delta }{⟨J^{-}⟩}}=-16B_{0}d{z}^2=-16B_{0}d\left[1-16d^2+32d{\left(J^{-}\right)}^2\right]$

or, when we make use of $d=d_0+\epsilon$,

(44)${\displaystyle \frac{\Delta }{⟨J^{-}⟩}}=-\left|U\right|+16d\left|U\right|\left[d-2{\left(J^{-}\right)}^2\right]-16B_0\epsilon \left[1-16d^2+32d{\left(J^{-}\right)}^2\right].$

Thus,

(45)$V^{eff}-{\displaystyle \frac{\Delta }{⟨J^{-}⟩}}=-32B_0{\left(J^{-}\right)}^2\left[1-16d^2\right]-256B_0{d}^3\approx -256B_0{d}_0^3=-{\displaystyle \frac{1}{16}}{B}_0{\left({\displaystyle \frac{\left|U\right|}{B_0}}\right)}^3<0.$

i.e., $V^{eff}$ is more negative than $\Delta /J^{-}$ in the small coupling limit so that, according to the analysis related to Figure 2, one obtains the pole for the Higgs mode below the spectral gap $2\Delta$.

Expanding the amplitude correlation function Equation (42) around the spectral gap

(46)${\chi }_{q=0}^{\left(0\right)AA}(\omega =2\Delta -\nu )\approx {\chi }_{q=0}^{\left(0\right)AA}(\omega =2\Delta )+{\displaystyle \frac{\pi }{4B_0}}\sqrt{{\displaystyle \frac{\nu }{\Delta }}}$

yields the pole from the resonance condition $1+|V^{eff}|{\chi }_{q=0}^{\left(0\right)AA}(\omega =2\Delta -\nu )=0$ at

(47)$\nu =2\Delta -{\omega }_H\sim \Delta {\displaystyle \frac{\left|U\right|}{B_0}}\sim \left|U\right|exp\left(-{\displaystyle \frac{2B_0}{\left|U\right|}}\right)$

i.e., the shift of the Higgs mode inside the spectral gap is exponentially small in weak coupling.

5. Conclusions

We have shown that the TDGA applied to the weak-coupling limit of the attractive Hubbard model leads to amplitude (‘Higgs’) modes which are split off from the quasiparticle continuum at $2\Delta$ but the corresponding binding energy is exponentially small in the attractive interaction. This is different from the BCS + RPA approach, where the ‘Higgs’ excitation always appears at $2\Delta$ and, therefore, is damped due to the interference with quasiparticle excitations. The analysis in the present paper is based on a constant DOS, which allowed for an analytical treatment of the problem in the weak-coupling limit. Using a numerical evaluation, we have already shown [19] that the DOS has some influence on the position of the Higgs mode and it would definitely be interesting to examine this dependency in more detail, and, eventually, also analytically, in future work. When the interaction becomes comparable with the bandwidth, the Higgs binding energy becomes sizeable [19] and, thus, long-lived amplitude modes could be observed in the crossover regime from BCS to BEC. In fact, recent investigations on cold atomic fermion condensates [20] confirm this picture, although, in this experiment, the modulation of the confinement also seems to play a role. Our theory also goes beyond the applicability to superconductors and should also have relevance in magnetic systems where the collective modes comprise the amplitude excitation of the magnetic order. In fact, in Ref. [19], we have proposed an experimental setup where the predicted bound state of the antiferromagnetic amplitude mode in undoped cuprates may be observed via magneto-optical methods. Of course, the main difference between superconducting and magnetic systems is the absence of low-lying excitations in the former due to the Anderson–Higgs mechanism, which pushes the Bogoljubov–Goldstone mode to the plasma energy. On the other hand, in a magnetic system, the decay channel of the amplitude mode into low-lying magnon excitations persists. In this regard, it is interesting that strong correlations can also stabilize quantum matter [33] so that the proposed split-off magnetic amplitude modes may have a lifetime long enough to be detected via a frequency-dependent Faraday rotated optical signal [19].

Author Contributions

Conceptualization, J.L. and G.S.; methodology, G.S.; validation J.L. and G.S.; writing, J.L. and G.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

We thank Lara Benfatto, Claudio Castellani, Mattia Udina, Paolo Barone, O.P. Sushkov, and Dirk Manske for useful discussions.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. (a) Ginzburg–Landau (GL) potential landscape of a superconductor. Bogoljubov phase excitations ([Forumla omitted. See PDF.]) and amplitude ‘Higgs’ (H) modes are indicated by yellow and blue arrows, respectively. (b) Evolution of the spectral gap and the energy of the Higgs mode as a function of the interaction parameter [Forumla omitted. See PDF.] (in terms of the bandwidth) for a two-dimensional system. In the regime [Forumla omitted. See PDF.], the Higgs mode is significantly split off inside the gap, cf. inset.

Figure 2. The black line shows the amplitude correlation function Equation (23) times [Forumla omitted. See PDF.]. The intersection of the horizontal lines with the black curve yields the position of the Higgs pole in Equation (22). In the BCS + RPA approximation, quasiparticles interact with an effective matrix element [Forumla omitted. See PDF.] [Equation (19)] (red dashed line). This produces a pole in the denominator of Equation (22) at [Forumla omitted. See PDF.]. In the TDGA, it turns out that the effective interaction (exemplified by the blue dashed line) produces a pole within the spectral gap.

Figure 3. The individual terms 1, 2, 3, and 4 which are defined in Equation (36).

Appendix A. BCS Correlation Functions

The correlation functions in the non-interacting BCS limit read (A1)$$\begin{array}{ccc}\hfill {\chi }_q^0(\delta {\Delta }_{-q}^{†},\delta {\Delta }_q^{†})& =& {\chi }_q^0(\delta {\Delta }_q,\delta {\Delta }_{-q})=-{\displaystyle \frac{2}{N}}\sum _k{\alpha }_k^{+}{\beta }_k^{+}{\alpha }_k^{-}{\beta }_k^{-}{\displaystyle \frac{E_k^{+}+E_k^{-}}{{\omega }^2-{(E_k^{+}+E_k^{-})}^2}}\hfill \end{array}$$ (A2)$$\begin{array}{ccc}\hfill {\chi }_q^0(\delta {\Delta }_{-q}^{†},\delta {\Delta }_{-q})& =& {\displaystyle \frac{1}{N}}\sum _k\left\{{\displaystyle \frac{{\left({\beta }_k^{+}\right)}^2{\left({\beta }_k^{-}\right)}^2}{\omega -E_k^{+}-E_k^{-}}}-{\displaystyle \frac{{\left({\alpha }_k^{+}\right)}^2{\left({\alpha }_k^{-}\right)}^2}{\omega +E_k^{+}+E_k^{-}}}\right\}\hfill \end{array}$$ (A3)$$\begin{array}{ccc}\hfill {\chi }_q^0(\delta {\Delta }_q,\delta {\Delta }_q^{†})& =& {\displaystyle \frac{1}{N}}\sum _k\left\{{\displaystyle \frac{{\left({\alpha }_k^{+}\right)}^2{\left({\alpha }_k^{-}\right)}^2}{\omega -E_k^{+}-E_k^{-}}}-{\displaystyle \frac{{\left({\beta }_k^{+}\right)}^2{\left({\beta }_k^{-}\right)}^2}{\omega +E_k^{+}+E_k^{-}}}\right\}\hfill \end{array}$$ (A4)$$\begin{array}{ccc}\hfill {\chi }_q^0(\delta {\rho }_q,\delta {\rho }_{-q})& =& {\displaystyle \frac{2}{N}}\sum _k{\left({\alpha }_k^{+}{\beta }_k^{-}+{\alpha }_k^{-}{\beta }_k^{+}\right)}^2{\displaystyle \frac{E_k^{+}+E_k^{-}}{{\omega }^2-{(E_k^{+}+E_k^{-})}^2}}\hfill \end{array}$$ (A5)$$\begin{array}{ccc}\hfill {\chi }_q^0(\delta {\Delta }_{-q}^{†},\delta {\rho }_{-q})& =& {\chi }_q^0(\delta {\rho }_q,\delta {\Delta }_{-q})\hfill \\ & =& {\displaystyle \frac{1}{N}}\sum _k\left({\alpha }_k^{+}{\beta }_k^{-}+{\alpha }_k^{-}{\beta }_k^{+}\right)\left\{{\displaystyle \frac{{\beta }_k^{+}{\beta }_k^{-}}{\omega -E_k^{+}-E_k^{-}}}+{\displaystyle \frac{{\alpha }_k^{+}{\alpha }_k^{-}}{\omega +E_k^{+}+E_k^{-}}}\right\}\hfill \end{array}$$ (A6)$$\begin{array}{ccc}\hfill {\chi }_q^0(\delta {\Delta }_q,\delta {\rho }_{-q})& =& {\chi }_q^0(\delta {\rho }_q,\delta {\Delta }_q^{†})\hfill \\ & =& -{\displaystyle \frac{1}{N}}\sum _k\left({\alpha }_k^{+}{\beta }_k^{-}+{\alpha }_k^{-}{\beta }_k^{+}\right)\left\{{\displaystyle \frac{{\alpha }_k^{+}{\alpha }_k^{-}}{\omega -E_k^{+}-E_k^{-}}}+{\displaystyle \frac{{\beta }_k^{+}{\beta }_k^{-}}{\omega +E_k^{+}+E_k^{-}}}\right\}\hfill \end{array}$$

References

1. Weinberg, S. The Quantum Theory of Fields—Vol. 2: Modern Applications; Cambridge University Press: Cambridge, UK, 1996.

2. Higgs, P.W. Broken symmetries, massless particles and gauge fields. Phys. Lett.; 1964; 12, pp. 132-133. [DOI: https://dx.doi.org/10.1016/0031-9163(64)91136-9]

3. Nagaosa, N. Quantum Field Theory in Condensed Matter Physics; Springer: Berlin/Heidelberg, Germany, 1999.

4. Pekker, D.; Varma, C.M. Amplitude/Higgs Modes in Condensed Matter Physics. Annu. Rev. Condens. Matter Phys.; 2015; 6, pp. 269-297. [DOI: https://dx.doi.org/10.1146/annurev-conmatphys-031214-014350]

5. Anderson, P.W. Coherent Excited States in the Theory of Superconductivity: Gauge Invariance and the Meissner Effect. Phys. Rev.; 1958; 110, pp. 827-835. [DOI: https://dx.doi.org/10.1103/PhysRev.110.827]

6. Papenkort, T.; Axt, V.M.; Kuhn, T. Coherent dynamics and pump-probe spectra of BCS superconductors. Phys. Rev. B; 2007; 76, 224522. [DOI: https://dx.doi.org/10.1103/PhysRevB.76.224522]

7. Matsunaga, R.; Hamada, Y.I.; Makise, K.; Uzawa, Y.; Terai, H.; Wang, Z.; Shimano, R. Higgs Amplitude Mode in the BCS Superconductors Nb_1-xTi_xN Induced by Terahertz Pulse Excitation. Phys. Rev. Lett.; 2013; 111, 057002. [DOI: https://dx.doi.org/10.1103/PhysRevLett.111.057002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23952432]

8. Matsunaga, R.; Tsuji, N.; Fujita, H.; Sugioka, A.; Makise, K.; Uzawa, Y.; Terai, H.; Wang, Z.; Aoki, H.; Shimano, R. Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor. Science; 2014; 345, pp. 1145-1149. [DOI: https://dx.doi.org/10.1126/science.1254697]

9. Mansart, B.; Lorenzana, J.; Mann, A.; Odeh, A.; Scarongella, M.; Chergui, M.; Carbone, F. Coupling of a high-energy excitation to superconducting quasiparticles in a cuprate from coherent charge fluctuation spectroscopy. Proc. Natl. Acad. Sci. USA; 2013; 110, pp. 4539-4544. [DOI: https://dx.doi.org/10.1073/pnas.1218742110]

10. Krull, H.; Bittner, N.; Uhrig, G.S.; Manske, D.; Schnyder, A.P. Coupling of Higgs and Leggett modes in non-equilibrium superconductors. Nat. Commun.; 2016; 7, 11921. [DOI: https://dx.doi.org/10.1038/ncomms11921]

11. Shimano, R.; Tsuji, N. Higgs Mode in Superconductors. Annu. Rev. Condens. Matter Phys.; 2020; 11, pp. 103-124. [DOI: https://dx.doi.org/10.1146/annurev-conmatphys-031119-050813]

12. Kemper, A.F.; Sentef, M.A.; Moritz, B.; Freericks, J.K.; Devereaux, T.P. Direct observation of Higgs mode oscillations in the pump-probe photoemission spectra of electron-phonon mediated superconductors. Phys. Rev. B; 2015; 92, 224517. [DOI: https://dx.doi.org/10.1103/PhysRevB.92.224517]

13. Doniach, S.; Inui, M. Long-range Coulomb interactions and the onset of superconductivity in the high-Tc materials. Phys. Rev. B; 1990; 41, 6668. [DOI: https://dx.doi.org/10.1103/PhysRevB.41.6668] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9992916]

14. Damle, K.; Sachdev, S. Nonzero-temperature transport near quantum critical points. Phys. Rev. B; 1997; 56, 8714. [DOI: https://dx.doi.org/10.1103/PhysRevB.56.8714]

15. Podolsky, D.; Auerbach, A.; Arovas, D.P. Visibility of the amplitude (Higgs) mode in condensed matter. Phys. Rev. B; 2011; 84, 174522. [DOI: https://dx.doi.org/10.1103/PhysRevB.84.174522]

16. Sherman, D.; Pracht, U.S.; Gorshunov, B.; Poran, S.; Jesudasan, J.; Chand, M.; Raychaudhuri, P.; Swanson, M.; Trivedi, N.; Auerbach, A. et al. The Higgs mode in disordered superconductors close to a quantum phase transition. Nat. Phys.; 2015; 11, 188. [DOI: https://dx.doi.org/10.1038/nphys3227]

17. Cea, T.; Castellani, C.; Seibold, G.; Benfatto, L. Nonrelativistic Dynamics of the Amplitude (Higgs) Mode in Superconductors. Phys. Rev. Lett.; 2015; 115, 157002. [DOI: https://dx.doi.org/10.1103/PhysRevLett.115.157002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26550746]

18. Cheng, B.; Wu, L.; Laurita, N.J.; Singh, H.; Chand, M.; Raychaudhuri, P.; Armitage, N.P. Anomalous gap-edge dissipation in disordered superconductors on the brink of localization. Phys. Rev. B; 2016; 93, 180511. [DOI: https://dx.doi.org/10.1103/PhysRevB.93.180511]

19. Lorenzana, J.; Seibold, G. Long-Lived Higgs Modes in Strongly Correlated Condensates. Phys. Rev. Lett.; 2024; 132, 026501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.132.026501]

20. Cabrera, C.R.; Henke, R.; Broers, L.; Skulte, J.; Collado, H.P.O.; Biss, H.; Mathey, L.; Moritz, H. Effect of strong confinement on the order parameter dynamics in fermionic superfluids. arXiv; 2024; arXiv: 2407.12645

21. Seibold, G.; Lorenzana, J. Time-Dependent Gutzwiller Approximation for the Hubbard Model. Phys. Rev. Lett.; 2001; 86, 2605. [DOI: https://dx.doi.org/10.1103/PhysRevLett.86.2605]

22. Seibold, G.; Becca, F.; Lorenzana, J. Inhomogeneous Gutzwiller approximation with random phase fluctuations for the Hubbard model. Phys. Rev. B; 2003; 67, 085108. [DOI: https://dx.doi.org/10.1103/PhysRevB.67.085108]

23. Seibold, G.; Becca, F.; Rubin, P.; Lorenzana, J. Time-dependent Gutzwiller theory of magnetic excitations in the Hubbard model. Phys. Rev. B; 2004; 69, 155113. [DOI: https://dx.doi.org/10.1103/PhysRevB.69.155113]

24. Seibold, G.; Becca, F.; Lorenzana, J. Theory of Antibound States in Partially Filled Narrow Band Systems. Phys. Rev. Lett.; 2008; 100, 016405. [DOI: https://dx.doi.org/10.1103/PhysRevLett.100.016405]

25. Seibold, G.; Becca, F.; Lorenzana, J. Time-dependent Gutzwiller theory of pairing fluctuations in the Hubbard model. Phys. Rev. B; 2008; 78, 045114. [DOI: https://dx.doi.org/10.1103/PhysRevB.78.045114]

26. Ugenti, S.; Cini, M.; Seibold, G.; Lorenzana, J.; Perfetto, E.; Stefanucci, G. Particle-particle response function as a probe for electronic correlations in the p-d Hubbard model. Phys. Rev. B; 2010; 82, 075137. [DOI: https://dx.doi.org/10.1103/PhysRevB.82.075137]

27. Schiró, M.; Fabrizio, M. Time-Dependent Mean Field Theory for Quench Dynamics in Correlated Electron Systems. Phys. Rev. Lett.; 2010; 105, 076401. [DOI: https://dx.doi.org/10.1103/PhysRevLett.105.076401]

28. Schiró, M.; Fabrizio, M. Quantum quenches in the Hubbard model: Time-dependent mean-field theory and the role of quantum fluctuations. Phys. Rev. B; 2011; 83, 165105. [DOI: https://dx.doi.org/10.1103/PhysRevB.83.165105]

29. Bünemann, J.; Capone, M.; Lorenzana, J.; Seibold, G. Linear-response dynamics from the time-dependent Gutzwiller approximation. New J. Phys.; 2013; 15, 053050. [DOI: https://dx.doi.org/10.1088/1367-2630/15/5/053050]

30. Schrieffer, J.R.; Wen, X.G.; Zhang, S.C. Dynamic spin fluctuations and the bag mechanism of high-Tc superconductivity. Phys. Rev. B; 1989; 16, 11663. [DOI: https://dx.doi.org/10.1103/PhysRevB.39.11663]

31. Vilk, Y.M.; Tremblay, A.M.S. Non-perturbative many-body approach to the Hubbard model and single-particle pseudogap. J. Phys. I; 1997; 7, pp. 1309-1368. [DOI: https://dx.doi.org/10.1051/jp1:1997135]

32. Seibold, G.; Lorenzana, J. Nonequilibrium dynamics from BCS to the bosonic limit. Phys. Rev. B; 2020; 102, 144502. [DOI: https://dx.doi.org/10.1103/PhysRevB.102.144502]

33. Verresen, R.; Moessner, R.; Pollmann, F. Avoided quasiparticle decay from strong quantum interactions. Nat. Phys.; 2019; 15, pp. 750-753. [DOI: https://dx.doi.org/10.1038/s41567-019-0535-3]