Consequences of integrability breaking in quench dynamics of pairing Hamiltonians

We study the collisionless dynamics of two classes of nonintegrable pairing models. One is a Bardeen-Cooper-Schrieffer model with separable energy-dependent interactions, and the other is a two-dimensional topological superconductor with spin-orbit coupling and a band-splitting external field. The long-time quantum quench dynamics at integrable points of these models are well understood. Namely the squared magnitude of the time-dependent order parameter $\mathrm{\Delta }\left(t\right)$ can vanish (Phase I), reach a nonzero constant (Phase II), or periodically oscillate as an elliptic function (Phase III). We demonstrate that nonintegrable models, too, exhibit some or all of these nonequilibrium phases. Remarkably, elliptic periodic oscillations persist, even though both their amplitude and functional form change drastically with integrability breaking. Striking new phenomena accompany loss of integrability. First, an extremely long timescale emerges in the relaxation to Phase III, such that short-time numerical simulations risk erroneously classifying the asymptotic state. This timescale diverges near integrable points. Second, an entirely new Phase IV of quasiperiodic oscillations of $\left|\mathrm{\Delta }\right|$ emerges in the quantum quench phase diagrams of nonintegrable pairing models. As integrability techniques do not apply for the models we study, we develop the concept of asymptotic self-consistency and a linear stability analysis of the asymptotic phases. With the help of these new tools, we determine the phase boundaries, characterize the asymptotic state, and clarify the physical meaning of the quantum quench phase diagrams of BCS superconductors. We also propose an explanation of these diagrams in terms of bifurcation theory.

Figure 1

Illustration of the large timescale $\tau$ that emerges in Phase III quenches $g_i\to g_f$ of nonintegrable pairing models. In all plots, the equilibrium gap corresponding to the initial coupling $g_i$ is ${\mathrm{\Delta }}_{0i}=1.33\times {10}^{-3}W$ while that for the final coupling $g_f$ is ${\mathrm{\Delta }}_{0f}=0.4W,$ and we took $N=2\times {10}^5$ equally spaced single-particle energy levels on the interval $[-W/2,W/2]$. The lines in the plots on the right are the local minima and maxima of the oscillations. In terms of the single-particle level spacing $\delta$, the evolution in the right column goes out to $t_{max}=0.94{\delta }^{-1}$. In (a) and (b), we see that the persistent elliptic oscillations in the integrable $s$-wave case stabilize after a small number of oscillations. In (c) and (d), the amplitude of the oscillations takes roughly 1000 times longer to stop changing. In (e) and (f), integrability is strongly broken and it is not even clear whether the oscillations stabilize to a constant amplitude. The nonintegrable model used was the separable BCS model (2.9) with $f\left(\varepsilon \right)$ from Eq. (5.1). The nearly integrable version uses $\gamma =W,$ while the far from integrable one has $\gamma =1.33\times {10}^{-2}W.$

Figure 2

Ground-state order parameter ${\mathrm{\Delta }}_0$, chemical potential $\mu$, and $E_{\mathrm{gap}}=E_{k=0,+}={(\sqrt{{\mathrm{\Delta }}_0^2+{\mu }^2}-h)}^2$ as functions of the external field $h$ in the spin-orbit model. One simultaneously solves the fermion number equation (2.16) and the self-consistency relationship Eq. (4.5) with the ground-state configuration (4.4). The vanishing of $E_{\mathrm{gap}}$ corresponds to a topological quantum phase transition. The number of fermions is $N_f=0.65N$, where $N$ is the number of spins. We express energies in units of the bandwidth $W$, including the spin-orbit coupling ${\alpha }^2=0.1W$, the level spacing $\delta =W/(N-1)$, and the BCS coupling $g=0.9\delta$. The Fermi energy in these units is ${\varepsilon }_F=\frac{W}{2N}{N}_f=0.325W$. These spin-orbit model parameters remain the same for the remainder of this work, up to adjusting the value of $N$. We do not consider a similar plot for the separable BCS model because in the particle-hole symmetric case considered, the fermion number $N_f=N$ and thus $\mu =0$.

Figure 3

