Topological phase transition in the quench dynamics of a one-dimensional Fermi gas

We study the quench dynamics of a one-dimensional ultracold Fermi gas in an optical lattice potential with synthetic spin-orbit coupling. At equilibrium, the ground state of the system can undergo a topological phase transition and become a topological superfluid with Majorana edge states. As the interaction is quenched near the topological phase boundary, we identify an interesting dynamical phase transition of the quenched state in the long-time limit, characterized by an abrupt change of the pairing gap at a critical quenched interaction strength. We further demonstrate the topological nature of this dynamical phase transition from edge-state analysis of the quenched states. Our findings provide interesting clues for the understanding of topological phase transitions in dynamical processes, and can be useful for the dynamical detection of Majorana edge states in corresponding systems.

Introduction.-Topological phases and phase transitions [1] in quantum many-body systems have recently attracted much attention in various physical contexts. A particularly interesting topological state is the topological superfluid (TSF) phase, where Majorana edge states may exist at its boundary [2,3]. Besides systems with intrinsic chiral p-wave pairing superfluidity [4][5][6][7][8], TSF phases can also be induced from an s-wave pairing superfluid in the presence of spin-orbit coupling (SOC) and an effective Zeeman field [9,10]. While the search for Majorana zero modes is being actively pursued in condensed matter systems such as the semiconductor/superconductor heterostructures [11,12], the possibility of realizing topological superfluidity and hence Majorana zero modes in ultracold atomic gases has also been extensively studied, due to the recently implemented synthetic spin-orbit coupling in these systems [13][14][15][16][17].
An advantage of ultracold atomic gases is the controllability, which not only allows for experimental characterization of the system over a wide range of parameters, but also provides a convenient way to study the dynamical processes. In recent years, dynamical processes in ultracold atomic gases play an increasingly important role in revealing key properties of the system [18]. Particularly, quench dynamics, in which the system evolves from a pure state following a sudden change of parameters, has been widely used in cold atoms experiments [19][20][21][22][23][24]. In the long-time limit, the quenched state is typically quite different from an equilibrium state described by the ensemble theory [25,26]. More importantly, it has been shown that depending on the quench parameters, the quenched state in this limit can exhibit different dynamical phases [27,28].
With the prospect of realizing topological orders in ultracold atomic gases, it is timely to study dynami-cal processes in systems with topological order, and ask questions like: whether topological properties are robust against dynamical perturbations; how topology should be defined and probed in dynamical processes, etc. On the one hand, in an equilibrium state with topological order, features like edge states are protected by a bulk gap [29], while such protection may not be effective in dynamical processes, especially in quench dynamics, where the sudden change inevitably induces high-energy excitations [18]. On the other hand, despite many recent studies on topological phases and edge states in quantum quenches, topology in a dynamical process is still an on-going research that needs further characterization [27,28,[30][31][32][33][34].
In this work, we theoretically study the quench dynamics of a one-dimensional ultracold Fermi gas in an optical lattice potential under the synthetic spin-orbit coupling that has recently been realized experimentally [13][14][15]. Starting from a ground state of the system, we implement a sudden change of the interaction strength and focus on the dynamics of the pairing gap as the system evolves under the quenched Hamiltonian. At sufficiently long evolution times, the pairing gap either approaches a finite value or oscillates in time, depending on the quench parameters. Interestingly, we identify a critical interaction strength in the quenched Hamiltonian, at which the pairing gap in the long-time limit undergoes an abrupt jump across the critical value that corresponds to the emergence of the topological phase transition in the ground state of the pre-quench Hamiltonian. We associate this critical interaction strength with a dynamical phase transition. We further demonstrate that edge states exist in the long-time limit only for processes with quenched interaction strength below the critical value of the dynamical phase transition, which suggests that the dynamical phase transition is also of topological nature. The topological property of the quenched state is thus determined by the quenched Hamiltonian, regardless of the initial state. Our findings provide clues for the understanding of topological phase transitions in dynamical processes, and may have interesting implications for future experiments.
Model Hamiltonian and phases at equilibrium.-Under the single-band approximation, the tight-binding Hamiltonian can be written as: whereĉ † jσ is the creation operator of fermions on site j with spin σ =↑, ↓, and the on-site interaction is attractive with U > 0. The effective Zeeman field h and the spinorbit coupling strength α are related to parameters of the Raman process generating the synthetic spin-orbit coupling.
Following the stand BCS-type mean-field theory, we can Fourier-transform Hamiltonian (1) into momentum space, and write the effective mean-field Hamiltonian in the grand canonical ensemble: whereĉ † kσ is the creation operator for fermions of momentum k and spin σ, with k defined on the first Bril- louine zone. The dispersion of the lowest band ǫ k,↑/↓ = −2 cos k − µ ± h, the spin-orbit coupling α k = 2iα sin k, and the pairing gap ∆ = U k ĉ −k↓ĉk↑ . The chemical potential µ is related to the on-site particle density n, which is kept fixed throughout the lattice. Here, we have 0 < n < 1 for µ < 0, and 1 < n < 2 for µ > 0. Due to the particle-hole symmetry of the effective Hamiltonian, we will focus on the case of µ > 0 henceforth. At equilibrium, the ground state of the system can be determined by minimizing the thermodynamic potential, which at zero temperature is reduced to Ω = Ĥ ef f . Here · · · is taken with respect to the ground state. Previous studies show that the ground state ofĤ ef f is a TSF with Majorana edge states for (µ − 2g s ) 2 + ∆ 2 ≤ h ≤ (µ + 2g s ) 2 + ∆ 2 , and is a conventional superfluid otherwise [35]. In Fig. 1, we map out the ground state phase diagram at zero temperature, which serves as a basis for further studies of quench dynamics. Consistent with previous works, the system is in the TSF phase in the weakly-interacting regime, and undergoes a topological phase transition to become a conventional superfluid (SF) phase in the strongly-interacting regime with the phase boundary labeled by U c . With our choice of parameters (h = 2g s and α = g s ), the transition boundary in Fig. refphasediagram is of first-order. This is manifested not only in the discontinuity of the pairing gap ∆ at the transition boundary, but also in the shape of the thermodynamic potential Ω(∆) (see right panel of Fig. 1). For the convenience of later discussions, we label the critical pairing gap of the topological phase transition Quench dynamics.
-We now study the quench dynam- ics of the system where the interaction U undergoes a sudden change. We assume that the system is initially in the zero-temperature ground state with U = U i . At the start of the quench (t = 0), the interaction is suddenly switched to U = U f , driving the system out of equilibrium. The system then evolves under the Hamiltonian with the quenched interaction U f . To characterize the dynamics of the system, we introduce operators: Together with the momentum-space density operatorn kσ =ĉ † kσĉ kσ , we may express the timedependent mean-field effective Hamiltonian as: whereτ z k =n k↑ +n −k↓ − 1,σ z k =n k↑ −n −k↓ and ǫ k = −2 cos k − µ. From the Schrödinger's equations, one may derive a set of closed equations for the expectation values of these operators. The dynamical parameters of the system can be calculated numerically from these equations [36]. Note that similar approach has been applied to study the quench dynamics in a Fermi gas throughout the BCS-BEC crossover without SOC [37][38][39][40], and more recently, the quench dynamics of a chiral p-wave superfluid near the topological phase transition [27]. With this approach, we focus on the dynamics of the pairing gap ∆(t) = U f k ĉ −k↓ĉk↑ during the evolution. Interestingly, the dynamics and the quenched states are quite different depending on the topological nature of the initial state.
The right panel of Fig. 2 shows typical quench dynamics of the pairing gap when the initial ground state is a conventional superfluid, while the left panel illustrates the ground state pairing gap as a function of the onsite interaction. When the evolution time following the quench is sufficiently long, the pairing gap typically approaches a finite value, suggesting that the system relaxes to a steady state. For large U f , the pairing gap of the steady state is close to that in the ground state of the quenched Hamiltonian. For small U f , the steady state pairing gap and that of the ground state at U f are noticeably different. In comparison, we demonstrate in Fig. 3 the quench dynamics where the system is initially in a TSF state. Different from the previous case, strong oscillations of |∆(t)| are observed in the quenched state when U f is deep in the conventional superfluid regime. We observe no significant damping in the oscillations within the longest evolution time scale that we can numerically implement.
Despite the difference in the detailed dynamics, for both of these two cases, we may numerically identify a critical final interaction U d . As the quench parameter U f crosses U d , the pairing gap |∆| in the long-time limit features a jump. This is demonstrated in Fig. 4, where we plot the typical pairing gap |∆| inf in the long-time limit as a function of U f . Here, |∆| inf is the asymptotic value of the pairing gap if the system approaches a steady state in the long-time limit; and is the time average of pairing gap in many periods in the case of strong oscillations. We associate this critical parameter U d with a dynamical phase transition. For a non-topological initial state, U d is numerically close but not equal to the first-order phase boundary in the ground state; while for a topological initial state, U d lies in the conventional superfluid regime, whose deviation from the topological phase boundary is much clearer. We have further checked numerically that U d is not sensitive to the initial interaction U i .
A particularly interesting observation here is that in either case, for U f < U d , |∆| inf is smaller than the critical value ∆ c that lies on the boundary of the topological phase transition of the initial state. In contrast, when U f > U d , |∆| inf is larger than ∆ c . This strongly suggests that the dynamical phase transition is also a topological one.
Edge states in dynamical processes.-As the quenched state in the long-time limit is quite different from the ground state of the quenched Hamiltonian, the characterization of the topology of the quenched state is in general a highly non-trivial problem. Here, we propose to invoke the bulk-edge correspondence and examine the existence of edge states in the long-time limit [41]. In heterostructure nanowires, a widely studied detection method for the Majorana edge state is by coupling the nanowire to a normal metal lead and measuring the tunneling conductance at zero bias. In the presence of a Majorana edge state, the resulting conductance is quantized to 2e 2 /h [42]. Here, we borrow this idea to theoretically study the existence of edges states in the quench dynamics. In particular, we examined the topology of the quenched state in the longtime limit. In practice, the transport measurement on cold atoms has been realized in recent experiments [43], suggesting the possibility of carrying out similar measurements in the future.
We consider a spin-orbit coupled Fermi gas loaded into a lattice of L sites, labeled as 0 to L − 1, and with open boundaries. A non-interacting Fermi gas in a half-infinite chain serves as a reservoir (lead), and is coupled to the edge site of the lattice gas labeled 0. The Hamiltonian of the reservoir is then: , with the coupling Hamiltonian:Ĥ V = g V σ (ĉ † −1,σĉ 0σ + h.c.), whereĉ 0σ and c † −1,σ are field operators at the end of the interacting and noninteracting lattice respectively. The full Hamiltonian is thenĤ tot (t) =Ĥ ef f (t) +Ĥ res +Ĥ V . Setting the chemical potential of the reservoir to µ R while noticing that the gap of the superfluid is always symmetric with respect to the Fermi surface, we may define the zero-bias differential conductance at the junction as: . (4) Here, the evolution of field operatorsĉ † −1,σ (t) andĉ 0,σ (t) is governed byĤ tot (t), where the system and the reservoir are decoupled, and · · · denotes the expectation value with respect to the initial state. We calculate G(t) by the numerical operator method [36,[44][45][46], and take the limit L → ∞ before t → ∞. Typically, we find that G(t) relaxes to G very fast. Fig. 5 shows the zero bias conductance as a function of quenched interaction U f , where the unit of conductance is 2/h [43,47]. When the initial state is a normal superfluid (the left panel), the conductance drops continuously to zero at the dynamical critical point U d . When the initial state is a topological superfluid (the right panel), the conductance drops abruptly from a value close to 2/h to approximately zero at U d . In either case, the behaviors of the conductance are qualitatively different when U f is on different sides of U d . For U f < U d , the conductance is finite and has a peak which may touch 2/h at an appropriate U f . While for U f > U d , the conductance is always close to zero as U f varies. These results strongly suggest the existence of edge states for quenches with U f < U d , regardless of the initial state. Thus, U d marks the onset of a dynamical topological phase transition. Typically, the dynamical critical point U d is close to the equilibrium phase boundary U c : in the left panel of Fig. 5, U c ≈ U d = 7.25g s ; while in the right panel they are U d = 8.15g s and U c ≈ 7.6g s .
Although the existence of edge states in the quench is not sensitive to the initial condition, the conductance in Fig. 5 behaves differently for U i < U c and U i > U c . For U i < U c , G reaches the quantized value at U f = U i (no quench), and is close to 2/h in a wide range of U f near U i , suggesting that the edge modes in quenched states are well protected before U f crosses the critical point. In contrast, for U i > U c , the emergence of edge states is accompanied by the crossing of the critical U c of U f in the quenched Hamiltonian. As the bulk gap in the ground state closes at U c , the inevitable excitations in the quenched state make the signature of the edge states less obvious.
Conclusions.-In summary, we find a dynamical phase transition in cold fermions with synthetic spin-orbit coupling under an interaction quench. This transition is characterized by an abrupt change of asymptotic pairing gap, which also crosses the boundary distinguishing the topological and the trivial ground states. The edge state analysis strongly supports that this dynamical transition is a topological one, with the topological properties of the quenched state sensitively relying on the quench parameters. With the recent achievement of synthetic spin-orbit coupling and the high controllability of cold atomic gases, the dynamical topological phase transitions discussed here may be probed in future experiments.

Differential equations
In this section, we derive the equations of motion for the expectation values of the operators. As an example, we consider the expectation value ofτ k in the Schrödinger's picture: where |ψ(t) denotes the wave function at time t. The Schrödinger's equation gives: The equation of motion for τ k (t) is then: The commutator on the right-hand side consists of six possible operators: {ν k↑ ,ν k↓ ,n k↑ ,n k↓ ,τ k ,σ k }. By calculating the derivative of their expectation values to |ψ(t) , we get a set of closed-equations for the expectation values of the operators: The pairing gap ∆(t) in Eq. (8) can be expressed as which can be solved self-consistently from Eq. (8) and Eq. (9).

Calculation of tunneling conductances
To calculate the tunneling conductance, we employ the numerical operator method first introduced in Ref. [44] for the real time Kondo problem, which was then extended to time-dependent systems [45] and topological superconductors [46].
The full Hamiltonian describing the reservoir and the superfluid is expressed aŝ We define the synonymsd jσ1 =ĉ † jσ andd jσ0 =ĉ jσ and re-express the conductance as .
(11) We first calculate the operatorsd −1,σ,1 (t) andd 0,σ,0 (t) in the Heisenberg picture defined aŝ where α is the abbreviation of (j, σ, s) (s = 0, 1),Ĥ tot is the time-dependent Hamiltonian, N is the total number of steps, t 0 = 0 is the initial time, and the intermediate time t n = nτ , where the time step τ = t/N . The calculation of Eq. (12) is divided into N steps with the time evolving from t N −1 to t 0 . In the n-th step, we undress the pair of operators e iĤtot(tN−n)τ and e −iĤtot(tN−n)τ fromd α . SinceĤ tot is quadratic, in principle, the intermediate results at each step can be written asd withd α (t n ) defined aŝ d α (t n ) =e iĤtot(tN−n)τ · · · e iĤtot(tN−1)τd α × e −iĤtot(tN−1)τ · · · e −iĤtot(tN−n)τ , and W α,α ′ (t n ) denoting the propagator. At n = N , we getd To calculate the propagator W α,α ′ (t N ), we derive an iterative relation. According to Eq. (14), we havễ Substituting Eq. (13) into the equation above, we find The calculation of these commutators is straightforward with the results Substituting Eq. (18) into Eq. (17) and comparing the left and the right sides, we obtain the iterative relation Here we drop terms of the order O(τ 3 ). In principle, the error caused by the discretization of time can be made arbitrarily small by letting τ → 0. In practice, we adaptively choose τ to make the discretization error negligible. According to Eq. (19), we start from n = 0 when W α,α ′ (t 0 ) = δ α,α ′ , and iteratively work out W α,α ′ (t n ) until t N = t. This approach typically works well, due to the critical fact that the matrix G α,β has only very few non-zero elements, such that the number of non-zero W α,β which needs to be kept at each step increases slowly according to Eq. (19). However, when N is very large, which is inevitable since we must choose a small τ and at the same time a large t to obtain the asymptotic limit of G(t), there will be too many non-zero W α,β which cannot be simultaneously kept. To work around this, we apply a truncation scheme and keep only a fixed number of nonzero propagators (the number is denoted by M ) with the largest magnitudes at each step. This truncation scheme is necessary for obtaining the steady limit of G(t). We decide the value of M adaptively, i.e., set an original M and increase it until the desired precision is obtained.
The logic in doing so is that the truncation error goes to zero in the limit M → ∞.
Substituting the expression ofd α (t) into Eq. (11), we have Next we calculate the derivative of correlation functions. The correlation functions of the initial state are non-zero only if the two field operators are both in the superfluid or both in the reservoir, since the superfluid and the reservoir are decoupled at the initial time. Considering that the initial state of the superfluid located on sites from 0 to L − 1 is independent of µ R , we immediately get as j, j ′ ≥ 0. As j, j ′ < 0, we notice that the reservoir is in the ground state of a free Fermi gas at the chemical potential µ R , and get with θ = arccos(−µ R /2g l ) and g l the hopping between two neighbor sites of the reservoir. Finally we find