Nonequilibrium dynamics of superconductivity in the attractive Hubbard model

# Nonequilibrium dynamics of superconductivity in the attractive Hubbard model

## Abstract

We present a framework of semiclassical superconductivity (SC) dynamics that properly includes effects of spatial fluctuations for the attractive Hubbard model. We consider both coherent and adiabatic limits. To model the coherent SC dynamics, we develop a real-space von Neumann equation based on the time-dependent Hartree-Fock-Bogoliubov theory. Applying our method to interaction quenches in the negative- Hubbard model, we show that the relaxation of SC order at weak coupling is dominated by Landau-damping. At strong coupling, we find a two-stage relaxation of the pairing field: a collapse of the synchronized oscillation of Cooper pairs due to spatial inhomogeneity, followed by a slow relaxation to a quasi-stationary state. SC dynamics in adiabatic limit is described by a quantum Landau-Lifshitz equation with Ginzburg-Landau relaxation. Numerical simulations of pump-probe excitation show that long time recovery of the pairing field is dominated by defects dynamics. Our results demonstrate the important role of spatial fluctuations in both limits.

The nonequilibrium dynamics of superconductivity (SC) subject to an external stimulation has been intensively studied for some time (1); (2); (6); (7); (3); (4); (5). This interest has recently been renewed by remarkable pump-probe experiments reporting the observation of the collective amplitude mode (8); (9). Prior work largely studies two physical limits determined by the relation between the quasiparticle relaxation time and the relaxation time of SC order parameter. The SC dynamics in the collisionless limit, , can be described using a time-dependent self-consistent field approach. An interesting phenomenon in this regime is the collective Rabi oscillations of the order parameter (2); (3). In the opposite, adiabatic limit, , the dynamics of the pairing field is usually described by time-dependent Ginzburg-Landau (TDGL) equation (10); (11). TDGL has recently been employed to simulate the out-of-equilibrium dynamics of superconductors in pump-probe experiments (12); (13).

Numerical simulations based on TDGL also provide useful insights on the dynamical inhomogeneity of nonequilibrium superconductivity (14); (15), for example, the formation of topological defects by rapid thermal quenches as described in the Kibble-Zurek scenario (16); (17). However, contrary to the numerous large-scale TDGL simulation studies, effects of spatial fluctuations in the collisionless limit is rarely investigated even though the emergence of inhomogeneous states is a rather common feature of systems far from the equilibrium.

In this paper, we present a theoretical framework for SC dynamics that properly includes the effects of spatial fluctuations in both the collisionless and adiabatic limits. To model the collisionless (coherent) limit, we develop a real-space formulation of the time-dependent Hartree-Fock-Bogoliubov (TDHFB) theory, which is particularly suitable for the negative- Hubbard model. Numerical solution shows two distinct dynamical regimes for the evolution of the pairing field. In particular, we demonstrate an inhomogeneity-induced collapsing of the synchronized oscillation of the SC order parameter at large coupling. To model the adiabatic limit, we use a real-space Anderson pseudo-spin representation (18) and show that the SC dynamics is described by a quantum Landau-Lifshitz equation with Ginzburg-Laudau relaxation.

We consider the Hubbard model with an attractive on-site interaction on the square lattice (20); (21); (19); (22):

 H=−∑ij,αtijc†i,αcj,α+U∑ini,↑ni,↓−∑iμni, (1)

Here is the creation operator of fermions with spin at site-, is the fermion number operator, and . We assume on-site attraction, . A nonzero chemical potential tunes the fermion density away from half-filling. The negative- Hubbard model provides a simple platform for investigating the crossover between the BCS and the Bose-Einstein condensation (BEC) regimes of fermionic superfluid. The model is known to have an enhanced O(3) symmetry at half-filling (), corresponding to the co-existence of SC and charge density wave (CDW) orders. Away from , this degeneracy is lifted and the SC state is energetically more favorable (20); (21).

