Quantum quench phase diagrams of an s-wave BCS-BEC condensate

Quantum quench phase diagrams of an -wave BCS-BEC condensate

E. A. Yuzbashyan, M. Dzero, V. Gurarie, and M. S. Foster Center for Materials Theory, Rutgers University, Piscataway, NJ 08854, USA
Department of Physics, Kent State University, Kent, OH 44240 USA
Department of Physics, University of Colorado, Boulder, CO 80309, USA
Department of Physics and Astronomy, Rice University, Houston, Texas 7700, USA

We study the dynamic response of an -wave BCS-BEC (atomic-molecular) condensate to detuning quenches within the two channel model beyond the weak coupling BCS limit. At long times after the quench, the condensate ends up in one of three main asymptotic states (nonequilibrium phases), which are qualitatively similar to those in other fermionic condensates defined by a global complex order parameter. In phase I the amplitude of the order parameter vanishes as a power law, in phase II it goes to a nonzero constant, and in phase III it oscillates persistently. We construct exact quench phase diagrams that predict the asymptotic state (including the many-body wavefunction) depending on the initial and final detunings and on the Feshbach resonance width. Outside of the weak coupling regime, both the mechanism and the time dependence of the relaxation of the amplitude of the order parameter in phases I and II are modified. Also, quenches from arbitrarily weak initial to sufficiently strong final coupling do not produce persistent oscillations in contrast to the behavior in the BCS regime. The most remarkable feature of coherent condensate dynamics in various fermion superfluids is an effective reduction in the number of dynamic degrees of freedom as the evolution time goes to infinity. As a result, the long time dynamics can be fully described in terms of just a few new collective dynamical variables governed by the same Hamiltonian only with “renormalized” parameters. Combining this feature with the integrability of the underlying (e.g. the two channel) model, we develop and consistently present a general method that explicitly obtains the exact asymptotic state of the system.

67.85.-d, 67.85.De, 34.90.+q, 74.40.Gh, 02.30.Ik, 05.45.Xt, 05.30.Fk, 05.45.Yv

I Introduction

The problem of a superconductor driven out of equilibrium by a sudden perturbation goes back many decades. Early studiesAnderson (1958); Galaiko (1972); Volkov and Kogan (1974); Galperin et al. (1981); Littlewood and Varma (1982); Shumeiko (1990) addressed small deviations from equilibrium using linearized equations of motion. An important result was obtained by Volkov and Kogan Volkov and Kogan (1974), who discovered a power law oscillatory attenuation of the Bardeen-Cooper-Schriffer (BCS) order parameter for non-equilibrium initial conditions close to the superconducting ground state.

In the last decade it was realized that even large deviations from equilibrium are within the reach of appropriate theoretical methods. Recent studies, motivated by experiments in cold atomic fermions, focused on quantum quenches – nonequilibrium conditions created by a sudden change in the superconducting coupling strength. Barankov et. al. Barankov et al. (2004) in a paper that set off a surge of modern research in this long-standing problem Amin et al. (2004); Andreev et al. (2004); Barankov and Levitov (2004); Szymanska et al. (2005); Warner and Leggett (2005); Yuzbashyan et al. (2005a, b, c); Yuzbashyan et al. (2006); Barankov and Levitov (2006); Yuzbashyan and Dzero (2006); Dzero et al. (2007); Barankov and Levitov (2007); Tomadin et al. (2008); Nahum and Bettelheim (2008); Gurarie (2009); Faribault et al. (2009) in the context of quantum gases, found that for initial conditions close to the unstable normal state, the order parameter exhibits large unharmonic periodic oscillations.

Subsequently, Yuzbashyan et. al.Yuzbashyan et al. (2006) developed an analytical method to predict the state of the system at large times based on the integrability of the underlying BCS model. This work extended Volkov and Kogan’s result to large deviations from equilibrium and showed that the oscillation frequency is twice the nonequilibrium asymptotic value of the order parameter, a conclusion confirmed by recent terahertz pump pulse experiments in films Matsunaga et al. (2013, 2014). Later studies Barankov and Levitov (2006); Yuzbashyan and Dzero (2006) mapped out the full quantum quench “phase diagram” for weakly coupled s-wave BCS superconductors finding that three distinct regimes occur depending on the strength of the quench: Volkov and Kogan like behavior, persistent oscillations, and exponential vanishing of the order parameter. Most recent researchPapenkort et al. (2007, 2009); Krull et al. (2013); Tsuji (2014) fueled by experimental breakthroughsMatsunaga and Shimano (2012); Beck et al. (2013); Matsunaga et al. (2013) investigates non-adiabatic dynamics of -wave BCS superconductors in response to fast electromagnetic perturbations. Closely related subjects developing in parallel are exciton dynamicsBrierley et al. (2011), collective neutrino oscillationsPehlivan et al. (2011); Raffelt et al. (2013), quenched p-wave superfluidsFoster et al. (2013, 2014) etc.

