A Memory of Majorana Modes through Quantum Quench

We study the sudden quench of a one-dimensional p-wave superconductor through its topological signature in the entanglement spectrum. We show that the long-time evolution of the system and its topological characterization depend on a pseudomagnetic field Reff(k). Furthermore, Reff(k) connects both the initial and the final Hamiltonians, hence exhibiting a memory effect. In particular, we explore the robustness of the Majorana zero-mode and identify the parameter space in which the Majorana zero-mode can revive in the infinite-time limit.

In this report we fuse two subjects by exploring the quench dynamics of Majorana zero-modes of a 1D p-wave superconductor under the entanglement measurement process, thereby offer a quantum information perspective of the manipulation of topological systems and the robustness of the Majorana zero-modes under sudden quench. We find that the topology of the infinite-time behavior can be determined by the properties of a pseudomagnetic field R eff , which connects both the initial and the final Hamiltonians, hence exhibiting a memory effect. In general, a quench across any phase boundary will not give rise to Majorana zero-modes. Surprisingly, the quench within the same topologically nontrivial phase may also lead to the loss of the Majorana modes. We provide an equation describing a critical surface in the parameter space of µ /2, µ′  /2, and ∆∆′ ∼∼ . For quenches whose parameters lie above the critical surface, the initial Majorana zero-modes can revive in the long-time limit.

Results
The 1D p-wave superconducting system of spinless fermions proposed by Kitaev 3  where σ = (α x , σ y , σ z ) are Pauli matrices, and R(k) = (0, − Δ sin k, t cos k + μ/2) is the pseudomagnetic field. The one-particle energy spectrum is simply ε µ = ± =± + + ∆ k R k t k k ( ) 2 ( ) (2 cos ) 4 sin 2 2 2 . The spinless p-wave superconductor (1) breaks the time-reversal symmetry but preserves the particle-hole symmetry therefore it belongs to the class D according to the classification of topological insulators and superconductors; it can be characterized by a Z 2 topological invariant 30,31 .