To model the coherent SC dynamics, we develop a real-space equation of motion approach based on the TDHFB theory, which is equivalent to the time-dependent Bogoliubov-de Gennes equations (23). By introducing a pairing field and on-site density , we perform the Hubbard-Stratonovich transformation to obtain the following BdG Hamiltonian:

 HBdG=−∑ij,αtijc†i,αcj,α+U∑i(Δic†i,↑c†i,↓+h.c.) +∑i(Uρi−μ)ni+U∑i(|Δi|2+ρ2i). (2)

Self-consistency requires that and for both spins , where the average is computed using the BdG Hamiltonian in Eq. (Nonequilibrium dynamics of superconductivity in the attractive Hubbard model). This auxiliary field Hamiltonian without the self-consistency constraint is usually the starting point of determinant quantum Monte Carlo (DQMC) method, which can be applied to the negative- Hubbard model thanks to absence of the sign-problem (20); (21); (22). Equation (Nonequilibrium dynamics of superconductivity in the attractive Hubbard model) is also the basis of classical Monte Carlo method assuming static auxiliary fields (24); (25); (26); (27). Although this method neglects quantum fluctuations in imaginary time, it does take into account thermal and spatial fluctuations of the order-parameter field, which are not included in a conventional mean-field treatment. Indeed, results obtained from the static auxiliary-field Monte Carlo agree remarkably well with that of DQMC simulations (26); (27).

The dynamics of the pairing field in the TDHFB theory is given by the Heisenberg equation of motion, . In fact, the on-site density and the SC order are essentially the diagonal elements of the normal and anomalous density matrices, respectively. Namely, and for a non-magnetic superfluid, which is our primary interest. A description of the SC dynamics requires the time evolution of both and , which is governed by the generalized von Neumann equation, , where is a generalized single-particle density matrix that includes normal as well as anomalous components. Explicitly,

 iℏdρijdt = ∑k(tikρkj−ρiktkj)+U(ρi−ρj)ρij (3) +U(ΔiΔ∗ij−ΔijΔ∗j), iℏdΔijdt = ∑k(tikΔkj+Δiktkj)+[U(ρi+ρj)−2μ]Δij (4) +UδijΔi−U(ρijΔj+ρjiΔi).

Similar equations have been used to model the superfluid dynamics of cold atoms (28); (29) using time-dependent density-functional theory (30).

The above TDHFB can be further simplified if we ignore the spatial fluctuations of SC and CDW order parameters; we shall refer to this approximation as the time-dependent mean-field (TDMF) method. The assumption of translation invariance allows for direct solution of Eqs. (3) and (4) in the Fourier representation. Furthermore, the number of dynamical variables reduces from to , where is the number of lattice sites. This TDMF method becomes exact in the special case of a BCS superconductor (2); (3); (4); (5), and has been widely employed to study quenched superconductors (31); (32); (33), including the negative- Hubbard model (34).

Here we apply the real-space TDHFB formulation to simulate quenches of the interaction strength in the negative- Hubbard model at zero temperature. We consider a time-dependent on-site attraction, for , which suddenly increases in magnitude, for . A chemical potential is used to stabilize a small SC order in the initial state. We have also included a tiny random on-site potential of order . This small initial perturbation is introduced to examine whether the out-of-equilibrium states are stable against inhomogeneities. The time dependence of the SC and CDW order parameters shown in Fig. 1(a) exhibits two distinct dynamical behaviors depending on the strength of the final . The quasi-stationary value at large is plotted in Fig. 1(b) as a function of ; also shown for comparison is the equilibrium at the corresponding . The crossover from the weak to strong coupling regimes roughly corresponds to the maximum of . We discuss the characteristic features of the two dynamical regimes below.

For small , such as the case shown in Fig. 1(c), the pairing order parameter oscillates with a frequency that is proportional to and exhibits collisionless dephasing. The dephased oscillation mainly results from the energy exchange between the collective mode, i.e. the SC order parameter , and the individual Cooper pairs in momentum space, as described in the Landau-damping mechanism. Moreover, we find that the pairing field in this weak-coupling regime shows weak inhomogeneity, which is dominated by long-wavelength fluctuations. The scenario described here is similar to that of interaction quenches of BCS superconductors (2); (3); (4); (5). This is also consistent with the fact that the SC order in the weak-coupling regime of the negative- Hubbard model is better described by a BCS-type model with a large coherence length. At longer time scales, the CDW order emerges as the pairing order parameter settles to its quasi-stationary value.