Most existing work addressed the dynamics in the BCS regime and, in particular, quenches such that the interaction strength is weak both before and after the quench. This was so that the system always remains in the BCS regime, since the physics of the condensate beyond this regime was not sufficiently well understood. However, a superfluid made up of cold atoms can be as well quenched from the BCS to the Bose-Einstein Condensation (BEC) regime, or within the BEC regime. With few exceptions Gurarie (2009); Foster et al. (2013, 2014), these types of quenches are not adequately studied in the existing literature.

Our paper aims to close this gap and analyze all possible interaction quenches throughout the BCS-BEC crossover in a paired superfluid, including BCS to BEC, BEC to BCS and BEC to BEC quenches. We fully determine the steady state of the system at large times after the quench – the asymptote of the order parameter as well as the approach to the asymptote, the many-body wave function, and certain observables, such as the radio-frequency absorption spectrum and the momentum distribution. In the BCS limit, we recover previous results. Beyond this limit the dynamics is quantitatively and sometimes qualitatively different. For example, the power law in the Volkov and Kogan like attenuation changes in the BEC regime, exponential vanishing is replaced with a power law, and persistent oscillations first change their form and then disappear altogether after a certain threshold for quenches from any initial (e.g. arbitrarily weak) to sufficiently strong final coupling. We believe an experimental verification of the predictions of this work is within a reach of current experiments in cold atomic systems.

The long time dynamics can be determined explicitly due to a remarkable reduction mechanism at work, so that at large times the system is governed by an effective interacting Hamiltonian with just a few classical collective spin or oscillator degrees of freedom. In a sense, the system “flows in time” to a much simpler Hamiltonian. This observation combined with the integrability of the original Hamiltonian (see below) lead to a method originally proposed in Ref. Yuzbashyan et al., 2006 for obtaining the long time asymptote (steady state) of integrable Hamiltonian dynamics in the continuum (thermodynamic) limit. Here we improve this method as well as provide its comprehensive and self-contained review including many previously unpublished results and steps. We do so in the context of the -wave BCS (one channel) and inhomogeneous Dicke (two channel) models, but with some modifications the same method also applies to all known integrable pairing models Gaudin (1983); Richardson (2002); Ibanez et al. (2009); Dunning et al. (2010); Rombouts et al. (2010); Dukelsky et al. (2004); Ortiz et al. (2005) , such as superfluidsFoster et al. (2013, 2014), integrable fermion or boson pairing models with nonuniform interactionsJ. Dukelsky (2001); L. Amico (2001), Gaudin magnets (central spin models), and potentially can be extended to a much broader class of integrable nonlinear equations.

The purpose of this paper is therefore twofold. First, it serves as an encyclopedia of quantitatively exact predictions, new and old, for the quench dynamics of real -wave BCS-BEC condensates in two and three spatial dimensions. Readers primarily interested in this aspect of our work will find most of the relevant information in the Introduction, Sect. VII, and Conclusion. In particular, Sect. I.4 concisely summarizes our main results and provides a guide to other sections that contain further results and details. Our second goal is to develop and thoroughly review a method for determining the far from equilibrium dynamics in a certain class of integrable models. We refer readers interested in learning about the method to Sect. II. Also, from this viewpoint, Sects. III and IV should be considered as applications of our approach and Sect. V as a related development.

A major experimental breakthrough with ultracold atoms was achieved in 2004, when they were used to emulate -wave superconductors with an interaction strength that can be varied at will Regal et al. (2004); Zwierlein et al. (2004). Experimental control parameter is the detuning – the binding energy of a two fermion bound state (molecule). This parameter determines the strength of the effective interaction between fermions and can be varied both slowly and abruptly with the help of a Feshbach resonance. Moreover, it is straightforward to make time resolved measurements of the subsequent evolution of the system. Thus cold atoms provide a natural platform to study quenches in superfluids and in a variety of other setupsCalabrese and Cardy (2006); Kollath et al. (2007).

At large we have fermionic atoms with weak effective attraction that form a paired superfluid, an analog of the superconducting state of electrons in a metal. As is decreased, the atoms pair up into bosonic molecules which then Bose condense. It was argued for a long time that both the paired superfluid and the Bose-condensed molecules are in the same phase of the fermionic gas, named the BCS-BEC condensate Leggett (1980); Nozières and Schmitt-Rink (1985). As decreases, the strength of the effective interaction (coupling) between fermions increases from weak to strong and the system undergoes a BCS-BEC crossover. At , where is the Fermi energy, the system is deep in the BCS regime, while at large negative it is deep in the BEC regime. It is not known how to recreate such a crossover in a conventional solid-state superconductor since the interaction strength cannot be easily adjusted.

In a quantum quench setup the system is prepared in the ground state at a detuning . At the detuning is suddenly changed, . At the system evolves with a new Hamiltonian . The main goal is to determine the state of the system at large times, .

i.1 Models and approximations

We consider two closely related models in this paper in both two and three dimensions. The first one is the well-known two channel model that describes two species of fermionic atoms interacting via an -wave Feshbach resonance

It is convenient to think of the two types of fermions of mass and energy as spin up and down, created and annihilated by operators and . The interaction term converts two fermions into a bosonic molecule and vise versa at a rate controlled by the parameter . Molecules are created and annihilated by and and have a binding energy . The parameter is set by the type of atoms and the specifics of a particular Feshbach resonance and cannot be changed in a single experiment; can be varied at will by varying the magnitude of the magnetic field applied during the experiment. This model describes atoms in the BCS regime when is large, which undergo a crossover to the BEC regime as is decreased.

A parameter with dimensions of energy important for our analysis of this model is , where is the bulk density of states (proportional to the total volume) at the Fermi energy . A well known parameter


controls whether the resonance is narrow or broad . This parameter is the dimensionless atom-molecule interaction strength or, equivalently, the resonance width.

A very convenient feature of the narrow resonance is that regardless of the regime of the system, controlled by , the system is adequately described with mean field theory Gurarie and Radzihovsky (2007). This is already clear from the form of the Hamiltonian: small implies that interaction is small.

Broad resonances on the other hand correspond to large . Under those conditions it is possible to integrate out the molecules to arrive at a simpler Hamiltonian Gurarie and Radzihovsky (2007) describing fermions interacting via a short range attractive interaction with variable strength




This is the single (one) channel or BCS model, which is the second model we analyze in this paper. It also describes the BCS-BEC crossover as is decreased ( is increased). However, while in the BCS and (to some extent) in the BEC regimes corresponding to large and small , respectively, mean field theory holds in equilibrium, for the intermediate values of (neither large nor small) the mean field theory is known to break down. A special value of in the middle of the regime unaccessible to the mean field theory already in equilibrium is called the unitary point. It corresponds to the interaction strength where molecules are about to be formed. Non-condensed molecules play an important role in the description of the unitary point and its special properties are a subject of many studies in the literatureNozières and Schmitt-Rink (1985); Burovski et al. (2006).

Just as earlier work on the far from equilibrium superconductivity, we analyze the quench dynamics in the mean field approximation where no molecules are transferred into or out of the BCS-BEC condensate after the quench, i.e. the dynamics of the condensate is decoupled from the non-condensed modes. We analyze the validity of this approximation for nonequilibrium steady states produced by quenches in two channel model in Appendix A. We find that the situation is similar to that in equilibrium Gurarie and Radzihovsky (2007). In the case of a broad Feshbach resonance, mean field is expected to hold for quenches where both initial and final detunings are far from the unitary point. A quench into the unitary point is a very interesting problem addressed by some publications before Bulgac and Yoon (2009), but the method we employ here is not applicable to this case.

Nevertheless, a variety of quenches are still accessible to our description even when the resonance is broad, including BCS BCS, BCS BEC, BEC BCS, and BEC BEC, where BCS and BEC stand for the value of the interaction strength far weaker or far stronger than that at the unitary point. In the case of BCS-BEC superfluids formed with interactions generated by narrow Feshbach resonances, the mean field theory treatment is valid even at the threshold of the formation of the bound state and throughout the BCS-BEC crossover. Here we consider quenches of the detuning for both narrow and broad resonances within mean field. Note that in the case of the one channel model we expect the mean field on the BEC side to be valid only in the far BEC limit where the ground state essentially consists of non-interacting Bose condensed moleculesLevinsen and Gurarie (2006).

In the mean field treatment the condensate is described by the part of the Hamiltonian (I.1), which is decoupled from terms in this approximation. The Hamiltonian therefore becomes