Examples of Phase I quenches for separable BCS models. The equilibrium gaps ${\mathrm{\Delta }}_{0i}, {\mathrm{\Delta }}_{0f}$ and integrability breaking parameter $\gamma$ are given in units of the bandwidth $W$, and there are $N=5\times {10}^4$ (a) and $N=2\times {10}^5$ (b) spins. The initial rapid decay of $\mathrm{\Delta }$ is shown, but out of caution one must simulate to longer times (still smaller than the inverse level spacing) in order to verify that the phase is indeed stable.

Figure 4

Quenches in the spin-orbit model that lead to (a) Phase I and (b) Phase II. Here the number of single-particle energies is $N={10}^4$, and all other parameters are the same as given in the caption of Fig. 2.

Figure 5

Examples of Phase II quenches for separable BCS models. In (a) $N=2\times {10}^5$ spins, and the quench from ${\mathrm{\Delta }}_{0i}=0.15W$ is close to the Phase I-II boundary. In (b), $N=5\times {10}^4$. The oscillatory power-law decay to a constant value takes a rather long time, and we have verified out to $t\delta =2$ in (a) and $t\delta =0.5$ in (b) that the amplitude of the oscillations is indeed decreasing to zero with power-law decay. In both plots, ${\mathrm{\Delta }}_0$ and $\gamma$ are expressed in units of the bandwidth.

Figure 6

The quench in the Lorentzian separable BCS model (blue dots) from Figs. 1 and 1 $[\gamma =W]$ and the corresponding elliptic function fit (solid red) from Eq. (5.24) with $a\approx 0.868205, {\mathrm{\Delta }}_{+}\approx 0.941415, {\mathrm{\Delta }}_{-}\approx 0.501511, {\tilde{\mathrm{\Delta }}}_{+}={\tilde{\mathrm{\Delta }}}_{-}^{*}\approx 0.915740+0.002407i$, and $t_0=2.801929$. To obtain these parameters, we fit $\dot{\mathrm{\Delta }}$ to Eq. (5.22) and then shift by the appropriate $t_0$. If a fifth-order polynomial is used instead of $P_4\left[\mathrm{\Delta }\left(t\right)\right]$, the coefficient of the ${\mathrm{\Delta }}^5$ term is $-6.08\times {10}^{-5}$, providing further evidence that this asymptotic $\mathrm{\Delta }\left(t\right)$ is indeed an elliptic function. Although only a short time frame is shown, this fit works well for the entire time interval from $t{\mathrm{\Delta }}_{0f}={10}^4$, which is the timescale after which the oscillation amplitude stabilizes, to the times shown. In this fitting procedure, $\mathrm{\Delta }$ is given in units of ${\mathrm{\Delta }}_{0f}=0.4W$ and time is measured in units of ${\mathrm{\Delta }}_{0f}^{-1}$ as pictured. In terms of the level spacing $\delta =5\times {10}^{-6}W$, the time domain pictured is $0.73125<t\delta <0.731688$.

Figure 7

A Phase III quench in a $f\left(\varepsilon \right)=exp[-|\varepsilon |/\gamma ]$ separable BCS model (blue dots), where $\gamma =0.5W, N=2\times {10}^5, {\mathrm{\Delta }}_{0i}=0.04W, {\mathrm{\Delta }}_{0f}=0.8W$. The corresponding elliptic function fit (solid red) from Eq. (5.24) has $a\approx 0.821896, {\mathrm{\Delta }}_{+}\approx 1.075648, {\mathrm{\Delta }}_{-}\approx 0.566069, {\tilde{\mathrm{\Delta }}}_{+}={\tilde{\mathrm{\Delta }}}_{-}^{*}\approx 0.010686+1.327633i$, and $t_0=2.131916$. To obtain these parameters, we fit $\dot{\mathrm{\Delta }}$ to Eq. (5.22) and then shift by the appropriate $t_0$. If a fifth-order polynomial is used instead of $P_4\left[\mathrm{\Delta }\right]$, the coefficient of the ${\mathrm{\Delta }}^5$ term is $4.22\times {10}^{-9}$. In this fitting procedure, $\mathrm{\Delta }$ is given in units of ${\mathrm{\Delta }}_{0f}$ and time is measured in units of ${\mathrm{\Delta }}_{0f}^{-1}$ as pictured. In terms of the level spacing $\delta$, the time domain pictured is $0.405<t\delta <0.40525$.