is described by the Hamiltonian
The topological characterization has a simple graphical interpretation 26 . If the closed loop  of R(k) in the R y -R z plane encircles the origin, zero-energy edge states exist and the system is in the nontrivial phase; otherwise, the loop can be continuously deformed to a point with the bulk gap perserved, hence the system is trivial. We plot the phase diagram of the p-wave superconductor in Fig. 1 Fig. 2. λ m is known as the one-particle entanglement spectum (OPES). For total system AB the bulk correlation function matrix is a 2 × 2 matrix in the Fourier space 27 . Specifically for the Hamiltonian (2) one has Dashed arrow: quench process from (0.5, 2) to (− 0.5, 0.1) [to be discussed in Figs 4(b) and 5(c)]. Insets: Representative traces of R(k) in the R y -R z plane. In the topological phases I and II, R(k) encircles the origin (red dots), while in the trivial phases III and IV, R(k) does not encircle the origin.
Scientific RepoRts | 6:29172 | DOI: 10.1038/srep29172 where k ∈ (− π, π ]. To connect the CFM and the edge state recall that in one dimension the nontrivial Berry phase for the periodic systems in their thermodynamic limit confirms the topological nature of the systems. According to the bulk-edge correspondence 37 where i, j ∈ A. Therefore the entanglement Hamiltonian shares the same eigenvalues of the block correlation function matrix C. From the discussion above we know that the Berry phase of C in its thermodynamic limit is also determined by the pseud-magnetic field R. According the bulk-edge correspondence we conclude that there appear two edge states for the block correlation function matrix C due to the natural open boundaries of the block and therefore the entanglement spectra share the same edge states. On the other hand, by comparing H(k) (Eq. 2) and C(k) (Eq. 3) it is evident that they have the same topological characterization. By applying the bulk-edge correspondence to both the original and entanglement Hamiltonian one see that the zero energy edge modes of the original Hamiltonian correspond to the λ m = 1/2 eigenvalues of the correlation function matrix. Therefore, in the topological phases I and II, the signature of Majorana zero-modes is the degenerate eigenvalues λ m = 1/2 in the OPES.
The Majorana zero-modes play an important role in the entanglement between the subsystem A with its environment B. We can calculate the entanglement entropy E S for the partition as = ∑ E S S m m where S m = − λ m log 2 λ m − (1 − λ m ) log 2 (1 − λ m ). The pair of Majorana modes with λ m = 1/2 contribute the maximal entanglement S m = 1. Hence they are known as the topological maximally-entangled states (tMES) 28,29 .
Different from non-integrable models, which will be thermalized at infinite time, the quench dynamics of integrable models has become an important topic since such a bulk system will not be thermalized but reach a steady state decided by the initial condition described by general Gibbs ensemble (GGE). On the other hand, the Majorana edge modes are tMES, robust against perturbations, it is interesting to question how a sudden quench affect the Majorana zero-modes. Naively, if we quench from a topological phase to a trivial one, the Majorana modes may evolve into the bulk, mix with bulk modes, and disappear eventually. What happens, then, if we quench from a trivial phase to a topological one, or from a topological phase to a different one? Will the static information in the final Hamiltonian dictate the dynamical state in the long-time limit?
To answer these questions we perform a sudden quench to the p-wave superconductor at t = 0 by switching Δ and μ at t < 0 to Δ ′ and μ′ at t > 0. This changes the Hamiltonian from H to H′ and, correspondingly, R to R′ . We then calculate the time-dependent OPES λ m (t) by diagonalizing the time- We first consider the quench processes across phase boundaries and focus on λ m close to 1/2 as we are primarily interested in the fate of the Majorana zero-modes. Figure 3(a-c) show the time evolution of the OPES near 1/2 by suddenly quenching the systems from phases II, III, and IV, respectively, to phase I. We find that the Majorana zero-modes fail to appear after a sufficiently long time, regardless of the topological properties of the original state. In other words, the quench of the topological systems with the Majorana edge modes cannot be thermalized. For comparison, Fig. 3(d-f) show the time evolution of the OPES for the sudden quench from phase I to phases II, III, and IV, respectively. The degenerate eigenvalues λ m = 1/2 persist for some time before they split and relax to separate values, depending on the final Hamiltonian. In other words, the Majorana zero-modes before a sudden quench are destroyed eventually.
To further confirm the topology of the steady states of the quench process in the infinite-time limit, we analyse the time-dependent CFM in Fourier space G(k, t) as illustrated in the method section: Time-dependent correlation matrix. We find it can be described by a pseudomagnetic field R(k, t) (15) through the relation G(k, t) = [1 − R(k, t) · σ]/2. In the infinite-time limit the sinusoidal time dependence dephases away and G(k, t = ∞ ) depends only on the effective pseudomagnetic field eff Therefore, the topology of the steady state is determined by both the initial and final Hamiltonians: the quench dynamics has a memory of the initial Hamiltonian, albeit entangled with the final Hamiltonian. The insets in We confirm that all the traces pass through the origin, indicating that the Majorana zero-modes are not stable in the infinite-time limit. Will Majorana zero-modes persist if we quench between two Hamiltonians in the same topological phase? Given the fact that the edge states cannot thermalize, the naive answer of yes needs to be examined. In Fig. 4(a,b) we show the OPES evolution near 1/2 for two quantum quenches both within phase I, as indicated in the phase diagram in Fig. 1. Surprisingly, we find that the Majorana zero-modes reappear in the steady state in the infinite-time limit in Fig. 4(a), while disappear in Fig. 4(b). This contrasts to the persistence of the edge modes in a dimerized chain, where the edge mode is an electron instead of Majorana zero-modes 29 . We also plot the corresponding R eff (k) in the insets of Fig. 4(a,b). R eff (k) encircles the origin in the former case, which confirms the persistent memory of the Majorana modes after a sudden quench. In sharp contrast, R eff (k) passes through the origin in the latter case, which is consistent with the loss of the memory of the Majorana modes. Note that the quench processes within the phase II shall behave similarly.
To support our findings, we further analyze the eigenstates of the after a sudden quench. In Fig. 5 we plot the probability sum ψ = ∑ | | = P sum j j 1,2 2 of the two states ψ 1,2 whose eigenvalues are closest to 1/2 at t = 0, 9, and 99. As expected, P sum exhibits sharp peaks due to the presence of the Majorana edge modes. The only case that such peaks survive [ Fig. 5(b)] is the quench process within phase I as discussed in Fig. 4(a). For comparison, such peaks dissolve into the bulk when we quenches the system to a different topological phase

