Dynamical topological phases in quenched spin-orbit coupled degenerate Fermi gas

The spin-orbit coupled degenerate Fermi gas provides a totally new platform to realize topological superfluids and related topological excitations. Previous studies have mainly focused on the properties of the ground state. Here we consider a two-dimensional Fermi gas with Rashba spin-orbit coupling subject to a perpendicular Zeeman field. For this system, we have found that its ground state topological structure is captured by the spin texture, which is readily measurable in experiments. We show that, when the Zeeman field is suddenly quenched, dynamical topological phases can be realized. More specifically, three post-quench dynamical phases can be identified according to the asymptotic behavior of the order parameter. In the undamped phase, a persistent oscillation of the order parameter may support a topological Floquet state with multiple edge states. In the Landau damped phase, the magnitude of the order parameter approaches a constant via a power-law decay, and this system can support a dynamical topological phase with a pair of edge states at the boundary. In the over-damped phase, the order parameter decays to zero exponentially although the condensate fraction remains finite. These predictions can be observed in the strong coupling regime of ultracold Fermi gas.

In the realm of condensed matter physics, one mainly concerns the physical properties of the ground states. This is because the far-from-equilibrium coherent dynamics is generally inaccessible in experiments due to unavoidable relaxation and dissipation from interactions with the environment. Even a sudden change of external parameter generally plays the essential role of only slightly perturbing the ground state, which can be well captured by the perturbation theory. The linear response theory, in this sense, turns out to be extremely successful in understanding most of our observations. However, the dynamics of quantum systems far-from-equilibrium is in general of great interest from a fundamental viewpoint because it can provide us with the properties of the system beyond ground state, for instance, excitations, thermalization, (dynamical) phase transitions and related universalities [1][2][3][4][5][6]; see more details in a recent Review [7]. In this regard, ultracold atom systems provide a totally new avenue for the exploration of intriguing far-from-equilibrium coherent dynamics in both weakly and strongly interacting many-body systems [4,8,9]. This is made possible by the precise control of key parameters in cold atomic systems as well as the ideal isolation from environment [10]. Such controllability has been demonstrated in a number of experiments [11,12].
For these reasons, the coherent dynamics of the s-wave Bardeen-Cooper-Schrieffer (BCS) superfluid has been intensively studied over the past decade [13][14][15][16][17][18][19][20]. In these studies, via the so-called Anderson's pseudospin representation [21], the BCS model can be exactly mapped into a classical spin model, which is proven to be integrable and can be solved exactly using the auxiliary Lax vector method [22,23]. It has been shown that the quench dynamics of the system depends strongly on both the initial and the final values of the quenched parameter, which is often chosen to be the interaction strength. In general, three different phases can be identified according to the long-time asymptotic behavior of the order parameter: the undamped persistent oscillation phase (synchronization phase), the Landau damped oscillation phase, and the overdamped phase. Integrability of the Hamiltonian is essential to understand these results. The p-wave superfluid has the same mathematical structure as the s-wave superfluid, thus similar phases are observed in a recent study [24,25]. As the p-wave superfluid supports topological phases, a quenched p-wave superfluid is found to support dynamical topological phases within certain parameter regimes. However, Fermi gas with p-wave interaction, realized by tuning the system close to a p-wave Feshbach resonance, suffers from strong incoherent losses due to inelastic collisions [26]. We thus do not expect a p-wave superfluid and its coherent dynamics to be realized in the near future of ultracold atom systems .
In this Article, we study the quench dynamics and topological edge states in a spin-orbit (SO) coupled superfluid Fermi gas in two dimension (2D), motivated by the very recent realization of SO coupling in ultracold atoms [27][28][29][30][31][32]. The ground state of this system can be topologically nontrivial in some parameter regimes [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. This is because the SO coupling, Zeeman field and s-wave interaction together can lead to effective p-wave pairing under proper conditions. This system possesses several control parameters that can be readily tuned in experiments. As a result, it naturally provides a practical platform to study the far-from-equilibrium coherent dynamics and related topological phase transitions. However, both SO coupling and Zeeman field breaks the integrability of this model, and change the system from a single-band to a two-band structure. It is therefore natural to ask the fundamental question: "What types of post-quench dynamical phases this system will exhibit, and how do these dynamical phases differ from the ones supported by the integrable models?" In our study, we choose the Zeeman field as the quench parameter, which has been realized in recent experiments [27][28][29][30][31][32]. It is surprising that all the phases supported by the integrable model still exist in the non-integrable SO coupled superfluid Fermi gas. We provide a complete phase diagram and investigate each phases, including their topological properties, in great detail. We also show how dynamical topological phases, which can support topologically protected edge states, emerge in this model.

Model Hamiltonian
We consider a 2D system of uniform SO coupled degenerate Fermi gas with s-wave interaction confined in the xy-plane, whose Hamiltonian is written as where s, s =↑, ↓ label the pseudospins represented by two atomic hyperfine states, k is the momentum operator, ξ k = k − µ, where k = k 2 /2m denotes kinetic energy and µ is the chemical potential, α is the Rashba SO coupling strength, σ x,y,z are the Pauli matrices, h is the Zeeman field strength along the z-axis, and c ks is the annihilation operator for fermions. g represents the inter-species s-wave interaction strength. Notice that the 2D system is created from the three dimensional system by applying a strong confinement along the z-axis, thus g in principle can be controlled by both confinement and Feshbach resonance. This 2D degenerate Fermi gas has been realized in recent experiments [52][53][54][55]. We consider the quench dynamics of the Fermi gas at the mean field level, thus the potential defect production, i.e., the so-called Kibble-Zurek mechanism [56,57], after the quench is not considered. Imagine that we prepare the initial system in the ground state with Zeeman field h i . At time t = 0 + , we suddenly change the Zeeman field strength from h i to some final value h f . This scheme should be in stark contrast to previous literatures, in which the interaction strength g generally serves as the quench parameter [13-20, 24, 25]. We choose the Zeeman field strength as our quench parameter for the following reasons. (1) As we shall show in the discussion of the ground state properties of this system, the Zeeman field strength directly determines the topological structure of the ground state. (2) In SO coupled quantum gases, the laser intensity and/or detuning serve as effective Zeeman field, and these parameters can be changed in a very short time scale, satisfying the criterion for a sudden quench. By contrast, a change of the interaction strength is achieved by tuning of the magnetic field which usually cannot be done for too short of a time scale. (3) Moreover, quenching the effective Zeeman field in SO coupled quantum gases has already been demonstrated in recent experiments [27][28][29][30][31], in which both the magnitude and sign of Zeeman field can be changed.
The superfluid order parameter is defined as ∆ = gΣ k c −k↓ c k↑ and the interaction term can be rewritten as where the Bogoliubov-de Gennes (BdG) operator reads as [33,34] Assuming f k± = (u k± , v k± , p k± , q k± ) T are the two energy levels of M k with positive eigenvalues, the order parameter can be determined by solving the gap equation ∆ = −gΣ k (v k+ q * k+ + v k− q * k− ) self-consistently. Also, g should be regularized by 1/g = − k 1/(2 k + E b ) [34]. As a result, in the following, the interaction strength is quantitatively defined by the binding energy E b ∈ [0, ∞).
The coherent dynamics of this model cannot be solved by mapping it into the classical spin model, thus it should be computed with exact numerical method. As a reasonable but general assumption, we assume that the wavefunction, after quench, is still BCS-like, i.e., |ψ(t) = k,s=± Γ k± (t)|0 , where Γ k± (t) = f † k± (t)Ψ k , and the vector f k± (t) are determined by the following time-dependent BdG equation, Here M k is the time-dependent BdG Hamiltonian in which the order parameter ∆(t) now evolves in time after the quench.
Ground State Properties. Before we turn to the discussion of the quench dynamics, let us first briefly outline the ground state properties of the system. The most salient feature of the system is a topological phase transition induced by the Zeeman field [33][34][35][36][41][42][43][44][45][46][47][48][49][50][51]. Without the loss of generality, we  assume that h ≥ 0, in which case, the spin-up atom (spindown atom) represents the minority (majority) species. It is well known that the topology of the superfluid is encoded in the topological index W , which corresponds to the topological state with W = 1 if h > µ 2 + ∆ 2 , while h < µ 2 + ∆ 2 yields W = 0 and the state is non-topological. In Fig. 1(a), we show how the order parameter and the single-particle excitation gap changes as a function of the Zeeman field strength h. At a critical point, h c = µ 2 + ∆ 2 , the single-particle excitation gap vanishes, which represents a topological phase transition. This feature is also essential for the realization of the far-from-equilibrium coherent evolution; see Discussion. To see this phase transition more clearly, we also examine the spin population, as shown in Fig. 1(b). We find that S z (0) , the spin component along the z-axis (which is just the population difference between the two spin species at k = 0), jumps discontinuously when the Zeeman field crosses the critical value h c , while the total spin polarization which is defined as S p = (n ↑ − n ↓ )/n, changes smoothly with respect to h. In fact, as shown below, the jump of S z (0) just implies a change of the topology of the spin texture.
These two different topological regimes feature very different ground state spin textures as shown in Fig. 2. Here we define the atomic spin vector as S k↓ c k↓ , whose ground state expectation values are given by S x (k) = 2 η=± [p * kη q kη ], S y (k) = 2 η=± [p * kη q kη ], and S z (k) = η=± (|q kη | 2 − |p kη | 2 ), from which the Skyrmion number is defined as where s(k) = S(k) /| S(k) | is the normalized spin vector, which maps the 2D momentum space to a unit sphere S 2 . In this sense, Q is nothing but the number of times the spin vector wraps around the south hemisphere. Note that for momenta with fixed magnitude |k|, s x and s y always sweep out a circle parallel to the equator. In the topologically trivial regime, we always have s z (0) = 0, see Fig. 1(b). Hence as |k| increases from zero, s begins at the equator, descends toward the south pole and then returns to the equator as |k| → ∞. Thus in this regime, s(k) initially sweeps out the shaded region in the southern hemisphere of Fig. 2(c), but then unsweeps the same area, resulting in a vanishing winding number Q = 0. In contrast, in the topologically non-trivial regime, we have s z (0) = −1, see Fig. 1(b). Hence as |k| increases from zero to infinity, the spin vector covers the entire southern hemisphere exactly once, as shown schematically in Fig. 2(d), which leads to a non-trivial winding number Q = −1. The sudden change of spin polarization is due to band inversion transition across the critical point.
It is worth pointing out that, for any quench, we always have ∂ ∂t S z (0, t) = 0 (see Methods), regardless of the initial and final values of h. Therefore, Q is unchanged over time after the quench. We can see that, the Q and W are equivalent in equilibrium, but their dynamics after a sudden quench are different. As shown below, W , which describes the topology of the full spectrum of the Hamiltonian, will evolve in time, while Q, which reflects the topological nature of the the state itself, will not. A similar conclusion is found in the study of the quenched p-wave superfluid [24,25]. We emphasize that the momentum space spin texture studied here can be measured in cold atom experiments using the standard time-offlight technique [30,58,59].
Phase diagram of the quenched spin-orbit coupled superfluid condensate. The different phases in this figure is obtained by the long-time asymptotic behavior of the order parameter upon the quench of the Zeeman field from initial value hi to final value h f . The diagonal light blue line, with hi = h f , is the case without quench, thus the quantum state is unchanged. hc marks the quantum critical point separating the topological superfluid and nontopological superfluid in the equilibrium ground state, which is determined by h 2 c = ∆ 2 + µ 2 . Three different dynamical phases observed in this system are labeled with I, II, and III by green, white and purple shaded areas, respectively. In phase I, |∆(t)| shows persistent oscillations, which is from the collisionless coherent dynamics. The dark blue dashed line separates phase I into non-topological Floquet state denoted as NTSFloquet and topological Floquet state labeled as TSFloquet. In phase II, |∆(t)| → ∆∞, a nonzero constant value, which serves as the basic parameter to determine the longtime asymptotic behavior. The orange dashed lines are the nonequilibrium extension of the topological phase transition at h = hc, which separates phase II into two parts, NTS and TS accordingly. Inside NTS (TS) region, the quasi-stationary steady state is trivial (nontrivial) phase without (with) topologically protected edge modes. W = 0 or 1 marks the topological index at t = +∞. In phase III, |∆(t)| → 0 due to strong dephasing from the out-of-phase collisions . All other parameters are identical to Fig. 1.
Dynamical phase diagram. We now turn to our discussion on the quench dynamics. As in Ref. [24,25], we capture the dynamics using a phase diagram presented in Fig. 3. The phase diagram contains three different phases, identified by the distinct long-time asymptotic behaviors of the order parameter in the parameter space spanned by the initial and final values of the Zeeman field h i and h f . The three phases are labeled as phase I, II and III in Fig. 3. More specifically, in the undamped oscillation phase (phase I), the magnitude of the order parameter oscillates periodically without damping, although the wavefunction does not recover itself periodically. In this regime, the spin polarization decays very fast to an equilibrium value via interband Rabi oscillation. In the damped oscillation (phase II) regime, the order parameter exhibits damped oscillation with a power-law decay. In the overdamped phase (phase III), the order parameter decays to zero exponentially. We also show that within phase I and II, there exists dynamical topological regimes where topological edge states emerge in the asymptotic limit. Next, we shall discuss the properties of each phase in more details.
Phase I. In Fig. 4 we plot the dynamics of the magnitude of the order parameter for a typical point in phase I (point A in Fig. 3), from which we find that |∆(t)| oscillates asymptotically as [18] where dn[u, k] is the periodic Jacobi elliptic function, and ∆ + and ∆ − are the maximum and minimum value of |∆(t)|, respectively. This multi-soliton solution is first derived by Barankov et al. for the BCS superfluid without SOC and Zeeman field in the BCS limit [18]. The fitted result with κ = 0.9772 using this empirical formula are presented in Fig. 4, which perfectly agrees with the numerical results. Note that similar persistent oscillation has also been observed in conventional BCS sand p-wave superfluids, which is interpreted as a synchronization effect, in which each pseudospin rotates in an effective magnetic field defined by the order parameter and kinetic energy. We do not have such a simple picture here, although non-harmonic oscillation between ∆ + and ∆ − can still be observed. This basic observation indicates that the undamped oscillation of the superfluids should be a quite general feature, not depending on whether the model is integrable or not.
To get more insights into this persistent oscillating behavior, it is helpful to investigate the spin population dynamics. We found that, as shown in the inset of Fig. 4, after the quench S p oscillates but the the oscillation amplitude decays very fast. As a result, S p reaches an equilibrium value very quickly, although perfect periodic oscillation of the order parameter persists. Recently, the dynamics of the spin polarization after quench in a Fermi gas above the critical temperature has been measured in experiments [29]. The decay of the spin polarization can be attributed to the interband Rabi oscillation, in which, different momentum state has slightly different Rabi frequency, such that destructive interference gives rise to the damping phenomenon. This result indicates that the persistent oscillation of the order parameter are not accompanied by a similar oscillation in the spin population dynamics. We find that the phase of the order parameter also changes dramatically in this case, and the asymptotic order parameter can be represented as ∆(t) = |∆(t)|e −i2µ∞t+iϕ(t) [24,25], where the phase ϕ(t) (modulo 2π) is also a periodic function with commensurate period to |∆(t)|. The phase factor's linear piece in time, −2µ ∞ t, can be gauged out by a unitary transformation (see Methods) [25]. After the gauge transformation, we obtain a BdG HamiltonianM k (t) that is periodic in time. For such a periodic system, we may invoke the Floquet theorem to examine its Floquet spectrum [67,68] To determine whether the system is topological in the Floquet sense, we calculate the spectrum ε k± in a strip by adding a hardwall boundary condition in the x-direction. Two examples of the spectrum are plotted in Fig. 5. In the example shown in Fig. 5(a), the spectrum is gapped, corresponding to a topologically trivial Floquet state. On the other hand, the spectrum shown in Fig. 5(b), exhibits gapless edge modes and hence can be regarded as a topological Floquet state. In the phase diagram of Fig. 3, inside phase I, the two dark blue dashed lines characterize the topological boundaries, which separate the non-topological states denoted by NTS Floquet from the topological states denoted by TS Floquet .
Phase II. In this Landau damped phase, the magnitude of the order parameter undergoes damped oscillation and finally reaches finite equilibrium value. Two examples (corresponding to points B and C in Fig. 3) are shown in Fig. 6. Here the magnitude of the order parameter can be described by the following power-law decay function [13][14][15][16] where ∆ ∞ is the magnitude of the order parameter in the longtime limit and by no means it is equal to the order parameter ∆ f determined by a Zeeman field of strength h f in equilibrium. E ∞ is the minimal band gap of the effective Hamiltonian at t → ∞. It should be pointed out that, unlike the conventional BCS model, E ∞ does not necessarily equal to ∆ ∞ in the current model, because of the SO coupling and Zeeman field. One can see that the second term gives rise to the decay of the order parameter. The exponent α, characterizing the power-law decay, is not a universal constant in this model. This is in distinct contrast to the BCS model without SO coupling and Zeeman field, where α = 1/2 [13,15,16] in the BCS limit, and α = 3/2 [14] in the Bose-Einstein Condensation (BEC) limit.
In this phase, the order parameter behaves as ∆(t) = ∆ ∞ e −i2µ∞t in the asymptotic limit [24,25]. Again we can gauge out the phase factor linear in time and treat µ ∞ as an effective chemical potential [25]. We can therefore construct a time-independent BdG Hamiltonian by replacing the chemical potential and order parameter in Eq. (2) with µ ∞ and ∆ ∞ , respectively. This is still a dynamical phase because µ ∞ = µ f , and ∆ ∞ = ∆ f , where ∆ f and µ f are equilibrium order parameter and chemical potential with Zeeman field h f . For example, for point B, we have ∆ f = 0.662E F and µ f = 0.199E F , whereas numerically we obtain ∆ ∞ ≈ 0.456E F and µ ∞ ≈ −0.019E F . Given the asymptotic timeindependent BdG Hamiltonian, the region of the dynamical topological phase can be determined by the condition  with a topological index W = 1, otherwise we have a nontopological dynamical phase with W = 0. In the following we will show that dynamical edge state can indeed be observed in the topological regime. We have also calculated the Chern number C of the system, which shows that C = 1 (0) in the topological (non-topological) regime defined above. Thus W is sufficient to characterize all the dynamical phases in this model. These results pave a new way to realize topological phase in this model.
The topological nature of the system can also be manifested by examining the existence of the edge states. To this end, we obtain the BdG spectrum by adding a hardwall boundary along the x-direction. Examples are shown in Fig. 7. We show the energy gap at zero momentum, , as a function of the final Zeeman field h f in Fig. 7(a) for fixed initial Zeeman field. The closing and reopening of energy gap E g signals the dynamical topological phase transition. Indeed, we show that the dynamical edge state can be observed in the topological regime, see Fig. 7(d) where bulk spectrum is gapped and edge state remains gapless after Zeeman field exceeds the critical value.
Phase III. In phase III, the order parameter quickly decays to almost zero (see Fig. 8) according to ∆(t) ∼ exp(−t/T * ), where T * ∼ 1/∆, the decay time, is equal to the order parameter dynamical time (see Discussion) . However, we need to emphasize that zero order parameter in the long-time limit does not mean the system has become a normal gas. We demonstrate this by showing the dynamics of both singlet and triplet condensate fraction, defined by n s = k | c k↑ c −k↓ | 2 /n and n t = k | c k↑ c −k↑ | 2 /n, where n is the total density. One can see that in the longtime limit, the condensate fraction remains non-zero even though the order parameter vanishes. Non-zero condensate fraction means that the system still contains non-trivial pair correlation. However, the pairing field for different momenta oscillates at different frequencies, which leads to dephasing and hence a vanishing order parameter.

Discussion
In this Article, we have demonstrated that dynamical topological phase can be realized in an SO coupled degenerate Fermi gas by quenching Zeeman field. The Zeeman field directly determines the topological properties of the ground state, which is completely characterized by the spin texture in z-direction at zero momentum S z (0). We want to emphasize that this quantity is directly measurable in cold atom experiments using the standard time-of-flight technique. We have further mapped out the post-quench phase diagram according to the asymptotic behavior of the order parameter. In the undamped phase, the persistent oscillation of the order parameter may support topological Floquet state with multiple edge states. In the Landau damped phase, the magnitude of the order parameter gradually approaches a constant via a powerlaw decay, and this phase contains a dynamical topological portion in certain parameter regions. One pair of edge modes can be observed in this case. In the over-damped phase, the order parameter quickly decays to zero exponentially while the condensate fraction remains finite.
The presence of the SO coupling and the Zeeman field breaks the integrability of our model. However, the same types of post-quench dynamical phases observed in our model are also present in integrable models studied previously. This raises the important question on the relationship between integrability and the long-time asymptotic post-quench behavior of the superfluid/superconducting system. This issue has been intensive studied in some other models regarding relaxation, thermalization and phase transitions [2,7,60,61], in which integrability plays the most essential role therein. Our work here shows that this is a rather subtle question and further studies are needed to provide a definitive answer.
We finally comment on the feasibility of observing the exotic dynamical topological phases unveiled in this Article. The dynamics of the superfluids are mainly determined by two characteristic time scales, that is, the energy relaxation time τ ∼ E F /E 2 gap , where E gap is the energy gap of the superfluids before quench, and the order parameter dynamical time τ ∆ ∼ 1/∆. The quench process from h i to h f can be realized via a frequency jump of the lasers at the time scale of several microseconds [31], which is much smaller than 1/E F . The far-from-equilibrium coherent evolution can be realized when τ > t ≥ τ ∆ [15,18]. The ultracold Fermi gas provides a natural system to explore the physics in the far-from-equilibrium condition at a time scale of 1/E F . In the BCS limit ∆ = √ 2E F E b approaches zero. Thus E gap = ∆ E F (µ > 0 in the BCS limit), and we immediately have τ τ ∆ (Using point A in Fig. 3 as an example, we have ∆ ∼ 0.013E F , µ ∼ 0.7E F , E gap ∼ 0.009E F , thus τ ∆ ∼ 70/E F and τ ∼ 10 4 /E F ∼ 160τ ∆ ). In this new system, the SO coupling and the Zeeman field can greatly change the band structure of the superfluids. For example, the energy gap is no longer determined solely by the order parameter and chemical potential, but instead, it is a very complex function of all parameters; see Methods. At the boundary of topological phase transition, we have E gap = E g = 0. In the vicinity of this boundary, τ ∼ E F /(h − ∆ 2 + µ 2 ) 2 , and we naturally expect that τ τ ∆ . We should emphasize that this condition, originally can only be realized in the BCS limit, now can be realized very easily in the strong coupling regime. Meanwhile, the temperature effect is also a critical issue in ultracold atomic system. In the BCS limit, we expect , where γ 0.577 is the Euler's constant. The required temperature should be very low in order to observe coherent dynamics of superfluids , which is of great challenge to the current experiments [52][53][54][55]. This dilemma can be greatly resolved in our model because of the dramatic change in band structure caused by SO couping and Zeeman field. In the strong coupling regime, we expect the critical temperature to be determined by the Kosterlitz-Thouless transition temperature T KT ∼ 0.1E F [34,64], which obviously can be experimentally accessed [52][53][54][55]. For these reasons, we expect that the relevant dynamics of order parameter and associated dynamical topological phase transitions in phase I and phase II regime can be realized using realistic cold atom setup at the currently achievable temperatures.

Methods
Equation of motion. The basic Hamiltonian in Eq. 1 is solved using standard mean field method by defining the or-der parameter ∆ = −g k c k↑ c −k↓ . Then the dynamics of the model is governed by the effective Hamiltonian M k in Eq. 2. The initial value of chemical potential and order parameter are determined by minimizing the thermal dynamical potential [33][34][35][36][41][42][43][44][45][46][47][48][49][50][51], which is equivalent to solve the coupled gap equation and number equation where E λ k± = λ ξ 2 k + α 2 k 2 + h 2 + |∆| 2 ± 2E 0 are the quasi-particle excitation energy, and λ = ± correspond to the particle and hole branches respectively, E 0 = h 2 (ξ 2 k + |∆| 2 ) + α 2 k 2 ξ 2 k . Throughout our numerical calculations, the energy unit is chosen as the Fermi energy 2πn is the Fermi vector for a noninteracting Fermi gas without SOC and Zeeman field in 2D. We only consider the physics at zero temperature. Throughout this work, we choose the binding energy E b = 0.2E F , and the corresponding scattering length k F a 2D = 2E F /E b 3.1623 and ln(k F a 2D ) 1.1513, which is around the BEC-BCS crossover regime. Since k F a 2D > 1, the superfluid is still made by weakly bound Cooper pairs. For a typical Fermi gas of 6 Li, k F ∼ 1/µm, E F ∼ 1kHz , and the basic time scale discussed in this work is 1/E F ∼ 1ms . The long time collective oscillations of degenerate Fermi gas that much longer than this basic time scale has been demonstrated in experiments [65,66].
The initial value of µ and ∆ can be directly determined by solving the gap equation and number equation selfconsistently, with which the initial wavefunction at t = 0 − is constructed. In this case, the topological phase can be realized when [33,34], However, this phase can only be realized in the BEC-BCS crossover regime. A simple argument is that, in both the BEC and BCS limit, µ 2 ∆ 2 , thus in these two limits, |h| |∆|, and the pairing is destroyed by Pauli depairing effect. As a result, the topological phase can only be realized in a small parameter window near strong coupling. In the following, we choose Zeeman field as the quench parameter instead of others, because it is the most easily controllable parameter in the current experiments [27][28][29][30][31][32]. Moreover, the Zeeman field is directly relevant to the topological boundary while many-body interaction is not; see Eq. 7 and Eq. 10.
Immediately after the sudden change of the Zeeman field, the system's wavefunction is assumed to keep the following BCS form where Γ k,± = f † k± Ψ k , with f k± satisfy the time-dependent BdG equation in Eq. 3. We can verify straightforwardly that Eq. 3 is equivalent to i∂ t |Ψ(t) = H|Ψ(t) .
The above semiclassical equation can be derived from δ Ψ * L = 0, where L = Ψ|i∂ t − H|Ψ . Spin texture at k = 0. The Hamiltonian at k = 0 can be reduced to We first consider the spin texture in the stationary condition. If |h| < µ 2 + ∆ 2 , the two eigenvectors with positive eigenvalues are f 0,+ = ( Then, from the expression of S z (k) from the main text, we can directly get that S z (0) = 0 and Q = 0, which means the system has topologically trivial spin texture. In contrast, if |h| > µ 2 + ∆ 2 , the two eigenvectors will become f 0, 2 , 0) T respectively. So, in this region, we will have S z (0) = −1 and Q = −1, which means the system has a topologically non-trivial spin texture. We see that the sudden changes of S z (0) is due to the band inversion transition across the critical point.
The time-evolution of Q can be discussed in a similar way. When k = 0, we have i∂ t q ± = ∆ * v ± − hq ± and i∂ t p ± = ∆ * u ± + hq ± . Substituting these equations intȯ S z (t) = η=±q * η q η + q * ηqη −ṗ * η p η + p * ηṗη , we immediately find thatṠ z (t) = 0. Thus we have the important conclusion that S z (0, t) = S z (0, 0) , which means Q remains unchanged over time. We need to emphasize that this spin texture is totally different from the pseudospin texture discussed in Ref. 24. The true spin texture discussed in this work can be directly probed in experiments from the time-of-flight imaging [30,58,59], and Q(t) = Q(0) can be directly verified.
Dynamical edge state in phase II regime. In this regime, the magnitude of the order parameter will gradually approach a constant, while its phase factor oscillates periodically, i.e., ∆(t) → ∆ ∞ e −iµ∞t up to a trivial constant, which is the only time-dependent parameter in the BdG Eq. 3. This oscillating phase can be gauged out by defining f k± = (ũ k± e −iµ∞t ,ṽ k± e −iµ∞t ,p k± e +iµ∞t ,q k± e +iµ∞t ) T e −iE k,± t , wheref k± = (ũ k± ,ṽ k± ,p k± ,q k± ) T . Inserting this wavefunction to Eq. 3, we find that whereM k is identical to Eq. 2 except that µ = µ ∞ and ∆ = ∆ ∞ . We immediately see that µ ∞ is the effective chemical potential of the model in the quasi-equilibrium condition. Note that this phase is still dynamical phase because the µ ∞ = µ f and ∆ ∞ = ∆ f , where µ f and ∆ f are equilibrium chemical potential and order parameter with Zeeman field h f ; See our numerical results in the main text. This model can support dynamical edge state in the topological regime defined by Eq. 7. Similar to the analysis in Ref. 33, we can prove exactly that the bulk system is always fully gapped except at the critical point of h c for k = 0. Thus the closing and reopening of the gap provides important indications for topological phase transition. To see the topological phase transition more clearly, we consider a strip superfluids with length L by imposing hard wall boundary at the x direction. To this end, we replace k x → −i∂ x , while k y remains as a good quantum number. Along the x direction, we construct the wavefunction using plane wave basis [62,63] where N max is the basis cutoff. Upon inserting this ansatz into Eq. 14, we can convert the matrixM k into a 4N max by 4N max matrix, whose diagonalization directly leads to the protected modes of dynamical edge state dictated by non-trivial topological invariants. Empirically, we found N max = 200 is a good basis cutoff for a long strip with L = 200k −1 F . The numerical results are presented in Fig. 7.
In the long-time limit, the order parameter in phase I approaches ∆(t) = |∆ ∞ (t)|e −2iµ∞t+iϕ(t) , where |∆ ∞ (t)| is periodic in time, e.g. see Fig. 4. We make a gauge transformation, similar to that in Eq. 14, by identifying µ ∞ as the effective chemical potential, and we obtainM k (t) from Eq. 2 by replacing µ with µ ∞ and ∆ with |∆(t)|e iϕ(t) . Obviously, M k (t) =M k (t + T ), where T is the period determined by both |∆(t)| and ϕ(t). Now we assume the eigenvectors of the above effective Hamiltonian to bef k± (t) = Φ k± (t)e −iεk±t , where Φ k± (t + T ) = Φ k± (t). Then we have where ε k± is the quasiparticle spectra. Similar to the discussion of dynamical edge state in the previous section, we impose a hard wall boundary condition along x direction with length L. We expand the wavefunction in the following way, (18) where N max and M max are basis cutoff for the spatial and temporal expansion, and A = (2/LT ) 1/2 is the normalization constant. Then Eq. 17 can be recasted into a sparse selfadjoint complex matrix form of size (4N max ×(2M max +1))× (4N max ×(2M max +1)). The direct diagonalization of this matrix gives rise to Floquet spectrum. In practice, since we are only concerned with the eigenenergies close to zero, we could utilize the shift and invert spectral transformation and compute only a portion of eigenenergies using the ARPACK library routines. For instance, we choose cutoffs as N max = 200 and M max = 15 and only compute 500 eigenvalues around zero-energy out of total 24800 ones for a given k y . The results are presented in Fig. 5. The robustness of these protected edge states are also examined by slightly changing the model parameters, in which we find that the linear dispersions of these edge states are unchanged.