are Anderson pseudospin-1/2 operators Anderson (1958) and

Hamiltonian (6) is also known as inhomogeneous Dicke or Tavis-Cummings model. In a quantum quench problem we need to solve Heisenberg equations of motion for this Hamiltonian for given initial conditions


where , and are Hermitian and anti-Hermitian parts of the operator , and are coordinate unit vectors.

The second step in the mean field treatment of the two-channel model is to replace Heisenberg operator in the first equation of motion in (8) with its time-dependent quantum mechanical average, , which is expected to be exact in thermodynamic limit as long as the state is macroscopically occupied at all times. This replacement can be shown to be exact in equilibrium using the exact solution for the spectrum of inhomogeneous Dicke model Richardson (1977); Gaudin (1983) and numerically for the time dependent problem Sträter et al. (2012). Upon this replacement equations of motion become linear in operators and taking their quantum mechanical average, we obtain


where . These are Hamiltonian equations of motion for a classical Hamiltonian


which describes a set of angular momenta (classical spins or vectors) coupled to a harmonic oscillator. Here, denotes the complex conjugate of . These dynamical variables obey the following Poisson brackets:


where , , and stand for spatial indecies , , and .

Similar steps in the case of the single channel model (3) lead to a classical spin Hamiltonian


together with the corresponding equations of motion.

An important characteristic of the system both in and out of equilibrium is the superfluid order parameter or the gap function defined in the two channel model as


In the one channel limit this expression turns into


The magnitude of the order parameter is known as the Higgs or amplitude mode for its similarity with the Higgs bosonPekker and Varma (2014); Barankov and Levitov (2007) and its time-dependent phase represents a Goldstone mode. Note however that out of equilibrium the gap function does not entirely determine the state of the system. It specifies the effective magnetic field acting on each spin according to Eq. (9), but there is still a certain freedom in how the spin moves in this field. For example, even for a constant field the spin can precess around it making an arbitrary constant angle with its direction.

In above models we took a free single particle spectrum, and labeled states with momenta . This choice is not essential for our analysis. We can as well consider an arbitrary spectrum . The pairing is then between pairs of time reversed statesAnderson (1959); Kurland et al. (2000), see also the first two pages in Ref. Yuzbashyan et al., 2005a for more details. For example, in Hamiltonian (6) this results in relabeling , , etc., where the state is the time reversed counterpart of . Our results below depend only on the density of the single particle states in the continuum limit regardless of whether these states are characterized by momenta or any other set of quantum numbers .

Figure 1: (color online) Ground state chemical potential for two channel model in 3d in units of the Fermi energy as a function of the ground state gap for various resonance width . is calculated from Eqs. (26) and (27). Note that in the two channel model is bounded from above by .

i.2 Ground state

In the ground state


where the magnitude is time-independent. Apart from an overall rotation about the -axis with frequency , the ground state is a static solution of the equations of motion that minimizes . The minimum is achieved when each spin is directed against its effective magnetic field, i.e.




Note that the length of the spin . This is because the ground state is a tensor product of single spin-1/2 wave functions and .

The equation of motion (9) for yields


which implies a self-consistency equation for


Further, the Hamiltonian (10) conserves


which is the average total number of bosons and fermion pairs. This number is related to and the chemical potential as


The Fermi energy is the chemical potential of the fermionic atoms at zero temperature in the absence of any interaction, when only fermions are present. It provides an overall energy scale and it is convenient to measure all energies in units of the Fermi energy. Thus, from now on, we set everywhere below


Below we often switch from discrete to continuum (thermodynamic limit) formulations. In the former version, there are discrete single-particle energy levels with certain degeneracy each. Any quantity we consider in this paper depends on only through , . For example, all spins on a degenerate level are parallel at all times and effectively merge into a single vector. There are such vectors, so we count distinct classical spins.

In thermodynamic limit, energies form a continuum on the positive real axis, i.e. are described by a continuous variable with a density of states that depends on the dimensionality of the problem


where is the bulk density of states (proportional to the system volume) at the Fermi energy, in 2d, and in 3d. Summations over turn into integrations,


With only fermions present, the total particle number is


where is the number of spatial dimensions. Interaction redistributes this number between fermions and bosons as in Eq. (21). Combining Eqs. (21) and (25) and taking the continuum limit, we obtain


where is the dimensionless resonance width defined in Eq. (2).

Similarly, Eq. (19) becomes in the thermodynamic limit


where is the high energy cutoff. In 3d it can be eliminated by an additive renormalization of the detuning , see e.g. Ref. Gurarie and Radzihovsky, 2007. This however does not affect our results for the quench dynamics as they depend on the difference between the initial and final values of the detuning.

Eqs. (26) and (27) contain two independent parameters not counting the cutoff. For example, we can choose and and determine and from these equations, or choose and and determine and etc., see Fig. 1 for a plot of for various in 3d. Note also that is proportional to the number of bosons and is therefore limited by the total number of particles. Eq. (26) implies


i.3 Quench setup and initial conditions

In a quantum quench setup we prepare the system in a ground state at a certain detuning , i.e. the initial state is


where are the ground state values determined by Eqs. (26) and (27) with . We then quench the detuning and evolve the system with the two-channel Hamiltonian (10) starting from the initial state (29) at .

The state of the system is fully determined by the many-body wavefunction, which in the mean field treatment is at all times a product state of the form


where and is the fermionic part of the wave function:


Bogoliubov amplitudes obey the Bogoliubov de Gennes (BdG) equations


with the normalization condition . Apart from an overall time-dependent phase (which is important for certain observables), these equations are equivalent to the classical spin equations of motion (9) and spins are related to the amplitudes as


where is the length of the spin. For quench initial conditions , as explained below Eq. (17).

Each quench is uniquely characterized by three parameters – the resonance width , the initial and final values of the detuning in units of the Fermi energy. Indeed, and determine and and thus the initial condition, while the equations of motion (8) in the thermodynamic limit depend only on and . To see the latter, note that model parameters enter the equation of motion for spin only through , while the equation of motion for the bosonic field can be equivalently written as


Instead of we find it more convenient to characterize the quench by – the ground state gaps corresponding to these values of the detuning. As discussed below Eq. (27), for a given , the detuning uniquely determines and vice versa. Note that has nothing to do with the time-dependent gap function and in particular with the large time asymptote . Whenever goes to a constant at large times, we denote this constant .

i.4 Main results

Our main result is a complete description of the long time dynamics of two and one channel models (10) and (12) in two and three spatial dimensions following a quench of the detuning (coupling in the one channel model) in the thermodynamic limit. A key effect that makes such a description possible is a drastic reduction in the number of effective degrees of freedom as . It turns out that the large time dynamics can be expressed in terms of just a few new collective spins plus the oscillator in the two channel case that are governed by the same Hamiltonians (10) and (12) only with new effective parameters replacing and . The number of collective spins is or 2 and or 2 for one and two channel models, respectively, depending on the quench. The difference is due to the presence of the oscillator degree of freedom in the latter case. For example, means that the effective large time Hamiltonian not only has no spins, but also the oscillator is absent, i.e. . This reduction effect combined with integrability of classical Hamiltonians (10) and (12) allows us to determine the state of the system (its many-body wave function) at . We explain this method in detail in Sect. II. This subsection provides a summary of main results obtained with the help of this method.

In Sects. III and IV, we construct exact quench phase diagrams shown in Figs. 2 - 5. Depending on the values of and either system reaches one of three distinct steady states labeled by I, II (including subregion II’) and III that can be thought about as nonequilibrium phases with second order phase transition lines between them ( limit of the order parameter is continuous along lines separating different regions). These steady states correspond to or 2 collective spins, respectively, for one channel model and to or 2 in the case of two channels.

Each point in the quench phase diagrams represents a particular quench specified by a pair of values . Here is the gap that the system would have in the ground state at detuning , which is a known function of . Values and – ground state gaps for and , respectively – uniquely determine and at fixed resonance width . Note that is not the magnitude of the actual steady state gap function . Each quench (or ) therefore maps to a single point and vise versa.

Steady states I, II, and III reached by the system at can be described in terms of the superfluid order parameter . In region I of phase diagrams in Figs. 2 - 5 the gap function vanishes at large times, , see Fig. 6.

In region II (including subregion II’) the magnitude of the order parameter asymptotes to a nonzero constant as illustrated in Fig. 7,


where are functions of (or, equivalently, of and ), and to be determined below, and is a contstant phase. Plots of and as functions of for fixed are shown in Figs. 9, 18, and 19. The quantity plays the role of the out-of-equilibrium chemical potential. Subregions II and II’ of region II correspond to and , respectively.

In region III of quench phase diagrams the amplitude of the order parameter oscillates persistently at large times as shown in Fig. 8,




dn is the Jacobi elliptic function and is an integration constant. The magnitude of the order parameter oscillates periodically between and . The phase contains linear and periodic partsnot (a)


Constants and are known functions of (or ), and to be specified below, see also Figs. 9 and 10 and refer to Sect. II.4.2 for more information about the periodic solution.

Figure 2: (color online) Detuning quench phase diagrams for two-channel model (I.1) in 2d for an assortment of resonance widths . Each point represents a single quench labeled by (vertical axis) and (horizonal axis) – pairing gaps the system would have in the ground state for initial and final detunings. At large times the system ends up in one of three steady states shown as regions I, II (including II’), and III. For quenches in region I the order parameter vanishes, . In II and in III oscillates persistently. Subregions II and II’ differ in the sign of (out of equilibrium analog of the chemical potential): in II and in II’. The diagonal, , is the no quench line. To the left of it are strong to weak coupling quenches, to the right – weak to strong. in 2d is the maximum possible ground state gap and is the ground state gap corresponding to zero chemical potential, i.e. is given by Eq. (26) for .
Figure 3: (color online) Same as Fig. 2 but in three spatial dimensions.
Figure 4: (color online) Interaction () quench phase diagram for one-channel model (12) in two dimensions. Otherwise same as Fig. 2


Figure 5: (color online) Interaction quench phase diagram for one-channel model (12) in 3d. Otherwise same as Fig. 2. Consider e.g. quenches from fixed infinitesimal coupling (small ) to various final couplings . Increasing () we move through gapless (I), gapped (II), then oscillating (III) steady states. As increases further oscillations disappear and we again end up in a steady state characterized by constant asymptotic (II’).
Figure 6: (color online) in region I for 3d two-channel model, , obtained from numerical evolution of spins following a detuning quench . Here , [cf. Fig. 3(b)]. From these two values all other parameters obtain, e.g. , and .
Figure 7: (color online) in regions II (top) and II’ (bottom) for 3d two-channel model, , obtained from numerical evolution of spins after quenching the detuning . in both panels, same as in Fig. 6. The final detuning corresponds to: (a) and (b) . See also Fig. 3(b).
Figure 8: (color online) Amplitude (Higgs mode) and phase of the order parameter in region III of Fig. 3(c) after detuning quench from deep BCS to BEC in 3d two-channel model for . Numerical evolution with 5024 spins vs. Eqs. (36) and (38). , and .
Figure 9: (color online) Limiting values of for 3d two-channel model at large times after a detuning quench as functions of (or, equivalently, of final detuning ) at fixed small (fixed initial detuning deep in the BCS regime). This corresponds to moving along a horizontal line (not shown) in Fig. 3(a) and (c) going through regions I where , II where , III where oscillates periodically between and , and into region II’ where again . Note that persistent oscillations appear and then disappear again as we decrease (i.e. increase at fixed ). The same behavior is observed in the 3d one-channel model, see Fig. 5.
Figure 10: (color online) Parameter in Eq. (36) for asymptotic in phase III as a function of at fixed small , same as in Fig. 9. For quenches within weak coupling limit , so nonzero quantifies deviations from this limit. Note that one must have , so that the expression under the square root in Eq. (36) is nonnegative.
Figure 11: (color online) Spin distribution as a function of (in units of Fermi energy) at large times after the quench in 3d two channel model. In phases I and II, is the projection of the spin onto its effective magnetic field (-axis in phase I) around which it precesses. In equilibrium (1 in the ground state) for all momenta and in phase I and 1 correspond to doubly occupied and unoccupied states, respectively. Quench parameters are and: (a) (BCS to deep BCS quench in phase I); (b) (BEC to deep BCS quench in phase I). In both cases . Note the Fermi-like shape of the distribution function in (a). Note that as , as it should, indicating that states at very high energies are empty.

Previous studies of the BCS dynamicsVolkov and Kogan (1974); Barankov et al. (2004); Yuzbashyan et al. (2006); Barankov and Levitov (2006); Yuzbashyan and Dzero (2006); Dzero et al. (2007) were performed in the weak coupling regime when both and are much smaller than a characteristic high energy scale (Fermi energy for cold gases and Debye energy for conventional superconductors). This limit corresponds to an infinitesimal vicinity of the origin in our quench phase diagrams in Figs. 2 - 5 . The weak coupling limit is universal in that it is independent of the resonance width and dimensionality and thus is the same in all diagrams. Critical lines separating regions I and II, and II and III are straight lines in this case coming out of the origin with slopes


Further, in Eq. (36) and , take a simpler form given by Eqs. (151) – (153), and (155).

There are several qualitatively new effects beyond the weak coupling regime. At smaller resonance width , gapless region I terminates below at along the vertical axis in 2d. This means that as initial coupling gets stronger ( increases), even quenches to arbitrarily weak final coupling (small ) do not result in vanishing at large times in contrast to the weak coupling regime where quenches with sufficiently large always do. The I-II critical line also displays an interesting backwards bending behavior for , see the inset in Fig. 2(b) and Eqs. (162) and (158).

Region III of persistent oscillations terminates at a threshold value of in 3d, see Figs. 3 and 5. This means that even quenches from an infinitesimally weak initial coupling ( in the one channel model, which corresponds to a vicinity of the normal state) to final couplings stronger then a certain threshold value produce no oscillations and instead goes to a constant. At finite but small initial gap (e.g. along the dashed line in Fig. 3) there is a reentrant behavior in both 2d and 3d as the final coupling () increases when first there are no oscillations, then they appear, and then disappear again. The threshold value of where the critical line separating regions II and III terminates is given by Eq. (171) (plotted as a function of the resonance width in Fig. 22) and Eq. (194) for one and two channel models, respectively. For more details about quench diagrams, such as the shape of the critical lines, various thresholds and termination points, values of parameters (e. g. , and ) characterizing asymptotic , see Sects. III and IV.

The large time asymptote of does not fully specify the steady state. One also needs to know the Bogoliubov amplitudes . We calculate them in Sect. II.4 in all three steady states. In terms of spin vectors, this translates into steady state spin distribution. Even in regions I and II where goes to a constant, the steady state of the system is far from any equilibrium state. Time-independent means that in a frame that rotates around -axis with frequency the magnetic field that acts on spin in Eq. (8) is constant. In equilibrium aligns with or (ground state). In steady states I and II it instead rotates around making a constant angle with it. Let be the angle between and (negative -axis in steady state I), so that in the ground state . Out of equilibrium determines the steady state spin distribution function and is given by Eq. (135). This expression for applies in all three steady states, but its interpretation in region III is slightly different and will be explained below. A plot of the distribution function is shown in Fig. 11. We explore the asymptotic states produced by detuning or interaction quenches in detail in Sect. II.4. In Sect. VII we provide further insight into their physical nature and discuss their experimental signatures.

We perform detailed analysis of linearized equations of motion that goes much beyond previous work even in the weak coupling regime and yields a range of new results. Small quenches of the detuning correspond to a small neighborhood of the diagonal in quench diagrams in Figs. 2 - 5, i.e. they fall within region II where and Eq. (35) applies. We show that within linear approximation and , i.e. there are no corrections to these equations linear in the change of detuning or, equivalently, in . This is in fact a general result that has been overlooked by previous work – to first order in deviations from the ground state always asymptotes to its ground state form for the Hamiltonian with which the system evolves at . Note however that when quadratic correction is taken into account one gets . For example, in the weak coupling regime we find


We obtain an exact expression for , Eqs. (232), (233), and (234), valid at all times and arbitrary coupling strength for both one and two channel models. In weak coupling regime this expression simplifies so that


From here short and long times asymptotes follow. At short times the order parameter amplitude rises or falls sharply as


And the long time behavior in the weak coupling limit is


At stronger coupling in region II (but not II’) the long time asymptote is still given by Eq. (43), only the coefficient in front of the second term on the right hand side is more involved.

Regions II and II’ differ in the sign of the phase frequency , in II and in II’. We will see below that frequency (Fourier) spectrum of quench dynamics in regions II and II’ is , so that the Fourier transform of a dynamical quantity reads . For the phase has a stationary point on the integration path at , while for it is absent. As a result, the long time behavior in three dimensions in region II’ changes


where , is of order one, and . The same expression holds for the one channel model after a replacement . Oscillation frequency and decay are in agreement with Ref. Gurarie (2009) and reflect the fact that in the absence of a stationary point, the long time asymptote is dominated by the end point of integration at , , and the density of states in 3d vanishes as at small .

In two dimensions linear analysis yields a different approach to the asymptote in region II’


because of a constant density of states and divergence of the Fourier amplitude of at small (see below). We also determine the time-dependent phase of the order parameter in all cases corresponding to Eqs. (41) – (45), asymptotes of individual spins as , and many other new results for the linearized dynamics in Sect. V.

Finally, we extend some of the above results for the long time behavior of to the nonlinear regime, though unlike the linear analysis these results are not rigorous. In region II


where is a dimensionless coefficient. This answer holds for both one and two channel models in either dimension.

For region II’ we argue that the answer depends on dimensionality similarly to the linear analysis and


where .

The approach to the gapless steady state (region I) is expected to be


where or , and


We discuss these nonlinear large time asymptotes in more detail in Sect. VI.

Ii Method

Here we describe a method that allows one to determine the asymptotic state of the system at long times. Both the quantum (6) and classical (10) two-channel models are integrable meaning that there are as many nontrivial conservation laws as there are degrees of freedom. There is an exact Bethe Ansatz type solution for the quantum spectrum Gaudin (1983). In the classical case integrability implies a formal inexplicit solution of the equations of motion in terms of certain multivariable special (hyperelliptic) functionsYuzbashyan et al. (2005c) that can be helpful for understanding certain general features of the dynamics. Evaluating specific dynamical quantities of interest for realistic initial conditions with this solution is however roughly equivalent to just solving the equations of motion numerically. But the latter could be as well done directly without the formal exact solution. This is a typical situation in the standard theory of nonlinear integrable systems.

Fortunately, it was realized that at least for the BCS type models the large time dynamics dramatically simplifies in the thermodynamic limit, so that the number of evolving degrees of freedom effectively drops to just a few spins. Building on this insight, Yuzbashyan et. al. Yuzbashyan et al. (2006) were able to develop a method that goes beyond the standard theory and explicitly predicts the long time dynamics in the thermodynamic limit.

The main idea of this method is as follows. First, we construct a special class of reduced solutions of the classical equations of motion for the two-channel model such that the dynamics reduces to that of just few effective spins. Then, we choose a suitable reduced solution and fix its parameters so that its integrals of motion match those for a given quench in the thermodynamic limit. Reduced solutions have only few additional arbitrary constants and cannot generally satisfy all of the quench initial conditions (29). There are initial conditions (two angles per spin plus two initial conditions for the oscillator mode ) and only correspond to the integrals of motion.

Next, exploiting the fact that for fixed BdG equations (32) are linear in the amplitudes and , we derive the most general asymptotic solution that has the same as the reduced one. It has the same integrals as the quench dynamics by construction and, in addition, arbitrary independent constants to match the remaining initial conditions. We conjecture that so constructed asymptotic solution is the true large time asymptote of the actual quench dynamics. To verify this few spin conjecture it is sufficient to show that the large time asymptote of the actual matches that of the reduced (and therefore general asymptotic) solution. We do so numerically in the nonlinear case and analytically for infinitesimal quenches when the dynamics can be linearized.

We consider the two channel model in this and the following sections and then obtain similar results for the one channel (BCS) model in Sect. IV by taking the broad resonance, , limit.

ii.1 Integrability and Lax vector construction

An object called Lax vector plays a key role in our approach. It encodes all the information about the integrals of motion and turns out to be especially useful in analyzing the quench dynamics in the thermodynamic limit. The Lax vector is defined as


where is an auxiliary complex variable and . Poisson brackets of components of satisfy the following Gaudin algebra:


This implies an important equality


Explicit evaluation of yields




It follows from Eq. (53) that these spin Hamiltonians mutually Poisson commute, i.e.


Moreover, the two-channel Hamiltonian (10) is


This implies that and are conserved by and establishes the integrability of the two-channel Hamiltonian. Note that is also conserved for any value of and serves as a generator of the integrals of motion for the two-channel model. The same construction works in the quantum case as well; one only needs to promote classical dynamical variables to corresponding quantum operators and replace Poisson brackets with commutators.

Equations of motion can be conveniently and compactly written in terms of the Lax vector as


Comparing the residues at the poles at both sides of this equation, we see that it is equivalent the equations of motion for spins (9).

The square of the Lax vector is of the form


where is the total number of distinct single particle energies , the product is similarly over nondegenerate values of , and is a polynomial in of degree . The roots of this spectral polynomial (or equivalently of ) play an important role in further analysis of the asymptotic behavior. Note that since is conserved, so are its roots. They thus constitute a set of integrals of motion alternative to Eq. (55). Since for real , its roots come in complex conjugate pairs.

ii.2 Reduced solutions

Let us look for special solutions of equations of motion (58) such that the Lax vector factorizes into time-dependent and independent parts


where (not to be confused with Pauli matrices) denote spins in this solution that can have arbitrary length to distinguish them from spins for the quench dynamics that have length 1/2. Further, are time-independent constants to be determined later and is the Lax vector for an effective -spin system


Here are new collective spin variables placed at new arbitrary “energy levels” . Note that the bosonic field