As increases, the pairing order parameter initially exhibits an oscillation of relatively constant amplitude as demonstrated in Fig. 1(d). The spatial distribution of is rather uniform during this initial stage; see Fig. 2(a). This coherent oscillation in the short timescales resembles the phase-locked regime discussed in the quantum quench of BCS superconductors (2); (3); (4); (5). Importantly, we find that the synchronized oscillation of Cooper pairs is unstable against disorder: any small initial inhomogeneity is amplified after the interaction quench, giving rise to collapse of the synchronized oscillation as shown in Fig. 1(d). Furthermore, the collapse of the initial coherent oscillation is accompanied by the emergence of highly inhomogeneous SC and CDW fields; see Figs. 2(b)–(d).

The emergence of the spatial inhomogeneity can also be understood from the fact that the SC state of the negative- Hubbard model is better described by a BEC of preformed fermion pairs in real space. Consequently, Landau-damping of collective SC order through these local Cooper pairs naturally leads to an inhomogeneous state as each oscillates with its own amplitude and frequency.

It is worth noting that the above scenario of collapsed oscillation induced by the emergence of inhomogeneous has been overlooked by previous TDMF studies (4); (5) which assumed a spatially uniform pairing field. This assumption is valid as long as the coherence length of Cooper pairs is larger than the system size. However, in the large regime of the attractive Hubbard model, the coherent length of a BEC pair is of the order of a few lattice constants. The effects of spatial inhomogeneity then have to be properly taken into account, which we achieved in our real-space simulations. Interestingly, the disordering of the local pairing field shown in Figs. 2(b)–(d) is mainly through the formation of domain walls. This is rather surprising as topological defects, or vortices, are expected to be generated in such processes. We believe this is the artifact of the TDHFB method, which neglects quantum fluctuations. Generalizations of our approach that include the dynamics of higher-order correlations or stochastic extension of the TDHFB are required in order to include fluctuations in both real-space and imaginary-time.

Now we turn to the adiabatic limit, , of the TDHFB and show that the pairing order parameter follows a novel Ginzburg-Landau-Lifshitz dynamics. We first employ the well known canonical transformation and , to map Hamiltonian (1) to a positive- Hubbard model at half-filling in a uniform magnetic field  (20); (21). Here and the phase factor describes a checkerboard pattern on the square lattice. In the transformed representation, the order parameters can be conveniently grouped into a pseudo-spin with components:

 Txi+iTyi=ΔieiQ⋅ri,Tzi=(ρi−1)/2. (5)

This essentially maps the combined SC/CDW order into a spin-density wave (SDW) order in a repulsive Hubbard model. The vector is the real-space version of the pseudo-spin introduced by Anderson for the dynamics of BCS superconductors (18). In the large limit, the BdG Hamiltonian reduces to the Heisenberg exchange interaction, and the pseudo-spins satisfy the LL dynamics in real space:  (35).

Generalization of the semiclassical SDW dynamics to the intermediate regime has been recently developed in Ref. (36). Here we briefly outline the formulation in our context. The time evolution of the pseudo-spin is governed by the conservation of pseudo angular momentum:

 dTidt=−i2∑jtijσβα(~ρiα,jβ−~ρjα,iβ), (6)

where are the reduced density matrix elements for the transformed fermions, and are related to and in the original representation. The dynamics of the density matrix is again governed by von Neumann equation , which is equivalent to Eqs. (3) and (4). Here is the single-particle transformed Hamiltonian; symbolically, .

In the adiabatic limit, electrons are assumed to quickly relax to the equilibrium state of the instantaneous . Consequently, we approximate the density matrix on the right-hand side of Eq. (6) by that describes the equilibrium electron liquid of the instantaneous Hamiltonian, i.e. . Substituting in Eq. (6) leads to our quantum Landau-Lifshitz dynamics (QLLD) method (36). We include a Ginzburg-Landau type damping and stochastic driving force to obtain

 Extra open brace or missing close brace (7)