Figure 8

A Phase III quench in the spin-orbit model (blue dots), where in units of the bandwidth $W$: ${\varepsilon }_F=0.1, {\alpha }^2=0.9, gN=2.315999, h_i=1.998980, h_f=0.801020$. These parameters uniquely determine the initial and final equilibrium gaps and chemical potentials through the use of Eq. (2.16) and Eq. (4.5). The energies ${\varepsilon }_j$ are uniformly distributed in the interval $[0,W]$, and the number of pseudospins is $N=8\times {10}^4$. As particle-hole symmetry does not hold, we fit $\mathrm{\Omega }\equiv {\left|\mathrm{\Delta }\right|}^2$ to the elliptic function definition in Eq. (5.24). The fit is $a\approx 0.776633, {\mathrm{\Delta }}_{+}\approx 0.096608, {\mathrm{\Delta }}_{-}\approx 0.080316, {\tilde{\mathrm{\Delta }}}_{+}={\tilde{\mathrm{\Delta }}}_{-}^{*}\approx 0.873604+0.883872i$, and $t_0=3.033272$. The fit (solid red) is good for all $t>\tau$, where $\tau$ is the relaxation time defined in Sec. 5c2. Here $\tau {\mathrm{\Delta }}_{0f}\approx 3050$. In the fitting procedure, $\mathrm{\Delta }$ is given in units of ${\mathrm{\Delta }}_{0f}$ and time is measured in units of ${\mathrm{\Delta }}_{0f}^{-1}$ as pictured. In terms of the level spacing $\delta$, the time domain pictured is $1.472<t\delta <1.473$, shortly after which finite-size effects take over.

Figure 9

Example of a deceitful quench in the $f\left(\varepsilon \right)=\varepsilon$ ($d+id$) separable BCS model, which at short times seems to enter Phase III on a similar timescale as the corresponding integrable $s$-wave quench with the same parameters. Part (b) shows that minimum of the $d+id\left|\mathrm{\Delta }\right(t\left)\right|$ is actually evolving over the entire timescale considered, and it is not clear what the asymptotic phase is. For both models, we used $4\times {10}^4$ single-particle energies ${\varepsilon }_j$ uniformly distributed on the interval $[0,W], {\mathrm{\Delta }}_{0f}=0.00625W, {\mathrm{\Delta }}_{0i}=0.05{\mathrm{\Delta }}_{0f}, {\varepsilon }_F=0.25W$ [77]. In Fig. 10, we explore similar quenches in the $d+id$ model at larger energy scales, where the dynamics are faster.

Figure 10

Study of the long time dynamics of $d+id$ model quenches, continued from Fig. 9. We keep the same parameters and the same ratio ${\mathrm{\Delta }}_{0i}/{\mathrm{\Delta }}_{0f}=0.05$ while varying ${\mathrm{\Delta }}_{0f}$. Pictured are the maxima and minima of oscillations of $\left|\mathrm{\Delta }\right|$. Part (a) shows that below a certain critical ${\mathrm{\Delta }}_{0f}\sim 0.0845W$, the amplitude of $\left|\mathrm{\Delta }\right|$ oscillations evolves over an extremely long timescale. When ${\mathrm{\Delta }}_{0f}=0.05W$, there are also multiple incommensurate frequencies, and it is unclear whether the asymptotic state is Phase II, III, or something else entirely. When ${\mathrm{\Delta }}_{0f}=0.075W$, the decay in amplitude of $\left|\mathrm{\Delta }\right|$ resembles typical decays to Phase II seen in other models (see Fig. 5). At ${\mathrm{\Delta }}_{0f}=0.1W$, the system rapidly enters Phase II at a smaller ${\mathrm{\Delta }}_{\infty }$ than would be inferred from the other two cases, indicating that we have crossed a transition point. Part (b) shows a quench at this transition point, where the Phase II state seen for ${\mathrm{\Delta }}_{0f}=0.1W$ exhibits an exponential instability and moves to an oscillatory state with unknown asymptotic behavior. The integrable $s$-wave BCS model, $f\left(\varepsilon \right)=1$, is deep in Phase III for all these values of ${\mathrm{\Delta }}_{0f}$ and ${\mathrm{\Delta }}_{0i}$.

Figure 11

Nonintegrable pairing models exhibit an extremely long relaxation time $\tau$ when the asymptotic state is Phase III, which is most prominent in the evolution of the minima of the oscillations of $\mathrm{\Delta }\left(t\right)$. Pictured is a study of $\tau$ as a function of ${\mathrm{\Delta }}_{0f}$, at fixed ${\mathrm{\Delta }}_{0i}={10}^{-3}W$ [(a) and (b)], and $\tau$ as a function of ${\mathrm{\Delta }}_{0i}$ at fixed ${\mathrm{\Delta }}_{0f}=0.4W$ [(c) and (d)] in the Lorentzian model at $\gamma =0.8W$ in the particle-hole symmetric case. The time $\tau$ is not monotonic in either case, but it is generally a decreasing function of the initial and final coupling strengths $g_i$ and $g_f$. In all plots, ${\mathrm{\Delta }}_{0i}$ and ${\mathrm{\Delta }}_{0f}$ are given in units of the bandwidth $W.$ In (a) and (b) $2.4\times {10}^4>N>1.2\times {10}^4$ and in (c) and (d) $N=8400$.

Figure 12

Study of the relaxation time $\tau$, see Fig. 11, in the Lorentzian model as a function of the integrability breaking parameter $\gamma$ at fixed ${\mathrm{\Delta }}_{0i}=.005W$ and ${\mathrm{\Delta }}_{0f}=0.6W$, where $\gamma =\infty$ is the integrable $s$-wave model. For these quench parameters, both the Lorentzian and $s$-wave models enter Phase III. Parts (a)–(c) show how the minimum of $\mathrm{\Delta }\left(t\right)$ slowly evolves and reaches an asymptote, while part (d) gives $\tau$ near $\gamma =0.4W$, where the minimum satisfies ${\tau }_{min}{\mathrm{\Delta }}_{0f}\approx 89$. This minimum is still greater than the relaxation time of the $s$-wave case, where $\tau {\mathrm{\Delta }}_{0f}\approx 65$. The relaxation time increases sharply away from $\gamma =0.4W$, especially in the direction of decreasing gamma, where $\tau {\mathrm{\Delta }}_{0f}\approx 64500$ at $\gamma =0.11W$. In all plots, $\gamma$ is given in units of the bandwidth $W$ and $N=5500$.

Figure 13

Study of the relaxation time $\tau$ in the Lorentzian model as a function of the integrability breaking parameter $\gamma$ at fixed ${\mathrm{\Delta }}_{0i}=0.2W$ and ${\mathrm{\Delta }}_{0f}=0.6W$. For these quench parameters, the $s$-wave model enters Phase II, while the Lorentzian enters Phase III. Panel (a) shows how the minimum of $\mathrm{\Delta }\left(t\right)$ slowly evolves and reaches an asymptote, while part (b) gives ${\tau }_{\min }$ near $\gamma =0.6W$, where ${\tau }_{min}{\mathrm{\Delta }}_{0f}\approx 175$. The relaxation time increases away from $\gamma =0.6W$ in both directions. In all plots, $\gamma$ is given in units of the bandwidth $W$ and $N=2800$.

Figure 14

(a) Examples of exactly self-consistent, periodic $\mathrm{\Delta }\left(t\right)$'s for the Lorentzian separable BCS equations of motion for different values of $\gamma$ at fixed ${\mathrm{\Delta }}_{0f}=0.5W$, period $T=225/W$, and $N=500$. For these fixed parameters, below ${\gamma }_{\min }\sim 0.172W$ the only exactly self-consistent, periodic $\mathrm{\Delta }\left(t\right)$ is a constant in time equal to the equilibrium value. (b) Convergence of ${\gamma }_{\min }$ as a function $N$. In both plots, ${\mathrm{\Delta }}_{0f}$ and $\gamma$ are given in units of $W$ and $T$ in units of $W^{-1}$.

Figure 15

Evidence that the self-consistent periodic $\mathrm{\Delta }\left(t\right)$ from Fig. 14 are elliptic functions. Squared time derivatives $\dot{\mathrm{\Delta }}\left(t\right)$ as a function of $\mathrm{\Delta }\left(t\right)$ are given by solid blue lines. These lines overlap strongly with the dashed lines, which are the fits to the defining differential equation for elliptic functions Eq. (5.22). If a ${\mathrm{\Delta }}^5$ coefficient is included in the fits, then it is several orders of magnitude smaller than those for the fourth-order fit shown here, providing strong evidence that ${\dot{\mathrm{\Delta }}}^2\left(t\right)$ is indeed a fourth-order polynomial in $\mathrm{\Delta }\left(t\right)$. In this plot, $\gamma$ is given in units of $W$ and $\mathrm{\Delta }$ in units of ${\mathrm{\Delta }}_{0f}$.

Figure 16

Quenches of nonintegrable separable BCS and spin-orbit models that do not conform to Phases I, II, or III. This quasiperiodic dynamics of the order parameter emerge early and persist for the entire time of the simulation, see also Fig. 17. Panel (a) is the particle-hole symmetric separable BCS model with sine coupling from Eq. (5.2) and $N=4\times {10}^5$ spins. In units of the bandwidth, the integrability breaking parameter is $\gamma =0.075$, while ${\mathrm{\Delta }}_{0i}=0.05$ and ${\mathrm{\Delta }}_{0f}=0.5$. Panel (b) is the spin-orbit model with $N=2\times {10}^5$ spins. In units of the bandwidth: ${\varepsilon }_F=0.4, {\alpha }^2=0.4, gN=2, h_i=2$, and $h_f=0.514256$.

Figure 17

Darker blue points are local minima and maxima of the oscillations for the quenches from Fig. 16 for the entire time of the simulations. These plots suggest that there are regions of quasiperiodicity (Phase IV) in the quantum quench phase diagrams of nonintegrable pairing models. Panel (a) is the same quench as in Fig. 16, and panel (b) corresponds to Fig. 16. In terms of the inverse level spacing, the time evolution goes out to $t_{max}=0.625{\delta }^{-1}$ in (a) and to $t_{max}={\delta }^{-1}$ in (b).

Figure 18

Phase diagram for interaction quenches $g_i\to g_f$ in the integrable limit $h_f=h_i=0$ of the spin-orbit model. Apart from the varying coupling constant $g$, the model parameters are the same as found in Fig. 2. The black dotted lines ${\mathrm{\Delta }}_{0i}=e^{\pm \pi /2}{\mathrm{\Delta }}_{0f}$ indicate the weak-coupling limit ($\mathrm{\Delta }\ll W$) phase boundaries [53]. The thick blue lines mark the true phase boundaries, which are characterized by the appearance of a new pair of complex roots of Eq. (B9) when passing from Phase I to Phase II or Phase II to Phase III.

Figure 19

Behavior of the roots of $L^2\left(u\right)$ for quenches from the ground state of $h_i\ne 0$ to $h_f=0$ in the continuum limit with spin-orbit parameters as given in Fig. 2. Each solid line is the absolute value of the imaginary part of a pair of complex conjugate roots. Regions of $h_i$ with one such line indicate that the asymptotic state is Phase II, while the region where there are two separate lines indicate Phase III. The vertical dashed lines indicate various critical values of $h_i$ where the system undergoes a phase transition or crossover. From left to right, $h_1=0.7813{\varepsilon }_F$ is the topological transition of the ground state, $h_2=0.9938{\varepsilon }_F$ is a Phase II-III transition, $h_3=1.6625{\varepsilon }_F$ is the BCS-BEC crossover, and $h_4=2.2938{\varepsilon }_F$ is a Phase III-II transition which also appears to correspond to ${\mathrm{\Delta }}_{0i}=0$ being the only self-consistent initial equilibrium gap. These critical values of $h_i$ depend in general on the various spin-orbit model parameters.