Discussion
We have shown in the result section that the maximally entangled states (or Majarona zero modes) can still disappear at t → ∞ even for the quench between the same topological phase. As hinted from (5), this is due to the condition ⋅ ′ =k k R R ( ) ( ) 0 happens at some k. Explicitly, this means where µ µ µ = + ′    ( )/2 and δµ µ µ = − ′    . Figure 6(a) shows the critical surface below which the equality can be satisfied for some k; while above the surface R eff (k) encircles the origin (to ensure that both the initial and final Hamiltonians are in the topological phase, we also need µ  , µ′ ≤  2), hence the Majorana zero-modes are memorized in the long-time limit. To understand why ⋅ ′k k R R ( ) ( ) can vanish at some k, we find it is instructive to consider the velocity of the R-vectors. Intuitively, if the velocity profiles of two R-vectors are similar then it is impossible to have ⋅ ′ =k k R R ( ) ( ) 0. In contrast, if the maximal angular velocity of R(k) and R′ (k) occur at different k points then it is possible to have ⋅ ′ =k k R R ( ) ( ) 0. For example, if R(k) rotates rapidly while R′ (k) rotates slowly then they should become perpendicular at some k. As shown in the method section: Maximal angular velocity and band gap, we find the position of the band gap and the maximal of the angular speed ω(k) = dR(k)/dk are at the same k (can be slightly different for some special parameter region). Since the position of the band gap can be tuned by the chemical potential μ, one can make the maximal angular velocity of R(k) and R′ (k) occur at different k by tuning the chemical potential of the initial and final Hamiltonian. Thus, for the quench within the same topologically non-trivial phases (both R and R′ vectors rotate clockwise or counter-clockwise), it is still possible to make R rotates slowly (rapidly) while R′ rotates rapidly (slowly) at a particular k, leading to ⋅ ′ =k k R R ( ) ( ) 0. A simple physical picture to understand this phenomenon is as follows: if the superconducting  2 of the two eigenstates of the entanglement Hamiltonian whose eigenvalues are closest to 1/2 at t = 0, 9, and 99 after a sudden quench (a) from I (0.5, 2) to II (0.5, − 2), (b) from I (0.5, 2) to I (0.5, 1) as in Fig. 4(a), and (c) from I (0.5, 2) to I (− 0.5, 0.1) as in Fig. 4(b). gaps of the initial and final Hamiltonians have no energy overlap, the Majorana modes are suppressed due to the mismatch of the corresponding single-particle states.
For the early time behavior, one may ask what determines the finite survival time of the degenerate λ m = 1/2 in Fig. 3(d-f). We can imagine that upon quench quasiparticles are generate in the bulk and propagate at a maximum velocity v max = [∂ ε(k)/∂ k] max . Therefore, at = ★ T L v /(2 ) max the Majorana zero-modes at the boundaries of the entanglement cut can exchange information and hence the degenerate levels in the OPES start to split. Chaotic oscillations then emerge in the entanglement spectrum due to the complex processes of quasiparticle interference and decoherence.
Finally, we discuss the robustness of our result under the perturbative Coulomb interactions: = ∑ − . − .
for the final Hamiltonian. Within the mean-field approach, the effects of the interaction is to renormalize the hopping, paring and chemical potential, resulting an effective mean-field Hamiltonian of the form of the un-perturbed Hamiltonian. As a result, the renormalized hopping are the pseudomagnetic fields of ′ H M with the form given by (2). To understand the effect of the mean-field approach, the same example as Fig. 4(a) (same as Fig. 5(b)) is considered in the follow calculation. Considering V = 0.1 for the final Hamiltonian, we obtain µ ′ ′ ≈ . corresponding to the movement from the initial black point to the red one in Fig. 7. However this change under perturbation does not cross the critical surface which is plotted as the blue curve in Fig. 7. The results are confirmed by Fig. 8, where the two OPES ψ 1 and ψ 2 near 1/2 and the probability sum ψ = ∑ | | = P sum j j 1,2 2 at t = 0, 9, and 99 for the red point are shown. As far as the final Hamiltonian is far away from the critical surface, the perturbations of the interaction do not affect the revival of the tMES. The same considerations also apply for the perturbation of the pairing and chemical potential.
In summary, we study the quench dynamics of a 1D p-wave superconductor using the OPES. We find that the system reach a final steady state whose topology can be determined by an effective pseudomagnetic field R eff (k) (5). As expected, sudden quenches from a topological phase to a trivial phase destroy the Majorana edge modes. However, the memory of the Majorana modes will also be lost if we quench the system to a different topological phase, or if the superconducting gaps before and after a sudden quench do not match. When both topological and energetic criteria are satisfied, the Majorana zero-modes will return after sufficiently long time.

Time-dependent correlation function. To obtain the time-dependent correlation function matrix
, one has to diagonalize H (2) using iH t k k iH t can be calculated by canonically transforming α twice to the proper operators using (9) and obtain  For zero temperature, ρ is the density matrix of the ground state, therefore Substituting (13) and (9) into (12) and using the properties of Pauli matrices (R · σ) (R′ · σ) = R · R′ + i(R × R′ ) · σ, we end up with a simple form where the effective pseudomagnetic field

Maximal angular velocity and band gap.
Here we are going to determine the momentum k 0 at which ω(k) reaches global maxima and the momentum k g at which the band gap of the dispersion ε(k) is located. Then, we illustrate an interesting relation between angular velocity ω(k) and the dispersion ε(k): that one has k 0 = k g in some parameter regimes and k 0 ≈ k g in other parameter regimes as shown in Fig. 9. This provides an intuitive picture on why the topology of the steady state is determined by the initial and final Hamiltonian as discussed in the discussion section. In the following, we will use the dimensionless parameters µ µ ≡ t /2 , ∆ = ∆ t / and set t = 1 as our unit of energy. Furthermore, we will assume that µ < 1 since we only want to consider the region of topologically non-trivial phases.
Extrema of the angular velocity ω(k). To identify the momentum k 0 at which the angular velocity ω(k) reaches its maximum, we start from the pseudomagnetic field In terms of the dimensionless constants defined above it can be rewrite as (0, sin( ), cos( )) (0, sin( ), cos( )), (17) where α = − ∆ (1 )/2, and β = + ∆ (1 )/2. This shows that R(k) can be decomposed into as constant vector in the z-direction, and rotating and counter-rotating vectors in the y-z-plane. Let θ(k) be the angle between R(k) and the z-axis, then one has Define the angular velocity as ω(k) = dθ(k)/dk. By differentiating the equation above we obtain the expression of angular velocity The k 0 at which the angular velocity ω(k) reaches its maximum should satisfy = and ε(k) is the dispersion. We note that due to the finite band gap, the denominator is always non-zero. From Eq. (20) we find k 0 = 0 or k = π or Note that in Eq. (22) we have discarded the solution with positive sign in front of the square root because it does not lead to a valid solution.
Extrema of the dispersion ε(k). To identify the momentum k g at which the band gap is located we start from the dispersion The Solving for k g we find k g = 0 or k g = π or Relationship between k 0 and k g . With the solution of k 0 , the location of the maximal angular velocity ω(k) and k g the location of the gap of the dispersion ε(k) we are now is position to illustrate the close relation between k 0 and k g . In Fig. 9 we show the difference |k 0 − k g | as a function of ∆ and µ. We find that k 0 and k g are almost the same in most of the parameter regime. In the following we provide more insight on why this happens via analyzing the solutions in various limits.
• Case 1: k o = k g = 0 or π Consider first the case where both Eqs (22) and (27) do not lead to a valid solution. In this case one has k 0 = 0 or k 0 = π and similarly for the k g . Since and Scientific RepoRts | 6:29172 | DOI: 10.1038/srep29172 one finds that if µ < 0 then k 0 = k g = 0, while if µ > 0 then k 0 = k g = π. Therefore the location of the band gap and the maximum angular velocity are at the same k, as shown in Fig. 10(a,b) for µ > 0 and µ < 0 respectively.
• Case 2: Strong paring; ∆ > 1 2 Consider next the case of strong pairing, where ∆ > 1 2 . Due to the restriction of µ < 1 it is straightforward to show that in this case the square root of Eq. (22) is always positive. It is then instructive to consider the two limits of ∆: ∆ → 1 and ∆ → ∞.In limit of ∆ → 1, one has Consequently both Eqs (22) and (27) have no valid solution and one falls back to the scenario of case 1.
In the other limit of ∆ → ∞, Eq. (22) becomes  In Fig. 11 we plot the absolute value of Eq. (33) as a function of µ and show that there is always a solution for k 0 within the parameter regime. Due to the µ dependence of k 0 , in this limit one has k 0 ≠ k g . However, in this limit one also finds By comparing Eqs (29) and (35), one finds that in this limit Eq. (33) defines a minimum of the angular speed instead of a maximum. Similarly, by comparing Eqs (30) and (36) one finds that the band gap is not located at the momentum defined by Eq. (34). Therefore, we again go back to the scenario of case 1 as shown in Fig. 10(c). This means that the k 0 and k g are asymptotically close to each other when the paring goes to zero. On the other hand, in this limit we also have  (30) and (40). Thus, the k 0 and k g defined in Eqs (22) and (27) are almost the same as shown in Fig. 10(d).
On the other hand, for limit of µ ∆ → − 1 2 2 , Eqs (22) and (27) become which gives no valid solution for k 0 and k g due to the restriction of µ < 1 and thus falls back to scenario of case 1. Figure 11. The absolute value of Eq. (33) for the allowed parameter regime µ < 1. There is always a solution for 0 < |k 0 | < π.