where , denotes the components of pseudo-spin, and are the damping constants of SC and CDW order parameters, respectively, and is a -correlated fluctuating force satisfying and . Note that QLLD requires solution of the equilibrium density matrix at every time-step, in analogy to Born-Oppenheimer quantum molecular dynamics (37). Rather than direct diagonalization of , we use the kernel polynomial method (38); (39) with gradient-based probing (40); (41) to estimate , and thus effective forces, at a cost that scales linearly with system size.

We next apply the QLLD to investigate the ultrafast relaxation of SC order subject to a short laser pulse. For simplicity, we assume that the effect of the pump pulse is to inject energy to the electrons, which quickly equilibrate to a state characterized by temperature . This is consistent with our adiabatic approximation for the SDW dynamics. The time dependence of the effective electron temperature is governed by the rate equation  (42), where is the heat-capacity of the electron liquid, is the coupling to the lattice, is the lattice temperature, and is the heat source due to the pump pulse. We further assume that throughout the relaxation process. The resultant , shown in Fig. 3(a), then controls the magnitude of the stochastic noise in our QLLD simulations of Eq. (7). We use the parameters, , damping , , , and , in units of the NN hopping .

Fig. 3(a) shows the time dependence of the magnitude of the pairing parameter and the electron energy gap obtained from simulations. The energy gap is estimated from the instantaneous electron density of states , shown in Fig. 3(c) for a few representative simulation times. Interestingly, while both quantities exhibit significant drop after the pulse excitation, it seems a very long timescale is required for their recovery even when temperature returns to almost zero. This rather slow dynamics can be attributed to long-lived topological vortex defects of the phase field , as demonstrated in Fig. 4. Here is the phase angle of the SC order parameter, i.e. . The time dependence of the vortex density , shown in Fig. 3(b), displays a long tail after the pulse excitation. It is known that phase-ordering of a system when quenched into a symmetry-breaking phase is dominated by defect dynamics (43). For the model, mean-field analysis within the TDGL framework suggests that pair annihilation of defects follows a power-law (44), with exponent in 2D. Careful numerical analyses have found logarithmic corrections to this power-law scaling (45); (46). Our preliminary analysis of the SC dynamics suggests a power-law tail with , which is consistent with known results.

To summarize, we have demonstrated the importance of spatial fluctuations in the nonequilibrium dynamics of superconductivity, especially for pairing field with short coherence lengths. For dissipationless SC dynamics, we have developed a real-space von Neumann dynamics method within the TDHFB framework. Applying our method to interaction quenches of the negative- Hubbard model, we have demonstrated that the quench-induced synchronized oscillation of Cooper pairs is unstable against inhomogeneity in the large regime. Finally, we have shown that the SC pairing field obeys a Landau-Lifshitz dynamics in the adiabatic limit. By retaining the electron degrees of freedom, large-scale QLLD simulations provide a unique capability to investigate the intriguing interplay of topological defects of the pairing field and the underlying quasiparticles.

###### Acknowledgements.
The authors thank C. Batista, G. Kotliar, H. Suwa, and Z. Wang for collaboration on related works and insightful discussions. G.-W. C. acknowledges support from the Center for Materials Theory as a part of the Computational Materials Science (CMS) program, funded by the DOE Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. K. B. acknowledges support from the Laboratory Directed Research and Development (LDRD) program at LANL.

### References

