Consequences of integrability breaking in quench dynamics of pairing Hamiltonians

TL;DR

This paper investigates how breaking integrability affects the long-time dynamics of pairing models after a quantum quench, revealing new phases and long relaxation times.

Contribution

It introduces the concept of asymptotic self-consistency and stability analysis to study nonintegrable pairing models, identifying new dynamical phases and phenomena.

Findings

Persistence of elliptic oscillations despite integrability breaking

Emergence of a new quasiperiodic phase (Phase IV)

Diverging relaxation time near integrable points

Abstract

We study the collisionless dynamics of two classes of nonintegrable pairing models. One is a BCS model with separable energy-dependent interactions, the other - a 2D 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 $\Delta(t)$ can either 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 time scale emerges in the relaxation to Phase III, such that short-time numerical simulations risk erroneously classifying the asymptotic state. This time scale diverges near integrable points. Second, an entirely new Phase IV of quasiperiodic oscillations of $|\Delta|$ 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.

Equations

$$\begin{split}\hat{H}_{f}&=\sum_{j\lambda}\varepsilon_{j}\hat{c}^{{\dagger}}_{j\lambda}\hat{c}_{j\lambda}-\frac{1}{g}\hat{\Delta}^{{\dagger}}\hat{\Delta},\quad\hat{\Delta}\equiv g\sum_{j}f_{j}\hat{c}_{j\downarrow}\hat{c}_{j\uparrow},\\ \hat{H}_{so}&=\sum_{\mathbf{k}ab}\bigg{[}(\varepsilon_{k}\delta_{ab}-h\sigma^{z}_{ab})+\alpha(k_{y}\sigma^{x}_{ab}-k_{x}\sigma^{y}_{ab})\bigg{]}\hat{c}^{{\dagger}}_{\mathbf{k}a}\hat{c}_{\mathbf{k}b}-\\ &-\frac{1}{g}\hat{\Delta}^{{\dagger}}\hat{\Delta},\quad\hat{\Delta}\equiv g\sum_{\mathbf{k}}\hat{c}_{-\mathbf{k}\downarrow}\hat{c}_{\mathbf{k}\uparrow}.\end{split}$$

$$\begin{split}\hat{H}_{f}&=\sum_{j\lambda=\uparrow\downarrow}\varepsilon_{j}\hat{c}^{{\dagger}}_{j\lambda}\hat{c}_{j\lambda}-\sum_{j}f_{j}\bigg{[}\Delta^{*}\hat{c}_{j\downarrow}\hat{c}_{j\uparrow}+h.c.\bigg{]},\\ \hat{H}_{so}&=\sum_{\mathbf{k}\lambda=\pm}\varepsilon_{k\lambda}\hat{a}^{{\dagger}}_{\mathbf{k}\lambda}\hat{a}_{\mathbf{k}\lambda}-\bigg{(}\frac{\Delta}{2}\sum_{\mathbf{k}\lambda}e^{-i\theta_{\mathbf{k}}}\bigg{[}\lambda\sin\phi_{k}\hat{a}^{{\dagger}}_{\mathbf{k}\lambda}\hat{a}^{{\dagger}}_{-\mathbf{k}\lambda}+\\ &+\cos\phi_{k}\hat{a}^{{\dagger}}_{-\mathbf{k}\lambda}\hat{a}^{{\dagger}}_{\mathbf{k}\bar{\lambda}}\bigg{]}+h.c.\bigg{)}\\ \end{split}$$

$$\begin{array}{rl}& \cos {\varphi }_k=\frac{h}{R_k},\sin {\varphi }_k=\frac{\alpha k}{R_k},\\ & R_k=\sqrt{h^2+{\alpha }^2{k}^2},\\ & {\epsilon }_{k\lambda }={\epsilon }_k-\lambda R_k,\lambda =\pm ,\bar{\lambda }=-\lambda ,\\ & \mathbf{k}=k_x+i{k}_y=k{e}^{i{\theta }_{\mathbf{k}}}.\end{array}$$

$$\begin{split}\Delta&=g\sum_{j}f_{j}\langle\hat{c}_{j\downarrow}\hat{c}_{j\uparrow}\rangle,\\ \Delta&=\frac{g}{2}\sum_{\mathbf{k}\lambda=\pm}e^{i\,\theta_{\mathbf{k}}}\bigg{[}\lambda\sin\phi_{k}\langle\hat{a}_{-\mathbf{k}\lambda}\hat{a}_{\mathbf{k}\lambda}\rangle+\cos\phi_{k}\langle\hat{a}_{\mathbf{k}\lambda}\hat{a}_{-\mathbf{k}\bar{\lambda}}\rangle\bigg{]},\end{split}$$

$${\hat{s}}_j^{-}={\hat{s}}_j^x-i{\hat{s}}_j^y={\hat{c}}_{j\downarrow }{\hat{c}}_{j\uparrow },{\hat{s}}_j^z=\frac{1}{2}({\hat{c}}_{j\uparrow }^{†}{\hat{c}}_{j\uparrow }+{\hat{c}}_{j\downarrow }^{†}{\hat{c}}_{j\downarrow }-1).$$

$$\begin{array}{rl}{\dot{\mathbf{s}}}_j& ={\mathbf{b}}_j\times {\mathbf{s}}_j,{\mathbf{b}}_j=(-2f_j{\Delta }_x,-2f_j{\Delta }_y,2{\epsilon }_j),\end{array}$$

$$\Delta =g\sum _j{f}_j{s}_j^{-}={\Delta }_x-i{\Delta }_y.$$

$$\begin{array}{rl}{H}_f& =\sum _{j}2{\epsilon }_j{s}_j^z-g\sum _{j,k}{f}_j{f}_k{s}_j^{+}{s}_k^{-}\\ & =\sum _{j}2{\epsilon }_j{s}_j^z-|\Delta {|}^2/g.\end{array}$$

$$\begin{array}{rl}{H}_f& =\sum _{j}2{\epsilon }_j{s}_j^z-g\sum _{j,k}{f}_j{f}_k^{\ast }{s}_j^{+}{s}_k^{-}.\end{array}$$

$$\begin{array}{rl}{H}_f& =\sum _{j}2{\epsilon }_j{s}_j^z-g\sum _{j,k}|f_j||f_k|s_j^{+}{s}_k^{-},\end{array}$$

$$\begin{split}\dot{\mathbf{S}}_{\mathbf{k}\lambda}&=\mathbf{B}_{k\lambda}\times\mathbf{S}_{\mathbf{k}\lambda}+\mathbf{m}_{k}\times\mathbf{L}_{\mathbf{k}\lambda}-\mathbf{m}_{k}T_{\mathbf{k}},\\ \dot{L}^{x}_{\mathbf{k}\lambda}&=-2\varepsilon_{k}L^{y}_{\mathbf{k}\lambda}+\frac{m^{y}_{k}}{2}\big{[}S^{z}_{\mathbf{k}+}+S^{z}_{\mathbf{k}-}\big{]}+B^{x}_{k\lambda}T_{\mathbf{k}},\\ \dot{L}^{y}_{\mathbf{k}\lambda}&=2\varepsilon_{k}L^{x}_{\mathbf{k}\lambda}-\frac{m^{x}_{k}}{2}\big{[}S^{z}_{\mathbf{k}+}+S^{z}_{\mathbf{k}-}\big{]}+B^{y}_{k\lambda}T_{\mathbf{k}},\\ \dot{L}^{z}_{\mathbf{k}\lambda}&=-2R_{k}\lambda T_{\mathbf{k}}+\frac{m^{x}_{k}}{2}\big{[}S^{y}_{\mathbf{k}\lambda}-S^{y}_{\mathbf{k}\bar{\lambda}}\big{]}-\frac{m^{y}_{k}}{2}\big{[}S^{x}_{\mathbf{k}\lambda}-S^{x}_{\mathbf{k}\bar{\lambda}}\big{]},\\ \dot{T}_{\mathbf{k}}&=2R_{k}L^{z}_{\mathbf{k}+}-B^{x}_{k+}L^{x}_{\mathbf{k}+}-B^{y}_{k+}L^{y}_{\mathbf{k}+}+\\ &+\frac{1}{2}\mathbf{m}_{k}\cdot\big{[}\mathbf{S}_{\mathbf{k}+}+\mathbf{S}_{\mathbf{k}-}\big{]},\end{split}$$

$$\begin{split}\Delta&=\frac{g}{2}\sum_{\mathbf{k}\lambda}\big{[}\sin\phi_{k}S^{-}_{\mathbf{k}\lambda}+\cos\phi_{k}L^{-}_{\mathbf{k}\lambda}\big{]}\\ &=\Delta_{x}-i\Delta_{y},\\ \mathbf{B}_{k\lambda}&=(-2\sin\phi_{k}\Delta_{x},-2\sin\phi_{k}\Delta_{y},2\varepsilon_{k\lambda}),\\ \mathbf{m}_{k}&=(-2\cos\phi_{k}\Delta_{x},-2\cos\phi_{k}\Delta_{y},0).\end{split}$$

$$N^{2}_{\mathbf{k}}=2T^{2}_{\mathbf{k}}+\sum_{\lambda}\big{[}\mathbf{S}^{2}_{\mathbf{k}\lambda}+\mathbf{L}^{2}_{\mathbf{k}\lambda}\big{]}=\frac{1}{4}.$$

$$\begin{array}{rl}{H}_{so}& =\sum _{\mathbf{k}\lambda }2{\epsilon }_{k\lambda }{S}_{\mathbf{k}\lambda }^z-2|\Delta {|}^2/g.\end{array}$$

$$\begin{array}{r}{N}_f=\sum _j(2s_j^z+1),\end{array}$$

$$\begin{split}N_{f}=\sum_{\mathbf{k}\lambda}\bigg{(}S^{z}_{\mathbf{k}\lambda}+\frac{1}{2}\bigg{)},\end{split}$$

$$\begin{array}{rl}{s}_{j0}^{-}& =\frac{f_j{\Delta }_0}{2E_j},s_{j0}^z=-\frac{{\epsilon }_j-\mu }{2E_j},\\ E_j(\Delta )& \equiv \sqrt{({\epsilon }_j-\mu {)}^2+f^2({\epsilon }_j)|\Delta {|}^2}.\end{array}$$

$$\begin{split}\frac{1}{g}=\sum_{j}\frac{f_{j}^{2}}{2E_{j}},\quad N_{f}=\sum_{j}\bigg{(}1-\frac{\varepsilon_{j}-\mu}{E_{j}}\bigg{)}\end{split}$$

$$\begin{array}{r}{s}_j^z(t)\to s_j^z(-t),s_j^{\pm }(t)\to s_j^{\mp }(-t),\\ \Delta (t)\to {\Delta }^{\ast }(-t).\end{array}$$

$$\begin{split}S^{x}_{\mathbf{k}\lambda 0}&=\frac{\Delta_{0}\sin\phi_{k}}{D_{k}}\bigg{[}\Delta_{0}^{2}+\xi^{2}_{k\bar{\lambda}}+E_{k+}E_{k-}\bigg{]},\\ S^{z}_{\mathbf{k}\lambda 0}&=-\frac{1}{D_{k}}\bigg{[}\xi_{k\lambda}(E_{k+}E_{k-}+\xi^{2}_{k\bar{\lambda}}+\Delta_{0}^{2}\sin^{2}\phi_{k})+\\ &+\Delta_{0}^{2}\cos^{2}\phi_{k}\xi_{k\bar{\lambda}}\bigg{]},\\ L^{x}_{\mathbf{k}0}&=\frac{\Delta_{0}\cos\phi_{k}}{D_{k}}\bigg{[}\Delta_{0}^{2}+\xi_{k+}\xi_{k-}+E_{k+}E_{k-}\bigg{]},\\ L^{z}_{\mathbf{k}0}&=\frac{1}{D_{k}}\bigg{[}2R_{k}\Delta_{0}^{2}\cos\phi_{k}\sin\phi_{k}\bigg{]},\\ \xi_{k(\lambda)}&\equiv\varepsilon_{k(\lambda)}-\mu,\\ E_{k\lambda}(\Delta)&\equiv\bigg{[}\xi_{k}^{2}+\Delta^{2}+R_{k}^{2}-2R_{k}\lambda\sqrt{\xi_{k}^{2}+\cos^{2}\phi_{k}\Delta^{2}}\bigg{]}^{1/2},\\ D_{k}&\equiv 2E_{k+}E_{k-}(E_{k+}+E_{k-}),\end{split}$$

$$\begin{array}{r}\frac{2}{g}=\sum _{\mathbf{k}\lambda }\frac{E_{k+}{E}_{k-}+{\Delta }_0^2+{\sin }^2{\varphi }_k{\xi }_{\mathbf{k}\lambda }^2+{\cos }^2{\varphi }_k{\xi }_{\mathbf{k}\lambda }{\xi }_{\mathbf{k}\bar{\lambda }}}{2E_{k+}{E}_{k-}(E_{k+}+E_{k-})}.\end{array}$$

$$\delta =\frac{W}{N-1}.$$

$$\begin{array}{r}{f}_{\text{lor}}(\epsilon ,\gamma )=\frac{\gamma }{\sqrt{{\gamma }^2+{\epsilon }^2}},\end{array}$$

$$\begin{array}{r}{f}_{\text{sin}}(\epsilon ,\gamma )=1+{\sin }^2(\epsilon /\gamma ),\end{array}$$

$$\begin{array}{r}{f}_{\text{cub}}(\epsilon ,\gamma )=\frac{({\gamma }^3+|\epsilon {|}^3{)}^{1/3}}{\gamma }.\end{array}$$

$$t_{\mathrm{fs}}\approx \frac{1}{\delta }=\frac{N-1}{W},$$

$$\begin{array}{rl}{\dot{s}}_j^z& =0,\\ {\dot{s}}_j^x& =-2{\epsilon }_j{s}_j^y,\\ {\dot{s}}_j^y& =2{\epsilon }_j{s}_j^x.\end{array}$$

$$\begin{array}{rl}{s}_j^z& =z_j,\\ s_j^x& =x_j\cos (2{\epsilon }_{j}t),\\ s_j^y& =x_j\sin (2{\epsilon }_{j}t),\\ z_j^2& =1/4-x_j^2.\end{array}$$

$$\Delta =g_f\sum _j{f}_j{x}_j\cos (2{\epsilon }_{j}t)\ne 0,$$

$$\begin{array}{rl}{\dot{s}}_j^z& =-2f_j{s}_j^y{\Delta }_{\infty },\\ {\dot{s}}_j^x& =-2{\epsilon }_j{s}_j^y,\\ {\dot{s}}_j^y& =2{\epsilon }_j{s}_j^x+2f_j{s}_j^z{\Delta }_{\infty }.\end{array}$$

Full text

1. I Introduction

2. II Models and pseudospin representation

3. III Main results

4. IV Ground state and quench protocol

5. V Simulations of nonequilibrium phases and stability analysis

6. V.1 Phases I and II

7. V.2 Stability analysis

8. V.3 Phase III

9. V.3.1 Universality of elliptic oscillations

10. V.3.2 Relaxation time

11. VI Phase III asymptotic solution

12. VI.1 External driving

13. VI.2 Phase III spin solution in the separable BCS model

14. VI.3 Asymptotic self-consistency

15. VI.4 Self-consistent solutions in the separable BCS model

16. VII Quasiperiodic Phase IV

17. VIII Conclusion

18. A Mean-field equations of motion

19. B Integrable limit of spin-orbit quenches

20. C Integrability breaking forbids asymptotic reduction

21. C.1 Existence of reduced solutions implies integrability and vice versa

22. C.2 Asymptotic $\Delta(t)$ does not match the 2-spin solution in nonintegrable cases

23. D The link between Lax constructions and the stability analysis

24. D.1 Lax norms

25. D.2 Phase I-II transition

26. D.3 Phase II-III transition

27. D.4 Real parts of Lax roots at the transitions

28. D.4.1 $s$-wave, II-III

29. D.4.2 $p+ip$, II-III

I Introduction

The past fifteen years have borne witness to impressive advances in the ability to experimentally control many-body systems where dissipative and decoherence effects are strongly suppressed. Studies of cold atomic gaseskinoshita2006 ; ligner ; Hofferberth:2007 ; weiler ; widera ; Gring:2012 ; Langen:2015 ; tang ; norcia ; smale , solid state pump-probe experimentsftdk ; ktn ; gcff ; mhmutws ; mtfsmutwas and quantum information processingriberio ; gorshkov0 ; zhou ; barends ; nichol ; song ; lukin ; gorshkov can now explore coherent many-body dynamics for long time scales, paving the way for the characterization of new phenomena. In particular, cold atomic gases with tunable interactionsjin ; ketterle ; bdz ; zass ; gps ; bdn are an instrumental experimental tool in the quest to understand previously inaccessible aspects of far from equilibrium many-body dynamics.

A major focus of recent theory and experiment has been the unitary time evolution of a system, initially in the ground state, subject to a sudden perturbationpssv ; gogolin ; moore . This experimental protocol, known as a quantum quench, can induce long-lived states with properties strikingly different from those of equilibrium states at similar energy scales. In this work, we focus on the quench dynamics of various superconducting models, which is a modern reformulation of the longstanding problem of nonequilibrium superconductivity in the collisionless regimeanderson ; galaiko ; volkov ; galperin . A canonical result is that the infinitesimal perturbation of a Bardeen-Cooper-Schrieffer (BCS) $s$-wave superconductor leads to power law oscillatory relaxation of the order parameter amplitude $|\Delta|$ to a constant valuevolkov .

Decades later, it was discovered that larger deviations could give rise to different dynamical phases identified by the asymptotic behavior of the amplitude of the order parameterbarlev0 ; amin ; yuzbashyan-jpa ; Szymanska ; yuzbashyan2005 ; yta ; barlev ; yuzdze . Consider the dynamics of $\Delta$ after quenches of the coupling $g$ in various superconducting models. When the final coupling $g_{f}$ is small enough, $\Delta$ vanishes rapidly in time; this behavior characterizes what we call Phase I. For intermediate $g_{f}$, $|\Delta|$ exhibits oscillatory power law decay to a nonzero constant (Phase II). For larger $g_{f}$, $|\Delta|$ exhibits persistent periodic oscillations (Phase III) – a nonlinear manifestation of what is known in the literature as the Higgs or amplitude modebarlev3 ; pekvar ; pashleit ; barvar ; nmfkd ; hkcvak ; kbums ; mve .

The exact quantum quench phase diagrams of the $s$-wave superconductor were eventually constructed using a sophisticated analytical method that relies on the model’s integrabilityydgf . It turns out that the integrable $p+ip$ topological superconductor exhibits the same three phases, and similar analytical tools lead to the construction of its phase diagramsfdgy . Thus, there may appear to be some profound connection between integrability and these three dynamical phases, but nonintegrable models also have Phases I and IIbarlev3 ; Szymanska ; ddgp ; dky ; psc and Phase III-like behavior is thought to persist in some nonintegrable models as well. On the other hand, the existence of Phase III in such models has not been convincingly established beyond the linear regime and aspects of quench dynamics unique to the nonintegrable case have not been explored.

Overall, the description of these nonequilibrium dynamical phases lacks a unifying mechanism applicable to finite quenches of nonintegrable pairing models. Here we present an in-depth study of the nonequilibrium phases of various nonintegrable superconducting models with and without spin-orbit coupling. A common feature of models we consider is that the order parameter takes the form of a single complex number. We establish that Phase III persists when integrability is broken note and give strong numerical evidence that the persistent oscillations are always elliptic, which generalizes the known behavior of integrable modelsbarlev0 ; ydgf ; fdgy .

Although the integrable and nonintegrable phenomenology are similar, we find that integrability breaking has profound consequences. Unique to nonintegrable models is an extremely long relaxation time scale $\tau$ which diverges as one approaches integrable points and is most prominent in quenches to Phase III. One must analyze dynamics beyond $\tau$ to truly observe Phase III, which has not been done in other studies. As illustrated in Fig. 1, for $t<\tau$, $|\Delta|$ may oscillate with several frequencies and a slowly evolving amplitude, both of which undermine naive analyses restricted to $t<\tau$. One may incorrectly conclude from the transient dynamics that the asymptotic nonequilibrium phase has several undamped frequencies, or that $|\Delta|$ is oscillating periodically while in fact the amplitude is still changing. Nonintegrable Phase III oscillations further require comparatively more elaborate elliptic functions to describe the oscillations.

To complicate the picture even further, certain quantum quenches of nonintegrable pairing models genuinely do not fit into any of the Phases I, II and III. Here the asymptotic $|\Delta|$ is truly quasiperiodic, leading us to conclude that there are regions of quasiperiodicity – a new Phase IV – in the quantum quench phase diagrams of these models.

Another consequence of integrability breaking arises in the analytical description of the three nonequilibrium phases. In the integrable case, there is a dynamical reduction in the number of degrees of freedom of the systemydgf ; fdgy such that Phases I, II and III correspond to an effective classical spin Hamiltonian with 0, 1 and 2 spins, respectively. Phase III in the general case, however, does not admit such a 2-spin representation. As a surrogate to this analytical method, we propose a stability analysis of Phases I and II that applies generally to finite quenches. The stability analysis is based on linearizing around the asymptotic solutions to the equations of motion in each of the phases. We can then nonperturbatively determine the phase I-II boundary as well as the phase II-III boundary in nonintegrable pairing models. Finally, we return to Phase III and argue that the self-consistency condition (gap equation) is responsible not only for the existence of persistent periodic oscillations of $|\Delta|$, but also for selecting elliptic functions amongst all possible periodic functions.

II Models and pseudospin representation

In this paper, we consider quantum quenches in two types of nonintegrable pairing models

[TABLE]

The Hamiltonian $\hat{H}_{f}$ is a separable BCS Hamiltonian where the $\varepsilon_{j}$ are $N$ single-particle energy levels, $\hat{c}^{{\dagger}}_{j\lambda}$ ($\hat{c}_{j\lambda}$) is a fermion creation (annihilation) operator for an electron with energy $\varepsilon_{j}$ and spin index $\lambda$, $g>0$ is the pairing interaction strength and $f_{j}\equiv f(\varepsilon_{j})$ is a generic function of $\varepsilon_{j}$. The Hamiltonian $\hat{H}_{so}$ describes a 2D topological spin-orbit coupled superconductor with $s$-wave interactionsso1 ; so2 . Here $\mathbf{k}=(k_{x},k_{y})$ is a two-dimensional momentum vector, $\sigma^{j}$ are Pauli matrices, $h$ is a Zeeman field and $\alpha$ is the Rashba spin-orbit coupling. We will take the density of states to be constant for both models, which is the case in 2D or at weak coupling, so that the single-particle energy levels are distributed uniformly on an interval of length $W$, called the bandwidth.

Apart from certain choices of $f(x)$, the separable BCS Hamiltonian $\hat{H}_{f}$ is a toy model for breaking integrability. The choice of $f^{2}(x)=C_{1}+C_{2}x$ produces a quantum integrable Hamiltonianrichardson ; ortiz2005 ; for example, $f(x)=1$ and $f(x)=\sqrt{x}$ correspond to the $s$-waveyuzbashyan-jpa and $p+ip$dunning ; dukelsky BCS models, respectively. A notable nonintegrable case is the $d+id$ modelMarquette , where $f(x)=x$. The spin-orbit Hamiltonian $\hat{H}_{so}$, on the other hand, can be realized with cold Fermi gasesmamatu ; ffvkzk ; dkwhmhlhv ; wyfmhczz ; cshybz ; wbljs ; fhmwzzzzz ; mamatu2 ; hmwpzclzz ; huizhai .

As both Hamiltonians in Eq. (1) have infinite range interactions, the mean-field approximation is expected to be exact in the thermodynamic ($N\to\infty$) limit. We therefore replace 2-body operators as follows $\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\hat{c}\hat{c}\approx\langle\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\rangle\hat{c}\hat{c}+\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\langle\hat{c}\hat{c}\rangle-\langle\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\rangle\langle\hat{c}\hat{c}\rangle$ in the equations of motion. We also diagonalize the noninteracting part of $\hat{H}_{so}$ through a unitary transformation $U_{\mathbf{k}}$ which is detailed in Appendix A. Up to additive constants, the effective mean-field Hamiltonians of Eq. (1) are

[TABLE]

The new parameters in $\hat{H}_{so}$ are

[TABLE]

Note that both $\alpha=0$ and $h=0$ correspond to integrable points of the spin-orbit model; in both cases, $\hat{H}_{so}$ becomes a Hamiltonian for two bands of independent $s$-wave BCS models. Most importantly, the mean-field order parameters $\Delta\equiv\Delta(t)$ are defined in terms of expectation values

[TABLE]

for their respective models.

We will discuss the mean-field dynamics generated by Hamiltonians (2) in terms of Anderson pseudospins $\mathbf{\hat{s}}_{j}=(\hat{s}^{x}_{j},\hat{s}^{y}_{j},\hat{s}^{z}_{j})$ which will allow for intuitive visualizations of the dynamics of different nonequilibrium phases. The transformation from fermions to pseudospins is given by

[TABLE]

In the spin-orbit case the pseudospin representation requires an additional set of auxiliary variables. For the sake of brevity, we relegate the derivations of the pseudospin equations of motion to Appendix A and simply state them here.

In the mean-field equations of motion that follow, $\mathbf{s}=\langle\mathbf{\hat{s}}\rangle$ are to be understood as classical variables satisfying the angular momentum Poisson brackets $\{s^{a}_{j},s^{b}_{k}\}=-\delta_{jk}\epsilon_{abc}s_{j}^{c}$. In the separable BCS model, we have

[TABLE]

where self-consistency requires

[TABLE]

The spin-length $s_{j}=1/2$ is conserved by Eqs. (6), which together with Eq. (7) are the equations of motion of the following classical spin Hamiltonian:

[TABLE]

Note that without loss of generality, we can choose $f_{j}$ to be real and nonnegative as we have done above. Indeed, let $f_{j}=|f_{j}|e^{-i\theta_{j}}$ be general complex numbers and

[TABLE]

We redefine the spins by making local rotations around the z-axis, $s^{-}_{j}\to s^{-}_{j}e^{-i\theta_{j}}$. In terms of the new spins the Hamiltonian becomes

[TABLE]

and the order parameter is $\Delta=\sum_{j}|f_{j}|s^{-}_{j}$. This transformation does not affect spin (angular momentum) Poisson brackets and therefore the equations of motion retain their form. We thus arrive at the same problem only with $f_{j}\to|f_{j}|$.

We use capital letters $\mathbf{S}_{\mathbf{k}\lambda}$ to denote the classical pseudospins in the spin-orbit model and must introduce (see Appendix A) a set of auxiliary variables: the scalars $T_{\mathbf{k}}$ and vectors $\mathbf{L}_{\mathbf{k}\pm}$, where $\mathbf{L}_{\mathbf{k}+}$ and $\mathbf{L}_{\mathbf{k}-}$ differ only in sign of the z-component. The equations of motion are

[TABLE]

where the momentum dependent fields $\mathbf{B}_{k\lambda}$ and $\mathbf{m}_{k}$ are defined in terms of the order parameter $\Delta$

[TABLE]

The first of these equations is the self-consistency relationship for the spin-orbit model. The equation for $\dot{\mathbf{S}}_{\mathbf{k}\lambda}$ in Eq. (11) corrects an error in a previous paperdky , which is missing the last term. For each $\mathbf{k}$, there is a conserved quantity analogous to pseudospin length

[TABLE]

Similar to Eq. (8), the classical spin-orbit Hamiltonian in pseudospin notation has a simple and compact expression

[TABLE]

Because of the simple relationship connecting $\mathbf{L}_{\mathbf{k}+}$ to $\mathbf{L}_{\mathbf{k}-}$, each momentum vector $\mathbf{k}$ corresponds to ten dynamical variables $(\mathbf{S}_{\mathbf{k}+},\mathbf{S}_{\mathbf{k}-},\mathbf{L}_{\mathbf{k}+},T_{\mathbf{k}})$ constrained by Eq. (13). Note that $T_{\mathbf{k}}$ and $L^{z}_{\mathbf{k}\lambda}$ do not appear in (14), but as discussed in Appendix A, they are necessary for the closure of the equations of motion. From now on we simplify notation to $\mathbf{L}_{\mathbf{k}}\equiv\mathbf{L}_{\mathbf{k}+}$ and define the 10-dimensional vector $\mathbf{\Gamma}_{\mathbf{k}}\equiv(\mathbf{S}_{\mathbf{k}+},\mathbf{S}_{\mathbf{k}-},\mathbf{L}_{\mathbf{k}},T_{\mathbf{k}})$.

Finally, the conservation of the total number of fermions $N_{f}$ in each model corresponds to the conservation of total z-component in the pseudospin language

[TABLE]

for the separable BCS model and

[TABLE]

for the spin-orbit model.

III Main results

The main purpose of this work is to compare the nonequilibrium phases of quenches from the ground state of nonintegrable pairing Hamiltonians, such as those in Eq. (1), to those of the integrable $s$-waveydgf and $p$-wavefdgy models. Some qualitative aspects of the primary phases are independent of integrability insofar as the squared modulus of the order parameter $\Delta$ may exhibit any of three distinct asymptotic behaviors in the continuum limit: it can relax to zero (Phase I), relax to a nonzero constant value (Phase II), or display persistent periodic elliptic oscillations (Phase III).

We first show the existence of these three phases in Sects. V.1-V.3 through direct numerical simulation of the dynamics. In Sect. V.2 we present a stability analysis of the phases of the separable BCS models which leads to conditions for nonequilibrium phase transitions. The stability analysis applied to integrable cases reduces to the known results that relied on exact solvabilityfdgy ; ydgf . Our analysis provides a physical explanation for the transitions in terms of the frequencies of linearized perturbations $\delta\Delta(t)$ of the asymptotic $\Delta$. The transition from Phase I to Phase II occurs through an exponential instability characterized by a pair of conjugate imaginary frequencies in the linearization spectrum, while that of Phase II to III occurs either when small harmonic oscillations fail to dephase or when an exponential instability occurs.

The appearance of some or all of Phases I-III in nonintegrable models suggests an underlying universality to quench dynamics, but we show that the story is less straightforward. One the one hand, these phases are understood in the integrable casesydgf ; fdgy . There is a dynamical reduction of the number of effective degrees of freedom, so that at large times the dynamics are governed by a Hamiltonian of the same form, but which has just a few collective degrees of freedom. The three phases correspond to 0, 1 or 2 effective spins for each phase, respectively. On the other hand, the nonintegrable dynamics admit no known analogous reduction because the 2-spin solutions to the equations of motion do not reproduce the observed asymptotic behavior of $\Delta$ in Phase III. If such a reducing “flow” in time of the Hamiltonian occurs in the nonintegrable case, then the form of the Hamiltonian itself must change. For specifics on this latter point, see Appendix C.

Importantly, nonintegrable pairing models also display dynamics markedly different from those in the main three phases. We illustrate this behavior with two examples in Sect. VII – one for the spin-orbit Hamiltonian and one for a particle-hole symmetric separable BCS Hamiltonian – where the magnitude of the order parameter oscillates quasiperiodically. We interpret this observation as an indication of a new quasiperiodic phase (Phase IV) unique to quantum quench phase diagrams of these models.

More subtle details of the dynamics in the main three phases change drastically once integrability is broken. We show in Sect. V.3.2 that nonintegrable models take an extremely long time to relax to Phase III. This time scale is absent in the integrable case, yet it diverges when one approaches the integrable limit. One must take this time scale into account when studying Phase III on the basis of numerical simulation alone. For example, in the weak coupling regime, the nonintegrable $d+id$ model may appear to quickly enter Phase IIIdidShortTime while in fact the minima of $|\Delta|$ oscillations have not converged to a fixed value. The further into the weak coupling regime one explores, the longer the relaxation time. Quenches outside of weak coupling have faster dynamics, but exhibit behavior that markedly contrasts with Phase III, and above a certain energy threshold the asymptotic state collapses rapidly to Phase II. This long relaxation time is typical in the nonintegrable case.

Despite these consequences of breaking integrability, our mixed strategy of simulation and stability analysis applies to the two rather different classes of nonintegrable pairing models found in Eq. (1). The separable BCS permits a standard Anderson pseudospin representation and is a single band model, while the spin-orbit model requires an expanded pseudospin representation, has multiple bands and a topological quantum phase transition. Yet both models have a single complex order parameter, which we believe is the essential characteristic that leads to the three phases.

The self-consistency relationship (7) for the order parameter is central to both our stability analysis of Phases I and II in Sect. V.2 and our investigations of Phase III in Sect. VI. In the former case, the frequencies of harmonic perturbations of a given nonequilibrium phase are constrained by the self-consistency requirement. As for Phase III, we show in Sect. VI that there is always a periodic solution to the spin equations of motion when $\Delta(t)$ is periodic, and that the general spin solution precesses around the periodic one. We then argue through numerical examples that further imposing the self-consistency requirement on $\Delta(t)$ selects elliptic functions amongst all possible periodic $\Delta(t)$.

IV Ground state and quench protocol

In a quantum quench, we prepare the system in the ground state with an initial order parameter $\Delta=\Delta_{0}e^{-2i\mu t}$, which corresponds to system parameters such as the interaction strength $g$, the equilibrium chemical potential $\mu$, the magnetic field $h$ and the spin-orbit strength $\alpha$. The amplitude $\Delta_{0}$ is constant in the ground state. At time $t=0$, we suddenly change one of these parameters, which throws the system out of equilibrium. In the separable BCS model we will consider quenches $g_{i}\to g_{f}$, but we will label the initial and final states by the coordinates $\Delta_{0i}\equiv\Delta_{0}(g_{i})$ and $\Delta_{0f}\equiv\Delta_{0}(g_{f})$. In the spin-orbit model, we will consider quenches of the magnetic field $h_{i}\to h_{f}$. The fermion number $N_{f}$ is fixed across the quench in both cases, which implies that the equilibrium chemical potential $\mu$ changes with $h$.

For a given $\Delta_{0}$ and $\mu$, we express the ground state configuration of the separable BCS model in a frame that rotates around the z-axis with frequency $2\mu$. We then orient each $\mathbf{s}_{j}$ against the magnetic $\mathbf{b}_{j}$, the z-component of which is shifted by $2\mu$,

[TABLE]

The relationship between $\Delta_{0}$, $g$, $N_{f}$ and $\mu$ obtains from the application of the definition of $\Delta$ in (7) to (15) and the configuration in (17),

[TABLE]

We will assert without loss of generality that $\Delta_{0i}$ is real in both models, which can always be achieved by a time-independent rotation in the $xy$-plane in pseudospin space.

Unless otherwise stated, we will simplify the analysis of the separable BCS model by restricting ourselves to cases where the order parameter $\Delta$ remains real for all time, i.e., $\Delta_{y}(t)=0$. To achieve this, we will consider the particle-hole symmetric case where the energies $\varepsilon_{j}$ are symmetrically distributed around the chemical potential $\mu$, which is set to zero without loss of generality. We will also only consider even functions $f(x)=f(-x)$. Under these conditions, any initial spin configuration that satisfies the symmetry conditions $s^{z}(\varepsilon_{j})=-s^{z}(-\varepsilon_{j})$, $s^{+}(\varepsilon_{j})=s^{-}(-\varepsilon_{j})$, as does the ground state (17), will do so for all time. This fact can be verified with the equations of motion (6) by considering time derivatives of quantities such as $s^{z}(\varepsilon_{j})+s^{z}(-\varepsilon_{j})$, which vanish under the aforementioned assumptions. We will not use particle-hole symmetry in the $d+id$ model, where $f(x)=x$ and $\varepsilon_{j}$ will be distributed on a positive interval. Further, Eqs. (6) and (7) are invariant under the time-reversal transformation

[TABLE]

Since the initial conditions (17) at $t=0$ also have this property, it holds at all times.

The ground state of the spin-orbit model is less obviousdky

[TABLE]

while $S^{y}_{\mathbf{k}\lambda 0}=L^{y}_{\mathbf{k}0}=T_{\mathbf{k}0}=0$. The corresponding self-consistent equation relating $\Delta_{0}$ to $g$ is

[TABLE]

The quantities $2E_{j}(\Delta)$ and $2E_{\mathbf{k}\lambda}(\Delta)$ in (17) and (20) are the excitation energies obtained by diagonalization of the quadratic mean-field Hamiltonians in Eqs. (2) at a given $\Delta$.

For given values of $g$, $N_{f}$, $\alpha$ and $h$, one can simultaneously solve Eq. (16) and Eq. (21) using the ground state configurations to obtain the corresponding chemical potential $\mu$ and ground state gap $\Delta_{0}$. As the ground state is rotationally symmetric in $\mathbf{k}$, and the equations of motion preserve this symmetry, in our numerics we always replace sums over momenta with sums over energies with a flat density of states $\sum_{\mathbf{k}}\to\sum_{\varepsilon}$. The level spacing $\delta$ is related to the number of spins $N$ and the bandwidth $W$ through

[TABLE]

Formally, in 2D this means $N-1=\frac{W}{2\pi}A$, where $A$ is the physical area of the system. Fig. 2 shows an example of the relationship between different parameters for the spin-orbit model.

V Simulations of nonequilibrium phases and stability analysis

Now we numerically simulate the equations of motion (6) and (11) and plot the behavior of $\Delta(t)$ for each of the three phases in Sects. V.1 and V.3. In Sect. V.3, we also characterize the long time scale of nonintegrable models in Phase III. In Sect. V.2, we introduce a stability analysis for the Phases I and II that gives the conditions under which a nonequilibrium phase transition occurs.

We will consider several integrability-breaking functions for $f(\varepsilon)$, which appears in the separable BCS equations of motion Eq. (6). All $f(\varepsilon)$ considered here will be even functions, and as we discuss in Sect. V.2, the particular form of $f(\varepsilon)$ affects which phases occur. With this in mind, we consider the “Lorentzian” couplingbarlev3

[TABLE]

the “sine” coupling,

[TABLE]

and the “cube root” coupling,

[TABLE]

The parameter $\gamma$ is fixed for any particular Hamiltonian, and it characterizes how strongly integrability is broken. For $\gamma\gtrsim W$, we have $f(\varepsilon,\gamma)\sim 1$ in all three cases, which we consider to be “nearly integrable”. For $\gamma\ll W$, integrability is strongly broken.

We control for finite size effects in our simulations by increasing $N$ until $\Delta(t)$ in the time window of interest no longer changes when $N$ is doubled. In practice, we find that finite size effects become significant at times $t>t_{\mathrm{fs}}$, where

[TABLE]

is the inverse single-particle level spacing, see also Ref. ydgf, . To observe the asymptotic dynamics, $N$ has to be sufficiently large, so that the relaxation time $\tau<t_{\mathrm{fs}}$.

V.1 Phases I and II

Figs. 3-5 contain examples of Phase I and Phase II quenches in both the separable BCS and spin-orbit models.

To heuristically understand the emergence of these two phases, one can insert the prescribed behavior of $\Delta$ into the equations of motion (6) and (11). This examination of the asymptotic solutions to the equations of motion in each phase will be important for the stability analyses of Sect. V.2.

The following applies to the separable BCS models in the particle-hole symmetric limit, but the analysis is analogous when this symmetry is broken and in the spin-orbit case. In Phase I, we set $\Delta$ to zero

[TABLE]

The most general solution that conserves both $s_{j}^{2}=1/4$ and the time-reversal symmetry (19) is

[TABLE]

where $z_{j}$ is the Phase I steady state spin-profile. In order for (28) to make sense as a solution to the actual equations of motion, Eq. (7) must hold, i.e., we must have that $\Delta=g_{f}\sum_{j}f_{j}s_{j}^{-}$ equals zero, which is called the self-consistency condition. Strictly speaking, the solution (28) violates the self-consistency condition

[TABLE]

but as the number of single-particle energies $N$ goes to infinity, i.e., in the continuum limit when the sum in Eq. (29) turns into an integral, $\Delta$ from Eq. (29) vanishes through dephasing for $1\ll t\ll 1/\delta=(N-1)/W$. This description is invalid for $t\sim N/W$. In this sense, we refer to the solution (28) as asymptotically self-consistent, which is a concept we will often use in the remainder of this paper.

Let us now replace $\Delta$ with $\Delta_{\infty}\neq 0$ in Eq. (6) to examine the asymptotic solutions corresponding to Phase II

[TABLE]

The solution which preserves spin length and the time-reversal symmetry is then

[TABLE]

where $Z_{j}$ is the Phase II steady state spin profile, which, along with $\Delta_{\infty}$, determines the other constants

[TABLE]

The solution (LABEL:eomHFphase2sol) must be asymptotically self-consistent, i.e., for $N\to\infty$, $\lim_{t\to\infty}\Delta=\Delta_{\infty}$, which implies

[TABLE]

which is the nonequilibrium analogue of the ground state self-consistency requirement (18).

V.2 Stability analysis

Now we consider the stability of Phases I and II for the separable BCS model by linearizing the equations of motion (27) and (30) about the asymptotic states given in (28) and (LABEL:eomHFphase2sol), respectively. The main result is Eq. (42), which is the equation for frequencies of linearized perturbations to the asymptotic $\Delta(t)$ of either Phase I or Phase II. For Phase I, the appearance of a complex conjugate pair of imaginary frequencies signals an exponential instability. For Phase II, a solution $\omega_{0}$ to Eq. (42) may enter the band gap, or a complex conjugate pair of frequencies may appear. The former case, which occurs in the integrable $s$-wave and $p+ip$ models, signifies a transition to Phase III because the linearized gap $\delta\Delta(t)$ oscillates persistently, i.e., it does not dephase. In Appendix D, we show that the nonequilibrium phase transitions predicted by this stability analysis both match and give a physical interpretation to the results obtained in integrable modelsydgf ; fdgy using tools inextricably linked to exact solvability.

Although the final result (42) applies generally, we limit the discussion to the particle-hole symmetric case to simplify the presentation. Let $\mathbf{s}_{j}=\mathbf{s}_{j0}+\delta\mathbf{s}_{j}$, where $\mathbf{s}_{j0}$ is the Phase I asymptotic solution from Eq. (28). Neglecting second and higher order terms, the linearized equations for the spin components are

[TABLE]

Expanding $\mathbf{s}_{j}(t)$ in Fourier components

[TABLE]

and using the Fourier space version of the self-consistency relation in Eq. (34), we find the following equation for the allowable frequencies $\omega$

[TABLE]

The following discussion uses particle-hole symmetry along with the empirical fact that for quenches from the ground state, $z_{j}\varepsilon_{j}<0$ in Phase I. Upon inspecting Eq. (36), one determines that there are $N/2$ unique $\omega_{j}^{2}$, of which all but one lie between consecutive $4\varepsilon_{j}^{2}$. The remaining $\omega_{0}^{2}$ is less than the smallest $4\varepsilon_{j}^{2}$, and can therefore be negative. A negative $\omega_{0}^{2}$ corresponds to a pair of conjugate imaginary frequencies, and therefore an exponential instability in $\delta\mathbf{s}_{j}$. We thus determine the Phase I boundary in $(\Delta_{0i},\Delta_{0f})$ space to be those values for which $\omega_{0}^{2}$ passes through zero.

The stability analysis for Phase II follows a similar logic. Consider the linearized equations of motion

[TABLE]

where now $\mathbf{s}_{j0}$ is the Phase II asymptotic solution from Eq. (LABEL:eomHFphase2sol). Again changing to the Fourier basis, we solve for $\delta\widetilde{s}_{j}^{x}(\omega)$ and apply the self-consistency condition for $\delta\widetilde{\Delta}(\omega)$, which reads

[TABLE]

Although in principle Eq. (38) can be solved numerically with $Z_{j}$ and $\Delta_{\infty}$ as input, such an approach is needlessly complex and obscures the mechanism by which Phase II gives way to Phase III. The difficulty presented by Eq. (38) stems from the fact that we required exact self-consistency. It turns out that relaxing this requirement to asymptotic self-consistency, defined in Sect. V.1, suffices to understand the Phase II-III transition.

We return to Eq. (37) and solve it in the time domain under the assumption $\delta{\Delta}(t)=\delta_{+}e^{-i\omega_{0}t}+\delta_{-}e^{i\omega_{0}t}$. We neglect higher order harmonics because the Phase III oscillations near the II-III boundary are small. Under this ansatz, $\delta s_{j}^{x}(t)$ has six frequencies: $\pm\omega_{0}$ and $\pm\omega_{0}\pm b_{j}$. If $\omega_{0}$ is a real frequency isolated from the continuum of $b_{j}$ defined in Eq. (32), then the constant $\Delta_{\infty}$ of Phase II is “unstable” in the sense that oscillatory perturbations do not dephase. The self-consistent equation for this harmonic $\delta\Delta(t)$ is

[TABLE]

This relation cannot hold for arbitrary $t$, but it will in the continuum limit if we require $\omega_{0}^{2}<b^{2}_{\textrm{min}}$ and $t\to\infty$, which allows the harmonic ansatz to be asymptotically self-consistent due to dephasing. Thus the equation for $\omega_{0}$, the frequency of a harmonic perturbation to $\Delta_{\infty}$ in Phase II, is

[TABLE]

Eq. (40) generalizes the small quench linearization method developed in Ref. barlev3, , which we recover by replacing $Z_{j}$ of Eq. (40) with the z-component spin profile of the $g_{i}$ ground state. For the Lorentzian coupling, $\omega_{0}$ is in the band gap for infinitesimal quenches, so that linearized Phase III oscillations do not decaybarlev3 .

In order to understand whether the finite quench dynamics admit such an isolated $\omega_{0}$, consider the implications of (40) combined with (33) for the $\Delta_{\infty}$ of Phase II. We find

[TABLE]

It helps to analyze (41) under the simplifying assumption that $Z_{j}/\varepsilon_{j}<0$, which holds exactly for the integrable $s$-wave model, and is therefore applicable in the weak-coupling regime ($\Delta_{0i},\Delta_{0f}\ll W$) of the general separable caseydgf . With this restriction, Eq. (40) implies $\omega_{0}^{2}$ is real, while Eq. (41) requires $\omega_{0}^{2}>0$, i.e., the allowed frequencies $\omega_{0}$ are purely real. We now examine the effect of the function $f_{j}$ in determining whether solutions $\omega_{0}^{2}$ to Eq. (41) are isolated from the $b_{j}^{2}$ continuum.

If $f_{j}<f(0)$ for all $j$ and $b^{2}_{\textrm{min}}=4\Delta_{\infty}^{2}$, then Eq. (41) has a solution $0<\omega_{0}^{2}<b^{2}_{\textrm{min}}$, and oscillations of $\delta\Delta(t)$ do not dephase. In this scenario, Phase III is the asymptotic state due the presence of persistent periodic oscillations about the Phase II solution. If $f_{j}<f(0)$ for all $j$ and $b^{2}_{\textrm{min}}<4\Delta_{\infty}^{2}$, then the relationship between $\omega_{0}^{2}$ and $b^{2}_{\textrm{min}}$ is not immediately obvious from Eq. (41). The Lorentzian coupling, where $f_{j}=\gamma(\gamma^{2}+\varepsilon_{j}^{2})^{-1/2}$, allows for both possibilities: If $\Delta_{\infty}\leq\gamma$, then $b^{2}_{\textrm{min}}=4\Delta_{\infty}^{2}$ and Phase II is not the asymptotic state. If $\Delta_{\infty}>\gamma$, then $b_{\textrm{min}}^{2}=4\gamma(2\Delta_{\infty}-\gamma)$, and we cannot characterize solutions to Eq. (41) without detailed knowledge of $Z_{j}$ and $\Delta_{\infty}$.

If $f_{j}\geq f(0)$ for all $j$, then $b^{2}_{\textrm{min}}=4\Delta_{\infty}^{2}$ and we find that solutions $\omega_{0}^{2}$ to Eq. (41) are not isolated from the $b_{j}^{2}$ continuum. In this case, the harmonic ansatz for $\delta\Delta(t)$ is not asymptotically self-consistent, and there are no persistent small oscillations about Phase II. The integrable $s$-wave model is defined by $f_{j}=f(0)=1$, in which case $\omega_{0}^{2}=4\Delta_{\infty}^{2}$ is the only solution to Eq. (41), which is not isolated. On the other hand, Phase III exists in the $s$-wave caseydgf . Therefore, $f_{j}\geq f(0)$ does not imply that such models will always reach Phase II. Indeed, the relaxation to Phase II is always accompanied by nonperturbative oscillations which persist in the case of Phase III.

Thus, even under the simplifying assumptions of particle-hole symmetry and $Z_{j}/\varepsilon_{j}<0$, the stability analysis of Phase II reveals a variety of possible behaviors in the separable BCS models. The nature of $f(\varepsilon)$ near $\varepsilon=0$ (the Fermi surface) is especially crucial to determining whether oscillations fully dephase to Phase II – a statement which extends to the non-particle-hole symmetric case in the weak coupling regime.

Upon relaxing the restriction $Z_{j}/\varepsilon_{j}<0$, isolated solutions to Eq. (41) can have nonzero imaginary part, thereby allowing for the possibility of exponential instabilities to Phase II solutions (see Fig. 10). In the non-particle-hole symmetric case, $\Delta(t)=\Delta_{\infty}e^{-2i\mu_{\infty}t}$ in Phase II, and the equation for the frequencies of harmonic $\delta\Delta(t)$ can be expressed in the form

[TABLE]

where

[TABLE]

The self-consistency equation for $\Delta(t)$ in Phase II has the same form as Eq. (33), with the substitution $\varepsilon_{j}\to\widetilde{\varepsilon}_{j}$. In the particle-hole symmetric limit, $S_{2}(\omega_{0})=0$ and the correct solution to Eq. (42) solves Eq. (40). In the limit $\Delta_{\infty}\to 0$, (42) is also the stability equation for Phase I. In Appendix D, we show that the Phase I-II and Phase II-III transitions given by (42) are identical to those obtained using exact solvability in the integrable $s$-wave and $p+ip$ models.

V.3 Phase III

V.3.1 Universality of elliptic oscillations

The asymptotic Phase III solution is significantly more complicated than its Phase I and Phase II counterparts (28) and (LABEL:eomHFphase2sol). We derive this solution in Sect. VI. Presently we provide evidence that the asymptotic behavior of $\Delta(t)$ can always be described by Jacobi elliptic functions. Consider first the particle-hole symmetric limit, for which we find

[TABLE]

where $P_{4}[\Delta(t)]$ is a generic fourth-order polynomial in $\Delta(t)$. Now parametrize $P_{4}[\Delta(t)]$ as

[TABLE]

where the real coefficients $\Delta_{\pm}$ are the maximum and minimum values of $\Delta(t)$, while $\widetilde{\Delta}_{\pm}$ are either complex conjugate or independent real numbers. This parametrization leads to the following solution for $\Delta(t)$

[TABLE]

where $\textrm{dn}[t,m]$ is the Jacobi-dn function. When particle-hole symmetry does not hold, then one replaces $\Delta(t)$ with $|\Delta(t)|^{2}$ in Eqs. (44)-(46). In Figs. 6 and 7 we show that Phase III oscillations in separable BCS models satisfy Eq. (44) and Eq. (46), while Fig. 8 shows the same for the spin-orbit model.

As a general rule of thumb, most spin-orbit quenches that superficially appear to relax to Phase III really have not. Fig. 8 is the result of a thorough search of the parameter space in order to find a true Phase III quench within a computationally achievable time. On the one hand, the final field $h_{f}$ has to be large enough so as to nonperturbatively break integrability, for small perturbations lead to long relaxation times. On the other hand, the fields cannot be so large as to suppress the equilibrium gap $\Delta_{0}$ scale, which is the scale of the oscillation frequency. The value of $\alpha$ must also break integrability nonperturbatively, but a larger $\alpha$ also requires a larger number of spins to reach the thermodynamic limit. Finally, it turns out that a smaller Fermi energy relative to the bandwidth promotes a faster relaxation time. We discuss this Phase III relaxation time further in Sect. V.3.2 in the context of the separable BCS models.

For the integrable $s$-wave case it can be shown analyticallyydgf that $\widetilde{\Delta}_{\pm}=\Delta_{\pm}$ and $a=1$, which greatly simplifies $P_{4}[\Delta(t)]$ and $\Delta(t)\to\Delta_{+}\textrm{dn}[\Delta_{+}(t-t_{0}),1-\frac{\Delta^{2}_{-}}{\Delta^{2}_{+}}]$. The mechanism behind the emergence of the three phases in the $s$-wave Hamiltonian is a dynamical reduction in the number of degrees of freedom. The Phase III asymptotic solution for $\Delta(t)$ is identical with that of a 2-spin $s$-wave Hamiltonian, while Phases II and I correspond to 1-spin and 0-spin solutions, respectively. In Phase III, this technique does not work for the separable BCS models. In Appendix C, we show that the 2-spin solution for these nonintegrable models is identical to that of the integrable case, up to a rescaling of time, while the general asymptotic solution that we observe is Eq. (46). Thus, if a reduction mechanism exists in the nonintegrable cases, the form of the $m$-spin Hamiltonian must also change.

V.3.2 Relaxation time

In Sect. V.2 we saw that there are examples of nonintegrable separable BCS models where the constant $\Delta_{\infty}$ of Phase II is unstable to harmonic perturbations, and in Sect. 8 we gave evidence that the Phase III oscillations of these models are elliptic functions. This behavior is typical of integrable models as well, although the form of the elliptic functions changes once integrability is broken. A more important difference, however, is that a long relaxation time scale $\tau$ emerges before the system truly reaches Phase III.

Fig. 9 gives an example of the long relaxation time in the $d+id$ model, which is the separable BCS model with $f(\varepsilon)=\varepsilon$. The initial dynamics at weak coupling seem to indicatedidShortTime that $|\Delta(t)|$ oscillates with a single frequency reminiscent of Phase III. Upon closer inspection, however, the amplitude of the oscillations slowly changes with no indication of stabilizing. In Fig. 10, quenches at higher energies provide further evidence that the long-time asymptotic state is difficult to determine based on the short-time dynamics.

Let us now explore the dependence of the relaxation time $\tau$ on $\Delta_{0i}$, $\Delta_{0f}$ and $\gamma$ in the Lorentzian separable BCS model defined in Eq. (23). We define $\tau$ as the minimum time after which the minimum of $|\Delta(t)|$ oscillations stays within $\eta=10^{-4}$ of its asymptotic value. This definition of $\tau$ and the precise value of $\eta$ are somewhat arbitrary, but empirically we find that the minima of $|\Delta(t)|$ take longer to relax to the stationary value than the maxima. Typically, the minimum will increase for a time until it begins to oscillate with decreasing amplitude about a final value. Most importantly, this definition of $\tau$ delineates clearly the difference between integrable and nonintegrable behavior. Fig. 11 shows the dependence of $\tau$ on the values of $\Delta_{0i}$ and $\Delta_{0f}$, with one or the other fixed. Generally speaking, we find that quenches at lower energy scales increase $\tau$.

More interesting is the dependence of $\tau$ on $\gamma$, the integrability-breaking parameter, at fixed $(\Delta_{0i},\Delta_{0f})$. First, let us examine quenches that lead to Phase III in both the Lorentzian and integrable $s$-wave models. Fig. 12 shows that $\tau$ has single minimum for $\gamma\sim 0.4W$ and increases away from this point both as $\gamma\to 0$ and as $\gamma\to\infty$. In all cases, the relaxation time of quenches in the integrable $s$-wave model, which is the $\gamma\to\infty$ limit of our separable BCS Hamiltonians, is far smaller. We believe that the increase of $\tau$ as $\gamma\to\infty$ is indicative of nonperturbative behavior of the dynamics in the vicinity of the integrable limit, see Sect. VIII.

The behavior of $f(\varepsilon)$ as $\gamma\to 0$ is model dependent; in the case of the Lorentzian model, the stability analysis of Sect. V.2 indicates that Phase II is unstable to harmonic perturbations if $\gamma>\Delta_{\infty}$; otherwise, Phase II could be stable. We observe in Fig. 12 large oscillations in the evolution of the minimum of $\Delta(t)$ at $\gamma=0.2W$, behavior which occurs in the range $0.13W\lesssim\gamma\lesssim 0.26W$ For $\gamma\lesssim 0.13W$, the minima oscillations disappear and $\tau$ begins to dramatically increase. Despite this qualitative change in the evolution of $|\Delta(t)|$, down to at least $\gamma=0.11W$ we still find that the system eventually enters Phase III with a reduced amplitude of oscillation.

Fig. 13 is similar to Fig. 12, except we now choose $\Delta_{0i}$ and $\Delta_{0f}$ such that the (integrable) $s$-wave model enters Phase II. The behavior of $\tau$ with respect to $\gamma$ is qualitatively similar, except there is no regime where the minimum of $\Delta(t)$ undergoes large oscillations.

The spin-orbit model also has a very long relaxation time to Phase III. In order to observe this asymptotic state, as is shown in Fig. 8, one must carefully choose model and quench parameters, otherwise $\tau$ is simply too large for our present numerical study.

VI Phase III asymptotic solution

We now explore the structure of the Phase III asymptotic state. First, we treat $\Delta(t)$ as a periodic external driving and show that there is always a periodic solution for the classical pseudospins (and auxiliary functions in the spin-orbit model), and then we provide evidence that the class of periodic $\Delta(t)$ that are also self-consistent are elliptic functions.

VI.1 External driving

In the separable BCS model, the mean-field dynamics can be described alternatively by a Gaussian wave function with complex Bogoliubov amplitudes $u_{j}(t)$ and $v_{j}(t)$

[TABLE]

where normalization requires $|v_{j}|^{2}+|u_{j}|^{2}=1$. The equations of motion for $u(t)$ and $v(t)$ follow from the time-dependent Schrödinger equation $i\frac{\partial}{\partial t}|\psi\rangle=\hat{H}|\psi\rangle$ applied to (47) with the mean-field Hamiltonian from (2),

[TABLE]

where we shifted the Hamiltonian by a constant $\hat{H}=\hat{H}_{f}-\sum_{j}\varepsilon_{j}$ in order to make it traceless. The mapping to the classical pseudospins is

[TABLE]

We shall discuss the nature of the asymptotic Phase III $\Delta(t)$ in terms of $v(t)$ and $u(t)$. To do so, consider first Eq. (48) with a periodic $\Delta(t)=\Delta(t+T)$ that is not necessarily self-consistent, which decouples each pair of $(u_{j},v_{j})$ from all the others. The abstract form of Eq. (48) is

[TABLE]

with

[TABLE]

where $\mathbf{u}$ and $\mathbf{v}$ are $m$-dimensional vectors, $\mathbb{A}$ is a constant real symmetric $m\times m$ matrix, $\mathbb{B}(t)$ is a complex $m\times m$ matrix periodic in $t$ with period $T$, and we dropped the index $j$ for simplicity. The forthcoming discussion is valid for all systems of this form, see also Ref. ydgf, . For example, the spin-orbit dynamics admit such a representation with $m=4$, while $m=1$ in the separable BCS model.

As $\mathbf{h}(t)$ is periodic by assumption, the Floquet theorem applies. There are thus $2m$ independent solutions $\boldsymbol{\psi}_{i}(t)$ to Eq. (50) of the form

[TABLE]

where the $\mathbf{U}_{i}(t)$ and $\mathbf{V}_{i}(t)$ are periodic with the same period as $\mathbf{h}(t)$ and the $\delta_{i}$ are complex numbers known as Floquet exponents. The solutions $\boldsymbol{\psi}_{i}(t)$ therefore have the property

[TABLE]

where the $\rho_{i}$ are known as Floquet multipliers. Because $\mathbf{h}(t)$ is Hermitian, Eq. (50) conserves the norm of the solutions $\boldsymbol{\psi}_{i}(t)$, which implies $|\rho_{i}|=1$ and $\delta_{i}=i\nu_{i}$ for $\nu_{i}$ real. Furthermore, the particular form of $\mathbf{h}(t)$ implies that if $\boldsymbol{\psi}=(\mathbf{u},\mathbf{v})^{T}$ is a solution then so is $\widetilde{\boldsymbol{\psi}}=(\mathbf{v}^{*},-\mathbf{u}^{*})^{T}$. This pairing of solutions implies that if $\delta_{i}$ is a Floquet exponent, then so is $-\delta_{i}$. In Sect. VI.2, we will use this latter fact to prove that there is always a periodic spin solution to Eq. (6) for a given periodic $\Delta(t)$.

Before continuing, we note that the Phase III asymptotic $\Delta(t)$ is only periodic in the particle-hole limit of the separable BCS model. In the general case, $\Delta(t)=F(t)e^{-2i\mu_{\infty}t}$, where $F(t)$ is periodic. Nonetheless, we can still reduce this problem, where $\mathbf{h}(t)$ is not periodic, to the periodic case by absorbing the phase $2\mu_{\infty}t$ in the following manner:

[TABLE]

so that the time evolution of $(\mathbf{u}^{\prime},\mathbf{v}^{\prime})^{T}$ is described by Eq. (50) with periodic $\mathbf{h}(t)$ of the form given in Eq. (51) where $\mathbb{A}$ is replaced by $\mathbb{A}^{\prime}$. In terms of the pseudospin representation of the dynamics, this transformation amounts to an overall time-dependent rotation about the z-axis with frequency $2\mu_{\infty}$.

VI.2 Phase III spin solution in the separable BCS model

Now we draw our attention to the behavior of the spin solutions to the separable BCS model for the periodic external $\Delta(t)$ considered in the previous section. The dimension of the matrix $\mathbf{h}(t)$ is now $2m=2$ and there are two independent solutions to the Floquet problem

[TABLE]

where $U_{j}(t)$ and $V_{j}(t)$ are periodic and we restored the index $j$. Using $\boldsymbol{\psi}_{1j}(t)$ and Eq. (49), we can construct a periodic spin solution $\boldsymbol{\sigma}_{j}(t)$ [i.e., a periodic solution of Eq. (6) for the given external $\Delta(t)$ that does not necessarily satisfy Eq. (7)],

[TABLE]

We will now show that the most general spin solution $\mathbf{s}_{j}(t)$ precesses about the periodic solution $\boldsymbol{\sigma}_{j}(t)$ with a variable angular velocity. First we write the general solution $\boldsymbol{\Psi}_{j}(t)$ as a linear combination of $\boldsymbol{\psi}_{1j}(t)$ and $\boldsymbol{\psi}_{2j}(t)$

[TABLE]

Although the coefficients of linear combination are in principle complex, we can drop the constant overall phase of $\boldsymbol{\Psi}_{j}(t)$ as well as absorb $\frac{1}{2}\times$the remaining constant relative phase into the definitions of $U_{j}(t)$ and $V_{j}(t)$. Once again using (49), we now write $\boldsymbol{\Psi}_{j}(t)$ in terms of spin variables. It is helpful to first parametrize $U_{j}(t)$ and $V_{j}(t)$ as

[TABLE]

whence

[TABLE]

Note that $\theta_{j}$ is the only time-independent quantity in Eq. (58). A geometric interpretation of the motion of the general solution $\mathbf{s}_{j}(t)$ with respect to the periodic solution $\boldsymbol{\sigma}_{j}(t)$ becomes clear once we use Eq. (58) to express $\mathbf{s}_{j}(t)$ in the body coordinate system of $\boldsymbol{\sigma}_{j}(t)$. Let $\hat{\mathbf{z}}^{\prime}_{j}=\hat{\boldsymbol{\sigma}}_{j}$, while $\hat{\mathbf{x}}^{\prime}_{j}$ lies along the line defined by the intersection of the plane spanned by {$\hat{\mathbf{z}}^{\prime}_{j},\hat{\mathbf{z}}_{j}$} and that perpendicular to $\hat{\mathbf{z}}^{\prime}_{j}$. Finally, $\hat{\mathbf{y}}^{\prime}_{j}$ satisfies $\hat{\mathbf{y}}_{j}^{\prime}\cdot\hat{\mathbf{x}}_{j}^{\prime}=\hat{\mathbf{y}}^{\prime}_{j}\cdot\hat{\mathbf{z}}^{\prime}_{j}=0$ and $\hat{\mathbf{x}}_{j}^{\prime}\times\hat{\mathbf{y}}_{j}^{\prime}=\hat{\mathbf{z}}_{j}^{\prime}$. These definitions lead to

[TABLE]

The general spin solution $\mathbf{s}_{j}(t)$ in this new coordinate system is then

[TABLE]

where $\boldsymbol{\sigma}_{j}\cdot\boldsymbol{\sigma}_{j\perp}=0$ and $\boldsymbol{\sigma}_{j\perp}$ is not periodic. We see from Eq. (60) that $\mathbf{s}_{j}(t)$ makes a constant angle $\theta_{j}$ with the periodic solution and rotates about it with a variable angular frequency $\dot{\alpha}_{j}(t)$. From Eq. (57) and the periodicity of $U_{j}(t)$ and $V_{j}(t)$, we conclude that $\alpha_{j}(t)-2\nu_{j}t$ is also periodic with the same period as the external $\Delta(t)$ driving the system.

VI.3 Asymptotic self-consistency

Thus far, we have considered $\Delta(t)$ to be an external periodic driving that is not necessarily self-consistent. We showed for any such external driving, there is a corresponding periodic spin solution $\boldsymbol{\sigma}_{j}(t)$ with the same period as $\Delta(t)$. Furthermore, we derived in Eq. (60) that the general spin solution $\mathbf{s}_{j}(t)$ precesses in a simple manner about $\boldsymbol{\sigma}_{j}(t)$. In the true quench dynamics, however, $\Delta(t)$ must be self-consistent, and we now show that this requirement implies that there always exists a set of constants $\theta_{j}$, such that the following integral equation holds for the asymptotic periodic $\Delta(t)$:

[TABLE]

The notation $\boldsymbol{\sigma}_{j}=\boldsymbol{\sigma}_{j}[\Delta]$ emphasizes that the periodic spin solution is some complicated nonlocal function of $\Delta(t)$. An analogous expression to Eq. (61) exists for the spin-orbit model.

Eq. (61) is simply asymptotic self-consistency, as introduced in Sect. V, applied to the Floquet problem studied in Sects. VI.1 and VI.2. To see this, suppose that we observe some Phase III asymptotic periodic $\Delta(t)$ after a quench from the ground state of the separable BCS model, as discussed in Sect. V.3. This $\Delta(t)$ is self-consistent by definition, i.e.,

[TABLE]

which we write in terms of the underlying periodic spin solution $\boldsymbol{\sigma}_{j}$ by using Eq. (58)

[TABLE]

where $\nu_{j}$ is the imaginary part of the Floquet exponent as introduced in Eq. (52). As in our analysis of self-consistency in Phases I and II, Eq. (63) cannot hold exactly, this time because the sum over $\sigma_{j\perp}^{-}[\Delta]$ is the only non-periodic term. Nonetheless, under the reasonable assumption that $\nu_{j+1}-\nu_{j}\sim\delta$, where $\delta$ is the level spacing, the sum over $\sigma_{j\perp}^{-}[\Delta]$ dephases in $N\to\infty$ limit as $t\to\infty$ (the $N\to\infty$ limit comes first), leading to Eq. (61).

VI.4 Self-consistent solutions in the separable BCS model

We have seen that an asymptotically self-consistent periodic $\Delta(t)$ satisfies the functional equation (61) in the separable BCS model. We now will give evidence that solutions to Eq. (61) are elliptic functions. In order to generate such solutions, fix a period $T$ and write $\Delta(t)$ as a Fourier series

[TABLE]

which we truncate to some $n_{\textrm{max}}$, such that $c_{n}=0$ if $|n|>n_{\textrm{max}}$. In the particle-hole symmetric limit, $\Delta(t)$ is a real quantity that satisfies $\Delta(t)=\Delta(-t)$ [see Eq. (19)], so that $c_{n}$ is real and equals $c_{-n}$.

For a fixed set of coefficients $c_{n}$, we determine $\sigma^{x}_{j}[\Delta(t)]$ by solving the equations of motion (6) from $t=0$ to $t=T$ with $\Delta(t)$ given by (64). If the choice of $c_{n}$ produces a self-consistent $\Delta(t)$, then it will be equal to the quantity $\Delta_{\textrm{comp}}(t)$ defined as

[TABLE]

for some set of $\theta_{j}$. For most choices of $c_{n}$, however, Eq. (65) will not hold. As both $\Delta(t)$ and $\Delta_{\textrm{comp}}(t)$ are periodic functions of time with the same period, we define a distance $r(\{c_{n}\})$ as

[TABLE]

A given $\Delta(t)$ is asymptotically self-consistent if and only if $r(\{c_{n}\})=0$.

We now explore the results of this procedure for the Lorentzian coupling of the separable BCS model for various values of the integrability breaking parameter $\gamma$. It turns out that this procedure works when we fix $\cos\theta_{j}=1$, i.e., we find exactly (and not just asymptotically) self-consistent solutions. In order to find such solutions, we start from the known values of the Fourier coefficients of the $s$-wave ($\gamma=\infty$) solution, which are close to the Fourier coefficients of the $\gamma\gg 1$ solutions. We then progressively lower $\gamma$ while finding Fourier coefficients that minimize $r(\{c_{n}\})$. Typically we obtain values of $r\sim 10^{-12}-10^{-11}$ before declaring the solution self-consistent.

Fig. 14 gives of examples of such solutions at fixed $\Delta_{0f}$ and period $T$. Notably, there is a minimum $\gamma=\gamma_{\textrm{min}}$ below which the amplitude of oscillation vanishes. As $\gamma$ is increased from this minimum, the amplitude of oscillations increases to a maximum and then decreases to a nonzero limiting value as $\gamma\to\infty$. Fig. 14 also shows the fast convergence of $\gamma_{\textrm{min}}$ as a function of $N$ for two examples of this procedure.

In Sect. V.3, we argued through example quenches that the $\Delta(t)$ of Phase III are always elliptic functions, i.e., they satisfy Eq. (44). We show in Fig. 15 that the exactly self-consistent $\Delta(t)$ from Fig. 14 also satisfy Eq. (44) to a high degree of accuracy. The Floquet analysis of the equations of motion from Sect. VI.1 applies to any periodic $\Delta(t)$. From Fig. 15, we conclude that the self-consistency requirement (61) is essential to selecting elliptic functions amongst all possible periodic functions.

VII Quasiperiodic Phase IV

Quenches that do not conform to Phases I-III are another intriguing consequence of integrability breaking. We present two such examples in Figs. 16 and 17. Figs. 16a and 17a show a particle-hole symmetric quench of the separable BCS Hamiltonian with sine coupling from Eq. (24). Figs. 16b and 17b depict a quench of the Zeeman field in the spin-orbit model (14). The quasiperiodic behavior of $\Delta(t)$ in Fig. 16a sets in very early on, as corroborated by Fig. 17a, and persists with no appreciable changes at least until the times shown in the figure. Similarly, Fig. 16b is representative of the long-time spin-orbit $|\Delta(t)|^{2}$ as evidenced by Fig. 17b. Based on our preliminary analysis of the Fourier spectrum of $|\Delta(t)|^{2}$ for this quench and of the maximal Lyapunov exponent with the method of local divergence ratesHR , we believe that it too is quasiperiodic. However, a more careful study is needed to unambigously distinguish between quasiperiodicity and chaos in this case. Such a study is beyond the scope of the present paper, where we mainly focus on the properties of Phases I-III.

Note that the simulation times in Figs. 16 and 17 are enormous compared to the characteristic time of a single oscillation and even to typical Phase III relaxation times $\tau\Delta_{0f}\sim 10^{3}$ we observed in Sect. V.3.2, cf. Fig. 12 and the caption to Fig. 8. Thus, both of these examples do not belong to Phases I, II, or III. We therefore conclude that there are regions of quasiperiodicity in the quantum quench phase diagrams of nonintegrable pairing models, which we call Phase IV.

VIII Conclusion

The far-from-equilibrium steady states reached by nonintegrable pairing models after a quantum quench admit a similar taxonomy as do the integrable cases. We have shown that some or all of Phases I-III may occur in the separable BCS models and spin-orbit model defined in Eq. (1). The persistent periodic oscillations characterizing Phase III are always elliptic functions, regardless of whether the model is integrable. Moreover, we have developed a stability analysis of the three phases, summarized in Eq. (42), which generalizes known results in the integrable cases and elucidates the mechanism of nonequilibrium phase transitions using the language of linear analysis.

Despite these striking similarities, important consequences accompany integrability breaking. As argued in Sect. V.2, some nonintegrable models may not exhibit all three phases. At the same time, an entirely new quasiperiodic Phase IV emerges in certain models. Another key byproduct of integrability breaking is the emergence of a new, extremely long relaxation time scale $\tau$ when the asymptotic state either is or appears to be Phase III. For $t<\tau$, $\Delta$ can oscillate with more than one fundamental frequency and a slowly varying amplitude. This time scale is a generic feature of nonintegrable models, and its existence renders short-time analyses inadequate for determining the long-time dynamics. Moreover, $\tau$ diverges as we approach integrable points (e.g., as $\gamma^{-1}\to 0$ in the separable pairing models of Sect. V), and it is often too large for the practical determination of the true asymptotic state.

While the squared modulus of $\Delta(t)$ [and $\Delta(t)$ itself in the particle-hole symmetric case] is always an elliptic function in Phase III, its parametrization is more complicated in nonintegrable models. As a result, the reduction mechanism discussed in Appendix C, which explains how Phase III manifests itself in the integrable models, does not apply to nonintegrable models. Nonetheless, we demonstrated in Sect. VI that the common structure of the nonintegrable models implies the existence of a periodic solution to the classical pseudospin equations of motion if $\Delta(t)$ is taken to be a generic periodic external driving. Using numerical examples, we argued that further requiring $\Delta(t)$ to be self-consistent selects elliptic functions amongst all possible periodic functions.

It is instructive to discuss the BCS quench dynamics in terms of bifurcation theoryIntro_bifurc ; Kuznetsov ; Hilborn . For example, consider the particle-hole symmetric separable BCS models with real $\Delta$. For fixed initial conditions (17) and any function $\Delta(t)$ with fixed $\Delta(0)$, the equations of motion (6) have a unique solution ${\bf s}_{j}[\Delta(t)]\equiv{\bf s}[\varepsilon_{j},\Delta(t)]$. Eq. (7) then provides a closed nonlinear integral equation for $\Delta(t)$ [cf. Eq. (61)],

[TABLE]

Phase I is a fixed point, $\Delta=0$, of this equationnoteSep20_2 , while Phase II corresponds to two fixed points $\Delta_{\infty}$ and $e^{i\pi}\Delta_{\infty}=-\Delta_{\infty}$. In Phase III we end up on one of two limit cycles related to each other by a rotation by $\pi$ around the z-axis [change of sign of $\Delta(t)$]. The Phase I to II and II to III transitions correspond to supercritical pitchfork and Hopf bifurcations, respectively, in this languagenoteSep20_1 . The same results apply to the spin-orbit model (14). We also note that this quantum quench phase diagram is surprisingly similar to the nonequilibrium phase diagram of two atomic condensates coupled to a heavily damped cavity modeaniket1 ; aniket2 . The mean-field dynamics of the latter system are described by the driven-dissipative variant of the Bloch equations (6) for two classical spins representing individual condensates. Moreover, there are islands of quasiperiodicity in the phase diagram of the two coupled condensates, where the dynamics are very similar to that shown in Figs. 16 and 17.

Bifurcation theory also offers a plausible explanation for the divergence of the relaxation time $\tau$ near integrable points. Consider Phase III for an integrable pairing Hamiltonian, such as the particle-hole symmetric $s$-wave BCS model. Suppose the corresponding limit cycle loses stability as soon as integrability is broken and another limit cycle emerges as an attractor. An example of such behavior is the transcritical bifurcationIntro_bifurc ; Kuznetsov ; Hilborn . Because the instability is weak for weak integrability breaking and because the evolution starts near the unstable limit cycle, the system takes a very long time $\tau$ to reach the attractor. The weaker the integrability breaking, the closer we are to the bifurcation and the longer the time $\tau$.

An interesting open problem is to explore the newly discovered quasiperiodic Phase IV. In particular, one needs to investigate the possibility that asymptotic oscillations of $|\Delta(t)|$ for certain quenches may be chaotic, rather than quasiperiodic, i.e., the potential existence of a chaotic phase in addition to the quasiperiodic one. Let us also mention that quasiperiodic $|\Delta(t)|$ also occurs in integrable models, but only when the initial (pre-quench) state is a highly excited state instead of the ground statemulti .

In this paper, we employed reduced BCS Hamiltonians (1) to model pairing dynamics. This description is valid only at times $t\ll\Gamma^{-1}$, where $\Gamma$ is the highest among the rates of processes such Hamiltonians neglect. These processes include pair-breaking collisionsvolkov ; galperin ; barlev0 ; ydgf , three-body losses in ultracold gasesmatt , thermal fluctuationssasha , etc. Thus, to reach the asymptotic state before these effects influence the dynamics, we need $\Gamma^{-1}\gg\tau$. In Phases II and III, this requirement is much more stringent than $\Gamma^{-1}\gg T_{\Delta}$ typically quoted in the literature on collisionless pairing dynamics. Here $T_{\Delta}$ is the characteristic period of $\Delta(t)$ oscillations ($T_{\Delta}$ is of the order of the inverse equilibrium gap $\Delta_{0f}$ in our separable BCS models). Another limitation is the parametric instability of Phase III with respect to spontaneous eruptions of spatial inhomogeneitiesturbulence ; mattSep20 ; chern ; comment . To avoid this instability, the system size has to be smaller than the superconducting coherence length.

Acknowledgements.

We thank M. Dzero, A. J. Millis and A. Patra for helpful discussions. This work was supported by the National Science Foundation Grant DMR-1609829. J.A.S. was supported by a Rutgers University Louis Bevier Dissertation Completion Fellowship.

Appendix A Mean-field equations of motion

The pseudospin equations of motion for the separable BCS model (6) obtain simply from the Heisenberg equations of motion $\frac{d}{dt}\hat{A}=i\,[\hat{H},\hat{A}]$ applied to the mean-field $\hat{H}_{f}$ in Eq. (2) and the pseudospin operators $\hat{\mathbf{s}}$ defined in Eq. (5). The classical spin variables $\mathbf{s}$ are the expectation values of the pseudospin operators $\mathbf{s}=\langle\hat{\mathbf{s}}\rangle$, and the time-dependent order parameter $\Delta$ is determined self-consistently according to Eq. (7).

The generalized pseudospin representation of the spin-orbit Hamiltonian $\hat{H}_{so}$ from Eq. (1) requires more workdky . First, we diagonalize the kinetic part of $\hat{H}_{so}$ through the following unitary transformation to new fermionic operators $\hat{a}_{\mathbf{k}\pm}$

[TABLE]

where $\mathbf{k}=ke^{i\theta_{\mathbf{k}}}$ and $\phi_{k}$ is defined in terms of the model parameters in Eq. (3). One can check that the new elementary excitation energies are $\varepsilon_{\mathbf{k}\pm}\equiv\varepsilon_{\mathbf{k}}\mp R_{\mathbf{k}}$. Eq. (68) implies

[TABLE]

Upon summing over $\mathbf{k}$, the last two terms in parentheses cancel with momenta of opposite sign. Therefore, the interaction term of (1) in this new basis becomes

[TABLE]

and upon taking the mean-field approximation $\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\hat{c}\hat{c}\approx\langle\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\rangle\hat{c}\hat{c}+\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\langle\hat{c}\hat{c}\rangle-\langle\hat{c}^{{\dagger}}\hat{c}^{{\dagger}}\rangle\langle\hat{c}\hat{c}\rangle$, the interaction term becomes

[TABLE]

Neglecting the constant term $\Delta^{*}\Delta/g$, we arrive at the mean-field spin-orbit Hamiltonian $\hat{H}_{so}$ in the $\hat{a}$ basis found in (2). Similar to the separable BCS model, we now search for a set of quadratic fermionic operators whose equations of motion are closed. Define the following operators

[TABLE]

where $\eta_{\mathbf{k}}=e^{i\,\theta_{\mathbf{k}}}=-\eta_{-\mathbf{k}}$ and, as usual, $\hat{S}^{-}=\hat{S}^{x}-i\,\hat{S}^{y}$ and $\hat{L}^{-}=\hat{L}^{x}-i\,\hat{L}^{y}$.

One can check that $\hat{\mathbf{S}}_{\mathbf{k}\lambda}$, $\hat{\mathbf{L}}_{\mathbf{k}\lambda}$ and $\hat{T}_{\mathbf{k}}$ are Hermitian operators. There is reflection symmetry in $\mathbf{k}$-space: $\hat{A}_{-\mathbf{k}\lambda}=\hat{A}_{\mathbf{k}\lambda}$ for all operators $\hat{A}_{\mathbf{k}\lambda}$ in (72), as well as the following band symmetry for $\hat{\mathbf{L}}_{\mathbf{k}\lambda}$: $\hat{L}_{\mathbf{k}+}^{-}=\hat{L}_{\mathbf{k}-}^{-}$ and $\hat{L}_{\mathbf{k}+}^{z}=-\hat{L}_{\mathbf{k}-}^{z}$.

We apply the Heisenberg equations of motion to (72) and $\hat{H}_{so}$ from (2) and then take expectation values to arrive at the generalized pseudospin equations of motion (11). The time-dependent order parameter $\Delta=\langle\hat{\Delta}\rangle$ as a function of the new variables can be found in Eq. (12). The factor $\eta_{\mathbf{k}}$ does not appear in Eq. (11), implying that the dynamics preserve any radial symmetry found in the initial state. As all initial states considered in this work are radially symmetric, one can opt to label the generalized pseudospin variables by their single-particle energies rather than their momentum vector.

Appendix B Integrable limit of spin-orbit quenches

The authors of Ref. ddgp, created a full nonequilibrium phase diagram of the spin-orbit model for quenches of the magnetic field $h_{i}\to h_{f}$ as a function of $h_{i}$ and $h_{f}$. However, this phase diagram needs to be revised by running simulations to much longer times $t>\tau$, which, in particular, may modify the Phase II-III boundarynote22 . The phase diagram of Ref. ddgp, is also missing the quasiperiodic Phase IV discovered in the present work.

In Ref. dky, , an attempt was made to analyze interaction and external field quenches to the integrable limit $h_{f}=0$, but mistakes led to an incorrect phase diagram for the interaction quenches. Here we correct those mistakes and generate a correct phase diagram.

When the external field $h$ is set to zero, $H_{so}$ from (14) becomes equivalent to the integrable $s$-wave model with a dispersion relation $\varepsilon_{\mathbf{k}\lambda}=\frac{k^{2}}{2}-\lambda\alpha k$. This becomes clear in the equations of motion (11) with $\cos\phi_{k}=0$ and $\sin\phi_{k}=1$, where the spin degrees of freedom $\mathbf{S}_{\mathbf{k}\lambda}$ decouple from the others and $\Delta$ depends only on $\mathbf{S}_{\mathbf{k}\lambda}$. In what follows, the initial state of the system will be the ground state for some $h_{i}\geq 0$ given by (20), and the Hamiltonian for $t\geq 0$ is

[TABLE]

We use the integrability of $H$ to construct the exact phase diagram using a technique imported from Refs. ydgf, which we now summarize briefly. The analysis centers around a quantity $\mathbf{L}(u)$ called the Lax vector (not to be confused with the variables $\mathbf{L}_{\mathbf{k}\lambda}$)

[TABLE]

The integrability of $H$ follows from the fact that $L^{2}(u)$ is conserved by the time evolution for arbitrary $u$, which implies conservation of the $2N$ roots of $L^{2}(u)$, which we call $u_{j}$. As demonstrated in Ref. ydgf, , each of the asymptotic nonequilibrium phases corresponds a unique number of isolated complex pairs of $u_{j}$ in the continuum limit. Phase I corresponds to zero isolated $u_{j}$, Phase II corresponds to one pair, and Phase III corresponds to two pairs. As the $u_{j}$ are constants of the motion, we can evaluate $L^{2}(u)$ at $t=0$ to determine the number of isolated pairs of $u_{j}$ and thus generate the phase diagram for a given $h_{i}$.

Let us first start with the case when $h_{i}=0$ and we quench the interaction $g_{i}\to g_{f}$. In this case the ground state self-consistency relationship is

[TABLE]

Using Eq. (75) along with the initial state given by Eq. (17), we find that the initial Lax vector has the form

[TABLE]

If $g_{f}=g_{i}$, i.e., the zero quench, then $\widetilde{\beta}=0$ and the only complex pair of roots is $u_{\pm}=\pm i\,\Delta_{0i}+\mu$. This is the degenerate Phase II case, where $\Delta(t)=\Delta_{0i}$ identically. When $g_{f}\neq g_{i}$, $L^{2}(u)=0$ implies

[TABLE]

We now construct the phase diagram shown in Fig. 18 for the $h_{i}=h_{f}=0$, $g_{i}\to g_{f}$ quenches in the spin-orbit model. As we will not utilize particle-hole symmetry, the chemical potential $\mu$ must be calculated from the fermion number Eq. (16), which in the present case reads

[TABLE]

In the continuum limit, we have the following translation from sums over $\mathbf{k}\lambda$ to integrals over the continuum for arbitrary functions $F(\varepsilon_{\mathbf{k}\lambda})$

[TABLE]

Thus, the spin-orbit coupling $\alpha$ at $h=0$ has the simple effect of introducing a peculiar density of states $\nu_{\alpha}(x)$ to the $s$-wave problem. Let $\widetilde{B}=\lim_{N\to\infty}\widetilde{\beta}/N$ and $n=\lim_{N\to\infty}N_{f}/N$, the latter of which is fixed for the entire phase diagram. For a given pair $(\Delta_{0f},\Delta_{0i})$, we first solve for $(\mu_{f},\mu_{i})$ and then for $\widetilde{B}$ through the following integral equations:

[TABLE]

We then use $\widetilde{B}$ and $\mu_{i}$ as input for the following integral equation:

[TABLE]

which we solve for $u$. The number of complex pairs of roots to Eq. (81) determines which nonequilibrium phase the system enters.

Quenches from $h_{i}\neq 0$ to $h_{f}=0$ still undergo integrable dynamics, except now the initial state is no longer the $s$-wave ground state. We consider the behavior of the zeros of $L^{2}(u)$ with respect to $h_{i}$ in the continuum limit with the spin-orbit parameters given in Fig. 2. The Lax vector is still as defined in Eq. (74), but we now enter the spin-orbit ground state (20) into the equation $L^{2}(u_{j})=0$, which implies $L^{x}(u_{j})=\pm i\,L^{z}(u_{j})$. The spin components of the $h_{i}\neq 0$ ground state are functions of the form $F_{\lambda}(\varepsilon_{\mathbf{k}})$ instead of $F(\varepsilon_{\mathbf{k}\lambda})$; we therefore do not use (79) for the continuum limit, but rather

[TABLE]

The result of the root calculation is given in Fig. 19, where we plot the absolute value of the imaginary part of each root pair. For small $h_{i}$, there is only one pair of complex roots, i.e., the asymptotic phase is Phase II. At a certain critical $h_{i}$, a second pair of complex roots appears, and the system enters Phase III. For larger $h_{i}$, the two pairs of roots merge into one and the system reenters Phase II. Phase I does not occur in $h_{f}=0$ quenches for the parameters we used.

Appendix C Integrability breaking forbids asymptotic reduction

An important property of the quench dynamics of integrable $s$ and $(p+ip)$-wave Hamiltonians is the dynamical reduction in the number of degrees of freedom at $t\to+\infty$ in the thermodynamic limitydgf ; fdgy . In particular, Phase  III in these models corresponds to the motion of two collective classical spins ${\bf S}_{1}$ and ${\bf S}_{2}$ governed by a Hamiltonian of the same form. The asymptotic order parameter $\Delta(t)$ in Phase III coincides with that of the $2$-spin problem. Further, there are special reduced solutions of equations of motion with the same $\Delta(t)$ that are of the form

[TABLE]

where $\alpha_{j}$, $\beta_{j}$ and $\eta_{j}$ are time-independent and $\hat{\bf z}$ is a unit vector along the z-axis. These observations lead to an analytical expression for $\Delta(t)$ and, moreover, help to construct the full asymptotic spin configuration in Phase III. We note also that, as we will see below, for the $s$-wave BCS model in the particle-hole symmetric case, (83) is equivalent to the ansatz of Ref. barlev0, .

We will now show that the above reduction mechanism relies on integrability and breaks down for nonintegrable separable BCS models. We will prove two independent statements: (i) reduced solutions exist only when $f^{2}(x)=C_{1}+C_{2}x$, i.e. only when the Hamiltonian is integrablerichardson ; ortiz2005 and (ii) $\Delta(t)$ for a 2-spin separable BCS Hamiltonian with an arbitrary choice of new $\varepsilon_{1,2}$, $f_{1,2}$ and $g$ does not match the asymptotic $\Delta(t)$ we obtained in Sect. V.3.

C.1 Existence of reduced solutions implies integrability and vice versa

We will follow the same steps as in the derivation of the 2-spin solutions in Ref. ydgf, and show that it only works for special choices of $f(x)$. First, we treat the general non-particle-hole symmetric case.

Let

[TABLE]

The 2-spin (reduced) Hamiltonian is

[TABLE]

where $\Delta=\widetilde{g}(\widetilde{f}_{1}S^{-}_{1}+\widetilde{f}_{2}S^{-}_{2})$. We take both $\widetilde{f}_{k}$ to be nonzero, because otherwise the two spins simply decouple and rotate uniformly around the z-axis.

Energy and $S^{z}_{1}+S^{z}_{2}$ are conserved. Since there are two conservation laws and two degrees of freedom, $H_{\mathrm{red}}$ is integrable. For more than two spins, integrability persists only for special choices of $\widetilde{f}_{k}$. This fact alone already distinguishes the 2-spin problem from that of a generic $N$-spin separable BCS Hamiltonian.

Conservation of energy and $S^{z}_{1}+S^{z}_{2}$ read

[TABLE]

We need $\widetilde{\varepsilon}_{1}\neq\widetilde{\varepsilon}_{2}$ or $|\Delta|$ will be constant. We use Eq. (86) to express $S^{z}_{k}$ in terms of $\Omega^{2}$,

[TABLE]

where $\widetilde{a}_{k}$ and $\widetilde{b}_{k}$ are time-independent and $\widetilde{a}_{1}=-\widetilde{a}_{2}\neq 0$. Furthermore, Eq. (83) implies a similar expression for $s_{j}^{z}$ in terms of the order parameter amplitude,

[TABLE]

Conservation of the energy

[TABLE]

and of $J_{z}=\sum_{j}s_{j}^{z}$ require

[TABLE]

We write the Bloch equations for the separable BCS Hamiltonian as

[TABLE]

Since the equations of motion and Eqs. (87) and (88) for the reduced solution and the 2-spin problem have the same form, we can treat both of them simultaneously.

Substituting Eq. (88) into Eq. (91), we find

[TABLE]

Next, we multiply Eq. (92) by $e^{i\Phi}$ and add the resulting equation to its complex conjugate,

[TABLE]

where we made use of Eq. (93). Integrating and adding the resulting equation and Eq. (93), we obtain

[TABLE]

where $\frac{a_{j}c_{j}}{f_{j}}$ is the integration constant and $A=\int dt\dot{\Phi}\dot{\Omega}$. The self-consistency condition $\Delta=g\sum_{j}f_{j}s_{j}^{-}$, combined with Eq. (90), implies $\sum_{j}a_{j}c_{j}=0$.

The analogous expressions for the 2-spin problem are

[TABLE]

and $\widetilde{a}_{1}\widetilde{c}_{1}+\widetilde{a}_{2}\widetilde{c}_{2}=\widetilde{a}_{1}(\widetilde{c}_{1}-\widetilde{c}_{2})=0$. Therefore, $\widetilde{c}_{1}=\widetilde{c}_{2}$ and the last term in Eq. (96) can be absorbed into $A$, which is defined up to a constant anyway, i.e.,

[TABLE]

Since $s_{j}^{-}$ is related to $S_{1}^{-}$ and $S_{2}^{-}$ via Eq. (83), this also eliminates the last term in Eq. (95), i.e.,

[TABLE]

Combining the conservation of the spin norm, $s_{j}^{2}=(s_{j}^{z})^{2}+|s_{j}^{-}|^{2}$, with Eqs. (88) and (98), we derive the following differential equation for $\Omega$:

[TABLE]

or, equivalently,

[TABLE]

This equation implies, among other things, that $A$ is a function of $\Omega$. Indeed, consider a set of numbers $x_{j}$, such that $\sum_{j}x_{j}=0$. Multiplying Eq. (100) by $x_{j}$ and summing over $j$, we find

[TABLE]

where $\lambda,\mu$ and $\kappa$ are real constants. Substituting this back into Eq. (100), we obtain

[TABLE]

where $w=\Omega^{2}$ and $\xi_{j}=\varepsilon_{j}-\mu$. These equations are consistent only when the coefficients of powers of $w$ are $j$-independent. In particular, we must have $f_{j}^{2}=4\lambda\xi_{j}+\mathrm{const.}$, i.e.,

[TABLE]

where $C_{1}$ and $C_{2}$ are real constants. This is the most general form of $f_{j}$ for which the separable BCS Hamiltonian (8) is known to be integrable richardson ; ortiz2005 . In particular, $C_{2}$=0 corresponds to the $s$-wave and $C_{1}=0$ to the $(p+ip)$-wave models. Conversely, when Eq. (103) holds and the separable Hamiltonian is therefore integrable, the $j$-independence of coefficients at $w^{2}$ and $w$ determines $a_{j}$ and $b_{j}$, and Eq. (102) means that $w=|\Delta|^{2}$ is a certain elliptic function of time.

C.2 Asymptotic $\Delta(t)$ does not match the 2-spin solution in nonintegrable cases

In Sect. V.3 we numerically determined $\Delta(t)$ in two nonintegrable separable BCS Hamiltonians, see Eq. (46). Here we show that $\Delta(t)$ for the most general separable 2-spin Hamiltonian (85) cannot match Eq. (46).

Since $\Delta(t)$ in Eq. (46) is real, we take $\Delta$ in the 2-spin problem to be real as well, though we do not a priori assume particle-hole symmetry in the 2-spin problem. All we need is to specialize the derivation of the previous subsection to the case of real $\Delta$. Then, the Bloch equations become

[TABLE]

Substituting Eq. (87) into the first two equations of motion, we obtain

[TABLE]

and

[TABLE]

where $\frac{\widetilde{a}_{k}\widetilde{c}_{k}}{\widetilde{f}_{k}}$ is the integration constant. As before, the self-consistency condition $\widetilde{g}(\widetilde{f}_{1}S_{1}^{-}+\widetilde{f}_{2}S_{2}^{-})=\Delta$ together with $\widetilde{a}_{1}=-\widetilde{a}_{2}$ imply $\widetilde{c}_{1}=\widetilde{c}_{2}\equiv\widetilde{c}$, and the conservation of spin length $(S_{k}^{x})^{2}+(S_{k}^{y})^{2}+(S_{k}^{z})^{2}=S_{k}^{2}$ yields

[TABLE]

Equating the coefficients at different powers of $\Delta$ for $k=1$ and 2, we find $\widetilde{c}(\widetilde{\varepsilon}_{1}-\widetilde{\varepsilon}_{2})=0\Rightarrow\widetilde{c}=0$,

[TABLE]

and two more relationships that constrain $\widetilde{a}_{k}$ and $\widetilde{b}_{k}$. The constraint (108) is a consequence of the requirement that $\Delta$ be real. Now Eq. (107) is of the form

[TABLE]

This is the same as the equation for the asymptotic $\Delta(t)$ for the integrable $s$-wave BCS Hamiltonian in the particle-hole symmetric case up to rescaling $\Delta_{\mathrm{new}}=\widetilde{f}\Delta$. This is not surprising because $\widetilde{f}_{1}=\widetilde{f}_{2}=\widetilde{f}$ and the factor of $\widetilde{f}^{2}$ in Eq. (85) can be absorbed into the coupling constant, $\widetilde{g}_{\mathrm{new}}=\widetilde{f}^{2}\widetilde{g}$ resulting in an integrable $s$-wave BCS Hamiltonian for two spins with $\Delta_{\mathrm{new}}=\widetilde{g}_{\mathrm{new}}(S_{1}^{-}+S_{2}^{-})=\widetilde{f}\Delta$. The solution of Eq. (109) is $\Delta(t)=\Delta_{+}\mathrm{dn}[\widetilde{f}\Delta_{+}(t-t_{0}),1-\frac{\Delta_{+}^{2}}{\Delta_{-}^{2}}]$. As we saw in Sect. V.3, in the nonintegrable case we find instead a more general differential equation Eq. (44) with the solution given by Eq. (46).

Appendix D The link between Lax constructions and the stability analysis

As mentioned above, the separable BCS model is integrable when $f_{j}^{2}=C_{1}\varepsilon_{j}+C_{2}$. Two important cases are the $s$-wave model where $f_{j}=1$ and the $p+ip$ model where $f_{j}=\sqrt{\varepsilon_{j}}$. In past workydgf ; fdgy , integrability has been exploited to determine the nonequilibrium asymptotic phases through the use of Lax constructions. These techniques are useful for constructing phase diagrams, but the physical interpretation of the phase transitions is obscured by the use of exact solvability. We demonstrate here that the stability equation Eq. (42), which applies to the nonintegrable cases as well, both predicts the same transition points and clarifies the physical meaning of the Lax construction.

In the following, we will assume the quantities $Z_{j}$, $\Delta_{\infty}$ and $\mu_{\infty}$ are given. They are functions of the quench parameters $\Delta_{0i}$, $\Delta_{0f}$, the particle number $N_{f}$, and the Fermi energy $\varepsilon_{F}$.

D.1 Lax norms

In the $s$-wave model, the Lax vector isydgf

[TABLE]

while in the $p+ip$ model its components arefdgy

[TABLE]

where $u$ is a complex (spectral) parameter.

We focus on the norms of these quantities, defined as $L^{2}(u)=L_{x}^{2}(u)+L_{y}^{2}(u)+L_{y}^{2}(u)$ in the $s$-wave case and $L_{2}(u)=uL^{+}(u)L^{-}(u)+[L^{z}(u)]^{2}$ for $p+ip$. Integrability follows from the fact that the $L^{2}(u)$ and $L_{2}(u)$ are conserved by the time evolution for arbitrary $u$, which implies conservation of their roots $u_{j}$. As demonstrated in Refs. ydgf, and fdgy, and discussed in Appendix B, each of the asymptotic nonequilibrium phases corresponds a unique number of isolated complex pairs of $u_{j}$ in the continuum limit. Phase I corresponds to zero isolated $u_{j}$, Phase II corresponds to one pair, and Phase III corresponds to two pairs.

The main result of this Appendix is that the roots of the Lax norm $u$ and the frequencies $\omega$ of $\delta\Delta(t)$ are related by $u-u_{r}=\pm\frac{1}{2}\sqrt{\omega^{2}-b_{\textrm{min}}^{2}}$, where $u_{r}$ is the real part of the root (cf. Refs. yta, ; ydgf, ), and $b_{\textrm{min}}$ is the band edge in the frequency spectrum ($b_{\textrm{min}}=0$ in Phase I). Thus, the new pair of complex conjugate Lax roots appears at the same time that $\omega$ emerges into the band gap (i.e., $\omega^{2}<b_{\textrm{min}}^{2}$ in Phase II and $\omega^{2}<0$ in Phase I). Here and below in this Appendix, we use the same notation $u$ for the roots and for generic values of the spectral parameter.

One may plug into the Lax norms the asymptotic spin solution (LABEL:eomHFphase2sol) for Phase II, but we shall use solutions that do not impose particle-hole symmetry. Letting $\widetilde{\varepsilon}_{j}=\varepsilon_{j}-\mu_{\infty}$, and noting that sums over the time-dependent terms dephase in the $t\to\infty$ limit, we find

[TABLE]

Eq. (112) reduces to the Phase I Lax norms when $\Delta_{\infty}=0$ and by convention $Z_{j}\to z_{j}$. In Phase II, Eq. (112) is supplemented by the self-consistency relationship

[TABLE]

D.2 Phase I-II transition

In the $s$-wave case, and in Phase I, we compare the stability equation Eq. (42) to the vanishing of the Lax norm $L^{2}(u)=0$. After some algebra, Eqs. (42) and $L^{2}(u)=0$ become

[TABLE]

respectively. We argued in Sect. V.2 that the Phase I-II transition occurs when a purely imaginary pair of complex conjugate $\omega_{0}$ emerges as solutions to Eq. (114a), implying an exponential instability to Phase I. The Lax construction stipulates that the same transition occurs when an isolated pair of complex conjugate $u$ solve Eq. (114b). In order for these two methods to match, we must make the identification $u-\mu_{\infty}=\pm\frac{1}{2}\omega_{0}$, i.e., the real part of the emergent Lax norm pair of roots must be $\mu_{\infty}$. We prove this is the case in Sect. D.4.

The corresponding equations for Phase I in the $p+ip$ model are

[TABLE]

and the same identification reconciles the two approaches.

D.3 Phase II-III transition

In Phase II, one applies the self-consistency relationship (113) to the Lax norms (112). In the $s$-wave case, $L^{2}(u)=0$ becomes

[TABLE]

and we see the single pair of isolated conjugate roots are $u_{\pm}=\mu_{\infty}\pm i\Delta_{\infty}$. The equation for the second pair of isolated roots that would signal a transition to Phase III is therefore

[TABLE]

After applying Eq. (113) to the quantities $S_{j}(\omega_{0})$ in the stability equation (42), we find for the $s$-wave model

[TABLE]

This simplifies Eq. (42) to

[TABLE]

Matching (119) to (117), we make the correspondence $u-\mu_{\infty}=\pm\frac{1}{2}\sqrt{\omega_{0}^{2}-4\Delta_{\infty}^{2}}$. As we discussed in Sect. V.2, an $\omega_{0}$ emerging out of the continuum and into the band gap signals the transition to Phase III. The band edge in the $s$-wave model is precisely $2\Delta_{\infty}$. We show in Sect. D.4 that the new pair of conjugate Lax roots has real part $\mu_{\infty}$. Therefore, the two approaches predict the same phase transition.

In the $p+ip$ case, $L_{2}(u)=0$ couples with (113) to give

[TABLE]

The single pair of isolated roots of Phase II is then

[TABLE]

and the emergent pair of conjugate roots solves

[TABLE]

To show that the stability analysis reproduces Eq. (122), we will need two relations. The first holds in general by applying the self-consistency relation (113) to the sums in (42)

[TABLE]

while the second is specific to the $p+ip$ model

[TABLE]

We substitute Eqs. (123)-(124) into Eq. (42), which becomes a quadratic function of $S_{3}$ and $S_{4}$. The solution is

[TABLE]

where $B_{1}=\sqrt{4\mu_{\infty}\Delta_{\infty}^{2}-\Delta_{\infty}^{4}}$ is the band edge when $u_{c}\geq 0$. In this parameter range, we identity $u-u_{c}=\pm\frac{1}{2}\sqrt{\omega^{2}-B_{1}^{2}}$. We show in Sect. D.4 that the real part of the emergent Lax roots is $u_{c}$, and therefore the stability analysis and Lax constructions give the same Phase II-III transition. When $u_{c}<0$, the band edge is no longer $B_{1}$, and we believe there to be no Phase II-III transition in that case.

D.4 Real parts of Lax roots at the transitions

The equivalence between the Lax construction and the stability analysis relies on the fact that the real parts of the emerging Lax roots are equal to $\mu_{\infty}$ at the Phase I-II transition in both integrable models, $\mu_{\infty}$ at the Phase II-III transition in the $s$-wave model, and $\mu_{\infty}-\frac{\Delta_{\infty}^{2}}{2}$ at the Phase II-III transition in the $p+ip$ model. In other words, the emergent second pair of isolated roots has the same real part as the first pair of isolated roots.

The Phase I-II transition real parts can be understood by a continuity argument. In the $s$-wave model, Eq. (116) implies that the single pair of roots can be written as $u_{\pm}=\mu_{\infty}\pm i\Delta_{\infty}$. As we approach the I-II boundary, $\Delta_{\infty}$ decreases continuously to zero, which implies the real part of both roots at the boundary is $\mu_{\infty}$. In the $p+ip$ case, a similar argument follows from Eq. (121).

D.4.1 $s$-wave, II-III

We use results from the spin reduction mechanism, discussed in Appendix C, of the $s$-wave model to obtain the real parts of the Lax roots at the Phase II-III transition. This discussion quotes several results directly from Sect. II B 3 of Ref. ydgf, . The isolated roots in Phase III of $L^{2}(u)$ are given by the roots of the 2-spin spectral polynomialydgf $Q_{4}(u)$

[TABLE]

We determine the real parameters $\mu$, $\rho$, $\kappa$ and $\chi$ at the transition, which will then give the roots of $Q_{4}(u)$. To do so, we use the differential equation and solution for the 2-spin $\Delta$, which is identical to that of the Phase III asymptotic $\Delta$ of the many-body problem, which we write as $\Delta=|\Delta|e^{-i\Phi}$. Let $w=|\Delta|^{2}=\Lambda^{2}+h_{1}$, where $h_{1}$ is a constant. The differential equation for $w$ is

[TABLE]

while the equation for the phase $\Phi$ is

[TABLE]

Upon rewriting (127) as an equation for $\Lambda$, we find

[TABLE]

where the constants $\Lambda_{\pm}$ are the maximum and minimum of the $\Lambda$ oscillations which are functions of the constants $\rho$, $\chi$ and $\kappa$. The solution of interest to Eq. (129) is

[TABLE]

Near the II-III transition, the oscillations of $\Lambda$ are small and it sufficient to keep only the first harmonic of Eq. (130)

[TABLE]

As we approach the II-III transition, $\Delta\to\Delta_{\infty}e^{-2i\mu_{\infty}t}$. Because $|\Delta|^{2}=\Lambda^{2}+h_{1}$ has the same frequency as $\Lambda^{2}$, and the frequency of small oscillations of $|\Delta|^{2}$ at the II-III transition is $2\Delta_{\infty}$, we conclude $\Lambda_{0}=\Delta_{\infty}$ and $h_{1}=0$. Using Eq. (128), we also find $\kappa=0$ and $\mu=\mu_{\infty}$.

It remains to determine the constants $\rho$ and $\chi$, which we do by plugging (131) into (127) and considering the $\mathcal{O}(\delta^{0})$ and $\mathcal{O}(\delta)$ terms separately. The result is $\rho=-\frac{\Delta_{\infty}^{2}}{2}$ and $\chi=\frac{\Delta_{\infty}^{4}}{4}$. The roots of the spectral polynomial $Q_{4}(u)$ from Eq. (126) at the Phase II-III transition therefore solve

[TABLE]

One solution to (132) is $u_{\pm}=\mu_{\infty}\pm i\Delta_{\infty}$, which is the single isolated pair characteristic of Phase II. The other solution is a double root at $u=\mu_{\infty}$, i.e., the new pair of roots that emerges in Phase III has real part $\mu_{\infty}$.

D.4.2 $p+ip$, II-III

In order to prove that the Lax construction and stability analysis predict the same $p+ip$ Phase II-III transition, we needed to assume that the real part of the emerging second pair of roots equals that of the first pair of roots $u_{\pm}$ from (121). Using results from Ref. fdgy, , we now show that this is indeed the case.

For brevity, our derivation will use the conventions of Ref. fdgy, , where the definitions of some quantities differ by numerical factors. One redefines $\varepsilon\to 2\varepsilon$, $2G\to g$, $\sqrt{2}\Delta\to\Delta$ and $u\to 2u$ in order to translate quantities from Ref. fdgy, to those in this work. While some details of the derivation depend on such conventions, the conclusion does not. We also assume $u_{c}\equiv\textrm{Re}[u_{\pm}]\geq 0$, which is the parameter regime where we show the equivalence of the Lax construction and stability analysis for the $p+ip$ model.

Eq. (4.3) of Ref. fdgy, gives the isolated pair of roots in Phase II to be $u_{\pm}=u_{c}\pm 2iE_{\mathsf{min}}$, where $E_{\mathsf{min}}$ is the minimum of the asymptotic dispersion relation [see text below Eq. (5.29) in Ref. fdgy, ]. According to Eq. (4.39) in Ref. fdgy, the frequency of small oscillations in Phase III close to the Phase II-III boundary is

[TABLE]

where $u_{\mathfrak{r}}$ is the real part of the pair of roots absent in Phase II. The frequency $\Omega_{c}$ should match the frequency of dephasing oscillations in Phase II close to the boundary. The text below Eq. (3.53) in Ref. fdgy, says that the latter frequency is

[TABLE]

Setting $\Omega_{c}=\Omega$, implies that on the Phase II-III boundary

[TABLE]