Damping of the Collective Amplitude Mode in Superconductors with Strong Electron-Phonon Coupling
We study the effect of strong electron-phonon interactions on the damping of the Higgs amplitude mode in superconductors by means of non-equilibrium dynamical mean-field simulations of the Holstein model. In contrast to the BCS dynamics, we find that the damping of the Higgs mode strongly depends on the temperature, becoming faster as the system approaches the transition temperature. The damping at low temperatures is well described by a power-law, while near the transition temperature the damping shows exponential-like behavior. We explain this crossover in terms of a temperature-dependent quasiparticle lifetime caused by the strong electron-phonon coupling, which smears the superconducting gap edge and makes the relaxation of the Higgs mode into quasiparticles more efficient at elevated temperatures. We also reveal that the phonon dynamics can soften the Higgs mode, which results in a slower damping.
pacs:71.10.Fd, 74.40.Gh, 71.38.-k
The problem of how a superconducting (SC) state evolves in time after an external stimulation has attracted the interest of researchers for a long time. Volkov1974 (); Littlewood1982 (); Barankov2004 (); Yuzbashyan2006a (); Barankov2006 (); Yuzbashyan2006 (); Papenkort2007 (); Papenkort2008 (); Gurarie2009 (); Schnyder2011 (); Matsunaga2013 (); Zachmann2013 (); Matsunaga2014 (); Krull2014 (); Kemper2015 (); Sentef2015 (); Tsuji2015 (); Capone2015 (); Murakami2016 (); Krull2016 () On the theoretical side, static mean-field analyses, which neglect inelastic collisions, have been widely employed to study the coherent dynamics. In the conventional weak-coupling BCS regime, coherent oscillations of the SC order parameter are known to decay with a power law () regardless of temperature Volkov1974 () and even beyond the linear-response regime. Yuzbashyan2006a (); Barankov2006 (); Yuzbashyan2006 (); Papenkort2007 (); Papenkort2008 (); Zachmann2013 () The power-law damping is remarkable, since it indicates a relatively slow decay without a specific relaxation timescale, which is usually observed only in special situations, such as near a critical point. Mean-field analyses have further been applied to various situations to reveal changes in the properties of the coherent oscillations beyond the BCS regime for a bulk system. For instance, in the BEC regime the power law becomes ,Gurarie2009 () while in quasi-1D systems the power strongly depends on the thickness of the system. Zachmann2013 () The effect of a finite quasiparticle lifetime in the weak-coupling regime has been briefly addressed in Ref. Littlewood1982, , but the damping of the coherent oscillations in the correlated regime remains an interesting theoretical issue.
The field has recently been stimulated by the observation of the collective amplitude (“Higgs”) mode Anderson1958 (); Higgs1964 (); Anderson1963 () in conventional superconductors in pump-probe experiments with a strong THz laser pump.Matsunaga2013 (); Matsunaga2014 () While previous experiments had already observed the collective amplitude mode in a coexistence region of superconductivity and charge order with Raman spectroscopy,Sooryakumar1980 (); Littlewood1981 (); Littlewood1982 (); Measson2014 (); Pekker2015 () the results reported in Refs. Matsunaga2013, ; Matsunaga2014, demonstrate the novel possibility of studying the collective mode in ordinary superconductors without coexisting orders. The pump-probe experiments can reveal not only the characteristic frequency of the collective mode but also its damping behavior after the pump. One important finding is that the damping of the coherent oscillations induced by a strong THz laser becomes significantly faster with increasing excitation intensity,Matsunaga2013 () which is difficult to explain with the BCS dynamics without collisions.Yuzbashyan2006a (); Barankov2006 (); Yuzbashyan2006 (); Papenkort2007 (); Papenkort2008 (); Zachmann2013 () In addition, the temperature dependence of the damping is attracting experimental interest.Matsunaga_prep () We also note that the sample used in these experiments, NbN, has a strong electron-phonon (el-ph) coupling.Wolf1985 () It is therefore important to go beyond the BCS dynamics and to clarify the effects of strong el-ph couplings on properties of the amplitude Higgs mode.
Theoretical studies of the coherent oscillations in superconductors with strong el-ph couplings have only appeared recently, Kemper2015 (); Sentef2015 (); Murakami2016 () and many important questions remain to be clarified. In this paper, we consider the effects of strong el-ph couplings on the temperature dependence of the damping behavior of the amplitude Higgs mode, taking into account the inelastic collisions of quasiparticles. Our study makes use of the non-equilibrium dynamical mean-field theory (DMFT), which enables us to simulate the damping behavior during the first several cycles after a pump pulse, as in pump-probe experiments. We find that the finite electron quasiparticle lifetime resulting from strong el-ph couplings leads to a strong temperature dependence of the damping of the amplitude mode. We also reveal that the phonon dynamics can soften the amplitude Higgs mode and extend its lifetime.
Ii Model and Method
We focus on a basic model for el-ph coupled systems, the Holstein model, which is described by the Hamiltonian,
where creates an electron with spin at site , is the electron hopping, with , and the electron chemical potential. creates an Einstein phonon with a bare frequency , and is the el-ph coupling. We assume a semi-elliptic density of states for free electrons, . We take as the unit of energy, and focus on the half-filled case (). In this model, the phonon-mediated electron-electron attraction leads to an s-wave SC state, whose order parameter is with being the total number of sites. We choose the order parameter to be real without loss of generality. In order to study the damping of the Higgs oscillations, i.e., the coherent oscillation of the amplitude of the SC order parameter, we consider a field coupled to the pair potential as a pump, . To be precise, we directly simulate the dynamics after a pump using the non-equilibrium DMFT framework (see below), with a small enough (the linear response regime), and evaluate the dynamical pair susceptibility, Murakami2016 ()
where and is the step function. Let us comment on a few points. (i) Even though a pump in the form of a pair potential field and the measurement of the dynamical pair susceptibility are rather academic tools, they allow us to focus on the amplitude dynamics of the order parameter, which has also been considered in previous investigations of the Higgs amplitude mode.Kulik1981 (); Benfatto2014 (); Benfatto2015 () (ii) The pair potential field is related to more realistic excitations. For example, a modulation of the effective attractive interaction () can, within the BCS picture, be regarded as a pair potential field pulse, since the interaction term is . We can expect the same effect for the modulation of the hopping parameter, since by changing the measure of time, we can map the hopping modulation to an interaction modulation. We also note that the hopping modulation can be realized as a second order effect of the electro-magnetic field Tsuji2011 (); Tsuji2015 (), or by modulation of a certain phonon mode.Sentef2015 (); Subedi2011 (); Forst2011a (). (iii) Our goal here is to evaluate the linear response function, Eq. (2), from a simulation of the time-evolution after a pump. In principle, one can obtain the same quantity by solving the Bethe-Salpeter equation. However, the latter procedure usually involves a numerical analytic continuation, which introduces ambiguities. By calculating the real-time information directly, we can avoid the analytic continuation. (Another way to avoid solving the Bethe-Salpeter equation has recently been proposed in Ref. Tsuji2016, .) Even though we can only access the first several oscillations after the pump with the present approach, this is sufficient, since the pump-probe experiment also measure only the first several cycles after a pump.
The dynamics of the system is simulated using the framework of the non-equilibrium DMFT, Aoki2013 () which becomes exact in the limit of infinite spatial dimensions. In DMFT, the lattice model Eq. (1) is mapped onto a single-site impurity model, whose local Hamiltonian is . The effective bath of the impurity problems is determined in a self-consistent manner such that the local electron Green’s function and the local self-energy of the lattice problem coincide with the impurity Green’s function and the impurity self-energy , respectively. Here is a Nambu spinor, is the time ordering operator on the Kadanoff-Baym contour, and a Pauli matrix, where a quantity with a hat represents a Nambu-Gor’kov matrix. Similarly, the impurity phonon Green’s function, , is identified with the local one in the lattice problem, . Here .
We solve the non-equilibrium effective impurity problem with two types of diagrammatic approximations. The first is the self-consistent Migdal (sMig) approximation. Murakami2015 (); Murakami2016 (); Bauer2011 (); Freericks1994 (); Capone2003 (); Hague2008 (); Leeuwen2015 (); Pavlyukh2016 () Here, the electron self-energy () and phonon self-energy () are expressed as
The sMig approximation neglects vertex corrections for the self-energies, which is justified when the phonon frequency is sufficiently smaller than the electron bandwidth. The dimensionless el-ph coupling is defined as . Since we are interested in the strong-coupling regime, we choose . With both self-energies considered self-consistently, electrons and phonons are renormalized, and their dynamics, including collisions between them, are taken into account. The detailed diagrammatic expression for in this approximation has been presented in Ref. Murakami2016, . In addition to the ladder diagrams with electron legs, which are already taken into account in the BCS dynamics and include the effect of the relaxation of the Higgs mode into quasiparticle excitations, the self-consistent Migdal approximation includes the ladder diagrams with phonon legs and hybridizations between these two types of diagrams.
where is the bare phonon propagator and the dimensionless el-ph coupling is defined as , which we choose in this paper. In this approximation, while electrons are renormalized and their collisions with phonons are considered, the phonons are not renormalized, stay in equilibrium and act as a heat bath. The diagrammatic expression for in the unrenormalized Migdal approximation contains the same type of diagrams as the BCS theory.Murakami2016 () Thus the relaxation into quasiparticle excitations is included in this approximation, but the coupling to the phonon dynamics and possible relaxation channels to phonons are ignored. Neglect of the phonon renormalisation makes the approximation less accurate than the self-consistent Migdal approximation for describing the isolated Holstein model. Murakami2015 () However, the uMig approximation phenomenologically describes a situation where the phonons, which are equilibrated by other degrees of freedom than the focused system, act as a heat bath for the electrons (in an open system, beyond the pure Holstein model description). Considering this, we introduce a finite phonon lifetime in the unrenormalized Migdal approximation. The phonon part is expressed as , where stands for the retarded part. The other components of the Green’s function (lesser, greater etc.) are connected to the retarded part by the equilibrium fluctuation-dissipation theorem. Aoki2013 ()
In order to analyze the decay of the amplitude mode after a pump at , we employ two types of fitting functions,
where and are fitting parameters. We use a least-square fit in the time interval , where is chosen as the first time at which .
Since we want to clarify the difference to the BCS mean-field dynamics, we first recapitulate the dynamical pair susceptibility obtained within the BCS approximation. We assume that the attractive interaction is represented as with , where represents the attractive interaction, is the cutoff energy and is the bare electron dispersion with momentum . is expressed as
where the bubble contribution is the retarded part of . One can prove that approaches linearly in the limit of . Here the SC gap in the BCS approximation is . The explicit expression of is
where , and . From this expression, one can see that behaves in the vicinity of , see Fig. 1(a). Hence the denominator of Eq. (7) becomes zero at , which corresponds to the amplitude Higgs mode. We also note that if there were a pole at , we could interpret the quantity as the rate of the scattering from the collective mode to quasiparticle excitations, since it corresponds to the half-width of the peak in the spectrum.Littlewood1982 () Therefore, Eq. (8) implies that below the SC gap () damping channels are energetically unavailable, while they are available above the gap (). Because of the factor in the second line, the contribution from the energetically available channels becomes small near . This results in a decrease of and going to zero, see Fig. 1(a), which indicates a slower decay than an exponential decay of this mode. However, because of the rapid increase of and , this mode is not undamped. Indeed, the analytic expression for the asymptotic behavior of is proportional to Volkov1974 ()
In Fig. 1, we show along with the fitting function Eq. (5). Regardless of the temperature (both well below and near ), the fitting works very well from the first oscillation and the exponent is . This has been found in previous works for various types of excitations.Volkov1974 (); Yuzbashyan2006 (); Papenkort2008 ()
Now we turn to the results of the unrenormalized Migdal approximation (with equilibrium phonons) in order to grasp the effects of the finite quasiparticle lifetime resulting from strong el-ph couplings. Figure 2 shows obtained within this scheme, which we denote by . In contrast to the BCS dynamics, the damping exhibits a strong dependence on the temperature: When the temperature is much lower than , the damping is well described by a power law but with a different exponent from the BCS value, while the exponential fitting is inadequate, see Fig. 2(a). When the temperature is close to , an exponential fit (Eq. (6)) becomes better than a power-law fit (Eq. (5)), as in Fig. 2(b). In between these two regimes, neither of the two fitting forms can accurately describe the decay. We also note that the decay of the amplitude mode as discussed above is not limited to the present type of the pump protocol. The same damping behavior (not shown) is observed after a hopping modulation, , which mimics the effect of a strong laser. Murakami2016 ()
In Fig. 3(a), we show for each fitting function the coefficient of determination, which is defined for a set of data, , and a fitting function as . Here is the average of , and a value of closer to indicates a better fit. It is evident that a power-law (Eq. (5)) provides the better description than an exponential law (Eq. (6)) at low temperatures, while the opposite is true near . The damping coefficients are depicted in Fig. 3(b). In contrast to the BCS dynamics, the damping shows a significant dependence on temperature, i.e., the damping becomes faster with increasing temperature.
Since in some previous analyses of the damping of the Higgs modeZachmann2013 (); Matsunaga2013 () in Eq. (5) was treated as a fitting parameter, we also consider this case and denote the corresponding fitting function as . We note that when and are large enough, behaves for a finite range of , so that can be regarded as an interpolation between and . The result, shown in Fig. 3, shows that always provides a better than and . At low temperatures the fitted value of stays around zero which means that essentially behaves as , while and increase in a similar manner with increasing temperature near , which is consistent with an exponential behavior. We note that it is hard to provide a physical interpretation for large negative (long time before the pump) near and it would be more natural to consider that the fitting function is more essential there, even though gives the best fit. We also note that in the BCS case always stays near zero regardless of the temperature.
Now we discuss the origin of the different decay behaviors in the BCS dynamics and the Migdal dynamics. In the BCS theory, because of the relation Volkov1974 (); Littlewood1982 () and the fact that the lifetime of a quasiparticle is infinite, the relaxation channel of the amplitude mode () is limited to the quasiparticle excitations just at , see Fig. 4(a). However, as shown in Eq. (8), the contribution from this channel becomes 0 due to the factor . On the other hand, both Migdal approximations can take into account collisions, or the lifetime of the quasiparticles and the incoherent parts in the spectrum. In addition to this, as pointed out in our previous analysis,Murakami2016 (); Tsuji2016 () the Higgs energy sticks to the SC gap () even for strong el-ph couplings. Here we emphasize that the relation between the Higgs energy and the SC gap is not trivial in the strong coupling regime and that it strongly affects the damping behavior of the collective mode. Because of these situations, it becomes possible for the amplitude Higgs mode to decay into excitations from the lower band at various ((1) in Fig. 4(b)), and into excitations to the incoherent parts ((2) in Fig. 4(b)). Since the energetically available relaxation channels are no longer restricted to , one can expect finite contributions from these channels. Now the quasiparticle lifetime decreases with increasing temperature, and the process (2) requires thermally excited quasiparticles above the Fermi energy. Hence these decay processes of the amplitude mode should become more significant closer to and make the decay faster. We note that, in the present Holstein model, there is a phonon window, below which the quasiparticle lifetime becomes very long at low enough temperatures. Therefore, at low enough temperatures the process (1) can practically use only , and this channel should suppressed as in the BCS case. This situation should lead the distinct change of damping laws at low temperatures and temperatures around . In more realistic situations, one may need to consider acoustic phonons. In the presence of such phonons, the phonon window should become less clear and the quasiparticle lifetimes should decrease more quickly with increasing temperature. Hence it is expected that the damping laws at low temperatures and temperatures around becomes less distinct and that the damping tends to be faster and more exponential-like. We also note that the temperature dependence of the damping has been briefly addressed in Ref. Littlewood1982, , where they consider the effect of the quasiparticle lifetime on top of the BCS dynamics and suggest that the Higgs mode becomes overdamped near .
The decrease of the quasiparticle lifetime is indeed evident in the temperature dependence of the electron spectrum, see Fig. 4(c)(d). In Fig. 4(d), we show and the Higgs frequency derived from the fitting for various , where stands for the 11 component in the Nambu-Gor’kov form. The gap edge becomes smeared as we increase , while the Higgs frequency is always located near the edge, i.e., indeed holds. Hence we conclude that near the relaxation of the amplitude mode into quasiparticles becomes efficient due to the strong el-ph coupling. This enhances the damping of the oscillations, which is well described by the exponential fit, Eq. (6). Further support for this picture is obtained from , whose value at increases with increasing temperature, as shown by the green arrows in Fig. 4(d).111We note that changes smoothly. The rapid increase of above comes from . We note that, if the time dependence of the irreducible vertex can be approximately described by a delta function with a renormalized coefficient, a type of equation similar to Eq. 7 is obtained, in which case we can indeed interpret as the efficiency of the relaxation.
Now, we move on to the self-consistent Migdal results, which include effects of the phonon dynamics. In order to single out the effects of the phonon dynamics on the Higgs mode, let us first look at the behavior of calculated with renormalized phonons but without phonon dynamics. Namely, we study the time evolution with the self-energy
where the superscript “eq” indicates that the propagator is the equilibrium one. We denote the pair susceptibility evaluated in this way as , since its diagrammatic expression consists of ladder diagrams with electron legs. Murakami2016 () In the result shown in Fig. 5, the only difference from the uMig approximation is that the phonon propagator is renormalized through the el-ph coupling and depends on the temperature. Murakami2016 () As in the uMig case, at low temperatures a power-law fit, Eq. (5), describes the damping well (see Fig. 5(a)), while around an exponential fit, Eq. (6), becomes better (see Fig. 5(b)). We note that, even when the phonon renormalisation is included, sticks to the SC gap edge, , and the quasiparticle lifetime decreases with increasing temperature, which is the same as in the uMig approximation. In Fig. 6, we display (along with fits) and . As pointed out in Ref. Murakami2016, , there emerges another collective amplitude mode originating from the phonon dynamics. In order to deal with this, we have added a term to the fitting functions (Eqs. (5),(6)). Again, the damping becomes faster with increasing temperature, and one can see a crossover of the damping from power law to exponential law. When we compare and , we notice different behavior between them, which indicates that has a lower frequency 222The difference is a few percent and is valid regardless of the existence of the phonon dynamics., and that the oscillations in are more slowly damped. Indeed, the fittings for yield smaller exponents than those for , compare Fig. 5 and Fig. 6. The softening of can be attributed to the hybridization between the Higgs mode and the amplitude mode originating from phonon oscillations. The softening of the Higgs mode makes it longer-lived due to the suppression of the available relaxation channels to quasiparticles. Even though the possible decay of the Higgs mode into two phonons is also considered within the self-consistent approximation, this channel is energetically suppressed by the reduction of the renormalized single-phonon spectral weight at (phonon anomaly). Allen1997 (); Murakami2016 ()
We finally comment on the relation between uMig and sMig. Firstly, recall that these two methods describe different physical set-ups: in sMig the system is isolated, while in uMig the system is open and the feedback from the phonon dynamics is neglected. Our results show that in both cases, one observes a crossover of the damping law as we vary the temperature. The common origin is the decreased quasiparticle lifetime and subsequent enhancement in the number of the available relaxation channels to quasiparticles with increasing the temperature. On the other hand, we would like to note that the uMig and sMig descriptions can lead to very different dynamics since they approach final states with different temperatures after strong excitations.Murakami2015 () In realistic situations, we may need to include both, the energy dissipation from the system and the feedback from the phonon dynamics. Which description is more appropriate depends on the specific problem.
We have studied the damping of the amplitude Higgs mode in a strongly-coupled phonon-mediated superconductor described by the Holstein model. The non-equilibrium DMFT results show that, in a sharp contrast to the BCS dynamics, the damping exhibits a strong temperature dependence and becomes faster as . Specifically, we have revealed that at low temperatures the damping of the Higgs oscillations is well described by a power law with an exponent distinct from the BCS value, and that near the transition temperature the oscillations tend to decay with an exponential law. In addition, we have shown that the phonon dynamics can soften the Higgs frequency and extend its lifetime.
Our study has focused on the initial several cycles of the coherent oscillations after a pump. How the amplitude mode behaves in the long-time limit and how its damping is related to the behavior extracted here from the initial cycles are questions which should be clarified in the future. It is also important to understand the effects of el-ph couplings beyond the Holstein model, such as couplings with acoustic phonons, as well as the effects of the band structure, impurities and the system size.
Acknowledgements.The authors wish to thank R. Matsunaga and R. Shimano for fruitful discussions and for showing us the experimental result prior to publication. HA is supported by a JSPS KAKENHI (No. 26247057) and ImPACT project (No. 2015-PM12-05-01) from JST, while YM has been supported by JSPS Research Fellowships for Young Scientists. PW acknowledges support from FP7 ERC starting grant No. 278023. NT is supported by JSPS KAKENHI (No. 25800192).
- (1) A. Volkov and S. Kogan, Sov. Phys. JETP 38, 1018 (1974).
- (2) P. B. Littlewood and C. M. Varma, Phys. Rev. B 26, 4883 (1982).
- (3) R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Phys. Rev. Lett. 93, 160401 (2004).
- (4) E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler, Phys. Rev. Lett. 96, 097005 (2006).
- (5) R. A. Barankov and L. S. Levitov, Phys. Rev. Lett. 96, 230403 (2006).
- (6) E. A. Yuzbashyan and M. Dzero, Phys. Rev. Lett. 96, 230404 (2006).
- (7) T. Papenkort, V. M. Axt, and T. Kuhn, Phys. Rev. B 76, 224522 (2007).
- (8) T. Papenkort, T. Kuhn, and V. M. Axt, Phys. Rev. B 78, 132505 (2008).
- (9) V. Gurarie, Phys. Rev. Lett. 103, 075301 (2009).
- (10) A. P. Schnyder, D. Manske, and A. Avella, Phys. Rev. B 84, 214513 (2011).
- (11) R. Matsunaga et al., Phys. Rev. Lett. 111, 057002 (2013).
- (12) M. Zachmann et al., New Journal of Physics 15, 055016 (2013).
- (13) R. Matsunaga et al., Science 345, 1145 (2014).
- (14) H. Krull, D. Manske, G. S. Uhrig, and A. P. Schnyder, Phys. Rev. B 90, 014515 (2014).
- (15) A. F. Kemper et al., Phys. Rev. B 92, 224517 (2015).
- (16) M. A. Sentef, A. F. Kemper, A. Georges, and C. Kollath, Phys. Rev. B 93, 144506 (2016).
- (17) N. Tsuji and H. Aoki, Phys. Rev. B 92, 064508 (2015).
- (18) F. Peronaci, M. Schiró, and M. Capone, Phys. Rev. Lett. 115, 257001 (2015).
- (19) Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Phys. Rev. B 93, 094509 (2016).
- (20) H. Krull et al., Nat Commun 7, 11921 (2016).
- (21) P. W. Anderson, Phys. Rev. 112, 1900 (1958).
- (22) P. W. Higgs, Phys. Rev. Lett. 13, 508 (1964).
- (23) P. W. Anderson, Phys. Rev. 130, 439 (1963).
- (24) R. Sooryakumar and M. V. Klein, Phys. Rev. Lett. 45, 660 (1980).
- (25) P. B. Littlewood and C. M. Varma, Phys. Rev. Lett. 47, 811 (1981).
- (26) M.-A. Méasson et al., Phys. Rev. B 89, 060503 (2014).
- (27) D. Pekker and C. Varma, Annual Review of Condensed Matter Physics 6, 269 (2015).
- (28) R. Matsunaga and R. Shimano, private communication.
- (29) K. E. Kihlstrom, R. W. Simon, and S. A. Wolf, Phys. Rev. B 32, 1843 (1985).
- (30) I. Kulik, O. Entin-Wohlman, and R. Orbach, J. Low Temp. Phys. 43, 591 (1981).
- (31) T. Cea and L. Benfatto, Phys. Rev. B 90, 224515 (2014).
- (32) T. Cea, C. Castellani, G. Seibold, and L. Benfatto, Phys. Rev. Lett. 115, 157002 (2015).
- (33) N. Tsuji, T. Oka, P. Werner, and H. Aoki, Phys. Rev. Lett. 106, 236401 (2011).
- (34) A. Subedi and L. Boeri, Phys. Rev. B 84, 020508 (2011).
- (35) M. Först et al., Nat. Phys. 7, 854 (2011).
- (36) N. Tsuji, Y. Murakami, and H. Aoki, arXiv:1606.05024 (2016).
- (37) H. Aoki et al., Rev. Mod. Phys. 86, 779 (2014).
- (38) Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Phys. Rev. B 91, 045128 (2015).
- (39) J. Bauer, J. E. Han, and O. Gunnarsson, Phys. Rev. B 84, 184531 (2011).
- (40) J. K. Freericks, Phys. Rev. B 50, 403 (1994).
- (41) M. Capone and S. Ciuchi, Phys. Rev. Lett. 91, 186405 (2003).
- (42) J. Hague and N. dâAmbrumenil, J. Low Temp. Phys. 151, 1149 (2008).
- (43) N. Skkinen, Y. Peng, H. Appel, and R. van Leeuwen, The Journal of Chemical Physics 143, (2015).
- (44) M. Schüler, J. Berakdar, and Y. Pavlyukh, Phys. Rev. B 93, 054303 (2016).
- (45) P. B. Allen, V. N. Kostur, N. Takesue, and G. Shirane, Phys. Rev. B 56, 5552 (1997).