1. A. F. Volkov and Sh. M. Kogan, Collisionless relaxation of the energy gap in superconductors. Sov. Phys. JETP 38, 1018 (1974).
2. R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Collective Rabi Oscillations and Solitons in a Time-Dependent BCS Pairing Problem, Phys. Rev. Lett. 93, 160401 (2004).
3. E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler, Relaxation and Persistent Oscillations of the Order Parameter in Fermionic Condensates Phys. Rev. Lett. 96, 097005 (2006).
4. R. A. Barankov and L. S. Levitov, Synchronization in the BCS Pairing Dynamics as a Critical Phenomenon, Phys. Rev. Lett. 96, 230403 (2006).
5. E. A. Yuzbashyan and M. Dzero, Dynamical Vanishing of the Order Parameter in a Fermionic Condensate, Phys. Rev. Lett. 96, 230404 (2006).
6. A.V. Andreev, V. Gurarie, and L. Radzihovsky, Nonequilibrium Dynamics and Thermodynamics of a Degenerate Fermi Gas Across a Feshbach Resonance, Phys. Rev. Lett. 93, 130402 (2004).
7. M. H. Szymańska, B. D. Simons, and K. Burnett, Dynamics of the BCS-BEC Crossover in a Degenerate Fermi Gas, Phys. Rev. Lett. 94, 170402 (2005).
8. R. Matsunaga, Y. I. Hamada, K. Makise, Y. Uzawa, H. Terai, Z. Wang, and R. Shimano, Higgs Amplitude Mode in the BCS Superconductors NbTiN Induced by Terahertz Pulse Excitation, Phys. Rev. Lett. 111, 057002 (2013).
9. R. Matsunaga, N. Tsuji, H. Fujita, A. Sugioka, K. Makise, Y. Uzawa, H. Terai, Z. Wang, H. Aoki, and R. Shimano, Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor, Science 345, 1145 (2014).
10. K. Binder, Time-Dependent Ginzburg-Landau Theory of Nonequilibrium Relaxation, Phys. Rev. B 8, 3423 (1973).
11. P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49, 435 (1977).
12. I. Madan, P. Kusar, V. V. Baranov, M. Lu-Dac, V. V. Kabanov, T. Mertelj, and D. Mihailovic, Real-time measurement of the emergence of superconducting order in a high-temperature superconductor, Phys. Rev. B 93, 224520 (2016).
13. D. M. Kennes and A. J. Millis, Electromagnetic response during quench dynamics to the superconducting state: Time-dependent Ginzburg-Landau analysis, Phys. Rev. B 96, 064507 (2017).
14. G. Karra and R. J. Rivers, Reexamination of Quenches in He (and He), Phys. Rev. Lett. 81, 3707 (1998).
15. N. B. Kopnin and E. V. Thuneberg, Time-Dependent Ginzburg-Landau Analysis of Inhomogeneous Normal-Superfluid Transitions, Phys. Rev. Lett. 83, 116 (1999).
16. T. W. B. Kibble, Topology of cosmic domains and strings, J. Phys. A 9, 1387 (1976).
17. W. H. Zurek, Cosmological experiments in condensed matter systems, Phys. Rep. 276, 177 (1996).
18. P. W. Anderson, Random-Phase Approximation in the Theory of Superconductivity, Phys. Rev. 112, 1900 (1958).
19. R. Micnas, J. Ranninger, S. Robaszkiewicz, Superconductivity in narrow-band systems with local nonretarded attractive interactions, Rev. Mod. Phys. 62, 113 (1990).
20. R. T. Scalettar, E. Y. Loh, J. E. Gubernatis, A. Moreo, S. R. White, D. J. Scalapino, R. L. Sugar, and E. Dagotto, Phase Diagram of the Two-Dimensional Negative- Hubbard Model, Phys. Rev. Lett. 62, 1407 (1989).
21. A. Moreo and D. J. Scalapino, Two-Dimensional Negative- Hubbard Model, Phys. Rev. Lett. 66, 946 (1991).
22. T. Paiva, R. Scalettar, M. Randeria, and N. Trivedi, Fermions in 2D Optical Lattices: Temperature and Entropy Scales for Observing Antiferromagnetism and Superfluidity, Phys. Rev. Lett. 104, 066406 (2010).
23. P-G. de Gennes, Superconductivity of Metals and Alloys (Benjamin, New York, 1966).
24. M. Mayr, G. Alvarez, C. Sen, and E. Dagotto, Phase Fluctuations in Strongly Coupled -Wave Superconductors, Phys. Rev. Lett. 94, 217001 (2005).
25. G. J. Conduit and Y. Meir, First-principles calculation of electronic transport in low-dimensional disordered superconductors, Phys. Rev. B 84, 064513 (2011).
26. S. Datta, V. N. Singh, and P. Majumdar, Radio-frequency spectroscopy of the attractive Hubbard model in a trap, Phys. Rev. A 89, 053609 (2014).
27. S. Tarat and P. Majumdar, A real space auxiliary field approach to the BCS-BEC crossover, Eur. Phys. J. B 88, 68 (2015).
28. A. Bulgac, M. M. Forbes, and G. Wlazlowski, Towards quantum turbulence in cold atomic fermionic superfluids, J. Phys. B: At. Mol. Opt. Phys. 50, 014001 (2017).
29. A. Boudjemâa and M. Benarous, Anomalous density for Bose gases at finite temperature, Phys. Rev. A 84, 043633 (2011).
30. A. Bulgac, Time-Dependent Density Functional Theory and the Real-Time Dynamics of Fermi Superfluids, Annu. Rev. Nucl. Part. Sci. 63, 97 (2013).
31. M. S. Foster, V. Gurarie, M. Dzero, and E. A. Yuzbashyan, Quench-Induced Floquet Topological -Wave Superfluids, Phys. Rev. Lett. 113, 076403 (2014).
32. F. Peronaci, M. Schiró, and M. Capone, Transient Dynamics of -Wave Superconductors after a Sudden Excitation, Phys. Rev. Lett. 115, 257001 (2015).
33. M. A. Sentef, A. Tokuno, A. Georges, and C. Kollath, Theory of Laser-Controlled Competing Superconducting and Charge Orders, Phys. Rev. Lett. 118, 087002 (2017).
34. J. Bünemann and G. Seibold, Charge and pairing dynamics in the attractive Hubbard model: Mode coupling and the validity of linear-response theory, Phys. Rev. B 96, 245139 (2017).
35. A. A. Burkov and A. Paramekanti, Stability of Superflow for Ultracold Fermions in Optical Lattices, Phys. Rev. Lett. 100, 255301 (2008).
36. G.-W. Chern, K. Barros, Z. Wang, H. Suwa, C. D. Batista, Semiclassical dynamics of spin density waves, Phys. Rev. B 97, 035120 (2018).
37. D. Marx and J. Hutter, Ab initio molecular dynamics: theory and implementation, in Modern Methods and Algorithms of Quantum Chemistry, 2nd ed., edited by J. Grotendorst (John von Neumann Institute for Computing Julich, Germany, 2000).
38. R. N. Silver and H. Röder, Densities of states of mega-dimensional Hamiltonian matrices, Int. J. Mod. Phys. C 5, 735 (1994).
39. A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).
40. K. Barros and Y. Kato, Efficient Langevin simulation of coupled classical fields and fermions, Phys. Rev. B 88, 235101 (2013).
41. Z. Wang, G.-W. Chern, C. D. Batista, and K. Barros, Gradient-based stochastic estimation of the density matrix, J. Chem. Phys. 148, 094107 (2018).
42. S.I. Anisimov, B.L. Kapeliovich, T.L. Perel’man, Electron emission from metal surfaces exposed to ultrashort laser pulses, Sov. Phys. JETP 39, 375 (1974).
43. A. J. Bray, Theory of phase ordering kinetics, Adv. Phys. 43, 357 (1994).
44. H. Toyoki, Pair annihilation of pointlike topological defects in the ordering process of quenched systems, Phys. Rev. A 42, 911 (1990).
45. B. Yurke, A. N. Pargellis, T. Kovacs, and D. A. Huse, Coarsening dynamics of the model, Phys. Rev. E 47, 1525 (1993).
46. A. J. Bray, A. J. Briant, and D. K. Jervis, Breakdown of Scaling in the Nonequilibrium Critical Dynamics of the Two-Dimensional model, Phys. Rev. Lett. 84, 1503 (2000).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters