Quench dynamics and parity blocking in Majorana wires

Of late, two different concepts in quantum many-body theory have elicited a surge of active research, partly stemming from experimental advances in condensed matter and cold atomic systemsthe concepts of quench dynamics [1–51] and topological order [52–57]. Quenching, or ramping, concerns initializing a system in its equilibrium configuration at some point in parameter space followed by inducing non-equilibrium behavior via dynamic tuning of one of the parameters. When the tuning occurs through a critical point separating two phases of matter, no matter how slow the tuning rate 1/τ, the diverging time scale associated with the critical point and critical exponent z always results in out-of-equilibrium dynamics in its vicinity. The quantum version of such Kibble–Zurek physics, initially studied as a thermal quench during the formation of the early Universe [4–7], offers a mine of valuable information about the critical point in question.

In the realm of topological systems, while quantum Hall systems have been hailed for their topological properties for over three decades [52], the recent attention on other systems has also been spectacular [53–57]. On the theoretical front, among others, two paradigm low-dimensional models have been avidly studied for their topological propertiesthe Majorana wire proposed by Kitaev [58, 59], which is effectively a lattice version of a spinless p-wave superconducting wire, and the two-dimensional Kitaev honeycomb model [60].

Here we explore the synergy of these two concepts, namely quench dynamics⁴ and topological order, a rich study that is still in its infancy. Such a marriage is exciting from at least two perspectivescan quench dynamics act as a probe for topological order? Can the existence of topological order lead to a different realm in non-equilibrium dynamics? Studies of such synergy [30, 31, 38, 39, 61] have just begun to explore diverse and exciting phenomena with regards to dynamic evolution of topological features. In previous work by two of the authors of this article and co-workers, reference [62], the term 'topological blocking' was coined with regards to the role played by a highlighting feature of topological order—ground state degeneracies—in quench dynamics. It is known that a system can have several topological sectors which are associated with these degeneracies and are distinguished from each other by an invariant based on a discrete global symmetry [63, 64]. In quenching between a topological and non-topological phase, if the ground states in the two phases belong to different topological sectors and the Hamiltonian commutes with the global symmetry at all times, the system never reaches the ground state of the final phase. As a result, the usual expectation that if the quench is sufficiently slow (i.e., almost adiabatic), the system always remains in the instantaneous ground state is violated. This topological blocking effect was demonstrated in [62] for the Majorana wire and the Kitaev honeycomb model constrained to periodic boundary conditions.

The goal here is to explore quench dynamics in a finite-sized Majorana wire having open boundary conditions, a system that has come into the limelight for topological features that we expect to affect dynamics in a profound way. These features concern the presence of isolated zero energy Majorana fermion bound states at the wire ends within the topological phase; their possible experimental detection in the context of spin–orbit coupled wires [65–70] has garnered much attention in terms of fundamental physics as well as implications for topological quantum computation [71]. In the thermodynamic limit, these Majorana end modes together form a Dirac fermion state that can either be occupied or empty. Thus, degenerate topological sectors are identified by fermion parity. In the finite-sized system, tunnel coupling between these end modes splits the degeneracy in a manner that can be tuned by changing the parameters of the system. Here we explore the quench dynamics of tuning through a succession of parity flips of the ground state. The investigation involving open boundaries requires the formulation of new dynamic quantum many-body techniques, which we develop here. By investigating measures commonly studied in quench dynamics, we demonstrate that topological order drastically affects non-equilibrium behavior, the most dramatic signature stemming from quench-dependent switching of topological sectors.

Our presentation is as follows. In section 2, we present a brief description of the salient features of topological blocking and the highlights of this work, including our results regarding parity switching and blocking in a finite-sized Majorana wire. In section 3, we begin our detailed exposition by reviewing the p-wave superconducting wire given by Kitaev's lattice Hamiltonian. We outline a derivation of its bulk spectrum and description of its topological phase diagram based on the presence or absence of Majorana end modes. In section 4, based on a transfer matrix formalism, we analyze the fate of the Majorana end modes for a finite-sized lattice. We obtain an approximate form for the degeneracy splitting and exact solutions for contours in the phase diagram where the splitting vanishes. We then invoke Kitaev's argument based on Pfaffian methods to determine the ground state parity of the system and confirm that parity switches occur at these degeneracy points. In section 5, we begin our discussion of the quench dynamics associated with varying a parameter of the underlying Hamiltonian linearly in time. We summarize the known results in the case of periodic boundary conditions, in particular, the Kibble–Zurek scaling of the post-quench excitations and the topological blocking phenomenon. In section 6, we describe the real space time-dependent formalism that we use to study this problem numerically. We focus on two measures, namely the overlap between the time evolved state and the instantaneous ground state, which we refer to as adiabatic fidelity as in previous work [2], and the residual energy, which is the difference between the expectation value of the Hamiltonian in the time evolved state and the instantaneous ground state energy. In section 7, we extensively discuss our numerical results in detail, pinpointing the effect of topological blocking associated with multiple parity switches. In section 8, we present an overview of our study and connect it with related phenomena, such as the fractional Josephson effect in junctions of Majorana wires, as well as to experiments.

Here, we describe the essence of the topological features that we target with regards to quench dynamics and present the highlights of this work before embarking on a detailed exposition.

The stage is set by the concept of topological blocking, which, as mentioned in the previous section, was studied in [62]. The study involved quench dynamics in topological systems elicited by changing a parameter of the Hamiltonian to tune from one quantum phase to another. In going between a topological phase and a trivial phase, as in the Majorana wire, or between two topological phases, as in the Kitaev honeycomb model, the focus was on mismatch of degeneracies. It was shown that if the system was initialized in the ground state in a phase with higher degeneracy and tuned to one with a lower degeneracy, two to one in the former case, and four to three in the latter, the phenomenon of topological blocking would occur. Certain topological sectors characterized by topological quantum numbers, for instance, fermion parity, would inhabit the ground state in the initial phase but would have no partner in the ground state of the final phase. As a result, in tuning through a quantum critical point separating the two phases, these sectors would evolve so as to have null overlap with the final ground state no matter how slow the tuning rate, in stark contrast with Kibble–Zurek physics, where only a rate-dependent fraction of the time-evolved state overlaps with the excitation spectrum above the (gapped) final ground state. Moreover, it was shown that even if one took overlap with the instantaneous sectoral ground state, different topological sectors would show quantitatively different dynamic behavior, particularly in wave function overlap.

Here, we explore this notion of topological blocking with regards to a different but related aspect—the switching of topological sectors due to quench dynamics within a topological phase. In a Majorana wire with open boundaries, the topological degeneracy is associated with the presence of Majorana zero modes at the edges. The degeneracy is split due finite-size coupling of these edge modes, which induces fermionic parity sectors within the topological region of the phase diagram. Here we build on the notion of topological/parity blocking arising from tuning through these parity sectors.

Zoom In Zoom Out Reset image size

In figure 2, we present some of our key results. Initializing the system in the ground state corresponding to a specific on-site chemical potential, and thus some fixed parity, we sweep the chemical potential to undergo several parity switches. As a measure of how closely the time evolved state tracks the instantaneous ground state, we evaluate the wave function overlap (adiabatic fidelity) associated with these two states. As seen in figure 2, the adiabatic fidelity plummets down to zero in certain chemical potential intervals that exactly correspond to the parity switched regions. The initial ground state, while being able to track some of the dynamic evolution, is thus forced to remain within its parity sector, an attribute of the topological phase. This multiple parity blocked dynamics is a dramatic, topologically induced deviation from the continuous evolution expected in quench dynamics.

Zoom In Zoom Out Reset image size

In what follows, we detail several aspects leading up to this quench behavior, including the formulation of the Majorana wire model, parity switching due to coupling of Majorana wavefunctions and a real space formalism to compute the many-body dynamics.

The system that forms our subject of study is a lattice version of the spinless fermionic one-dimensional p-wave superconducting wire with spinless fermions, also referred to as the Kitaev chain [58]. This system can be mapped exactly to a spin-1/2 XY model in a transverse field (i.e., a magnetic field applied along the z direction) via the Jordan–Wigner transformation [72]. The parameters of the system are the nearest-neighbor hopping amplitude w, the superconducting pairing amplitude between nearest neighbors Δ, and a chemical potential μ. This model is a paradigm system for demonstrating numerous interesting topological properties, including the existence of Majorana modes at the ends of an open chain in the topological phase.

The Hamiltonian of such a system with N sites and open boundary conditions is given by

$$\begin{eqnarray}&&H=-\mathop{\sum }\limits_{n=1}^{N-1}(-wf_{n+1}^{\dagger }{{f}_{n}}+\Delta f_{n+1}^{\dagger }f_{n}^{\dagger }+{\rm h}{\rm .c}{\rm .})-\mu \mathop{\sum }\limits_{n=1}^{N}(f_{n}^{\dagger }{{f}_{n}}-1/2),\end{eqnarray} \tag{ 1 }$$

where f_m are Dirac fermion operators obeying the commutation relations $\{{{f}_{m}},{{f}_{n}}\}=0$ and $\{{{f}_{m}},f_{n}^{\dagger }\}={{\delta }_{mn}}$. We now introduce $2N$ Majorana fermion operators as ${{a}_{2n-1}}={{f}_{n}}+f_{n}^{\dagger }$ and ${{a}_{2n}}={\rm i}(f_{n}^{\dagger }-{{f}_{n}})$. These satisfy the relations $a_{m}^{\dagger }={{a}_{m}}$ and $\{{{a}_{{\rm L}}},{{a}_{m}}\}=2{{\delta }_{\operatorname{lm}}}$. In terms of these operators, the Hamiltonian takes the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} H & = & -\frac{{\rm i}}{2}\mathop{\sum }\limits_{n=1}^{N-1}\left[ (w-\Delta ){{a}_{2n-1}}{{a}_{2n+2}}-(w+\Delta ){{a}_{2n}}{{a}_{2n+1}} \right]-\frac{{\rm i}\mu }{2}\mathop{\sum }\limits_{n=1}^{N}{{a}_{2n-1}}{{a}_{2n}}. \\ \end{array}\end{eqnarray} \tag{ 2 }$$

The Hamiltonian in equation (1) can be diagonalized, up to a constant, to give the canonical form

$$\begin{eqnarray*}&&H=4\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}b_{j}^{\dagger }{{b}_{j}},\end{eqnarray*}$$

where the ${{\lambda }_{j}}$ are non-negative real numbers. This can be done through a transformation of the form

$$\begin{eqnarray}&&\bar{b}=B\bar{a},\end{eqnarray} \tag{ 3 }$$

where $\bar{a}={{({{a}_{1}},{{a}_{2}},\cdots ,{{a}_{2N}})}^{T}}$ and $\bar{b}={{({{b}_{1}},\cdots ,{{b}_{N}},b_{1}^{\dagger },\cdots ,b_{N}^{\dagger })}^{T}}$ are column vectors with $2N$ components. The $(2N)$-dimensional matrix B comprises of the eigenvectors of H and their Hermitian conjugates, belongs to the unitary group $U(2N)$, and has its determinant satisfying the property ${\rm det} (B)=\pm 1$. The energy eigenvalues of the Hamiltonian are 0 and $4{{\lambda }_{j}}$.

The phase diagram of this model is shown in figure 3. The phases I and II are topologically non-trivial and, in the thermodynamic limit, have zero energy Majorana modes bound to the ends of the wire, whereas such modes are absent in the topologically trivial phases III and IV. These Majorana modes have finite support at the ends and decay rapidly into the bulk with a decay length proportional to the reciprocal of the bulk gap. One can understand the existence of the Majorana end modes by considering the extreme limit of $w=\Delta$ and $\mu =0$. The Hamiltonian reduces to $H={\rm i}w{{\sum }_{n}}{{a}_{2n}}{{a}_{2n+1}}$. The Majorana operators a₁ and a_2N are not paired with any other operators in the system and therefore do not appear in the Hamiltonian. These isolated modes correspond to the zero energy eigenvectors localized at the ends. The existence of these modes is robust even away from this extreme limit and they only disappear with the closing of the bulk gap.

Zoom In Zoom Out Reset image size

The ground state of the system in the topological phase is thus doubly degenerate and has two zero energy eigenvalues corresponding to the Majorana modes. These Majorana modes can be combined to form a complex Dirac fermion state, which can be either empty or occupied. Hence, each of the degenerate ground states has a specific fermion parity and the system can be characterized by a related Z₂-valued topological invariant. This ground state parity will play an important role in the subsequent sections.

There are three phase boundaries, indicated by dark lines in figure 3, where the bulk gap vanishes. These are the quantum critical lines across which there is a topological phase transition. In the thermodynamic limit (infinite wire) or for a closed chain, one can transform the Hamiltonian into Fourier space, and the single particle energy spectrum takes the form

$$\begin{eqnarray}&&{{E}_{k}}=\pm \sqrt{{{(2w{\rm cos} k+\mu )}^{2}}+4{{\Delta }^{2}}{{{\rm sin} }^{2}}k}.\end{eqnarray} \tag{ 4 }$$

This spectrum has a finite superconducting gap in all the phases; the gap vanishes as one crosses one of the critical lines and reopens upon entering another phase. In the spin language, the topological phases correspond to the ferromagnetic phases of the transverse field XY model (where either the x or the y component of the spins has long range order), and the trivial phases are in the paramagnetic phase.

These characteristic features of the system, namely, the topological invariant, the spectrum of the bulk and end modes, and the wave functions of the Majorana end modes have been discussed extensively in previous work (see, for example, [34, 56, 58]). So far, all the above mentioned characteristics of the model assume the size of the system to be much larger than the decay length of the Majorana end modes. In the next section we consider the case when the Majorana modes at the two ends have a finite overlap, giving rise to consequences such as parity blocking in quench dynamics.

4.1. Tunneling between Majorana end modes

In a Majorana wire of finite length, the two Majorana end modes are no longer completely decoupled since there is some overlap between their wave functions (figure 1); the overlap shifts their energies slightly away from zero. In the extreme limit of the topological phase with $\Delta =w$, the end modes are exactly localized at the ends with no overlap between them; hence the effective Hamiltonian governing these modes is

$$\begin{eqnarray}&&{{H}_{f}}={\rm i}J{{a}_{1}}{{a}_{2N}},\end{eqnarray} \tag{ 5 }$$

with J = 0, But in general, the effective Hamiltonian has an expression in terms of the (almost) zero energy eigenvectors localized at the ends of a finite length wire [58]

$$\begin{eqnarray}&&{{H}_{f}}={\rm i}Jb^{\prime} b^{\prime\prime} ,\end{eqnarray} \tag{ 6 }$$

where $b^{\prime} ={{\sum }_{j}}(\alpha _{+}^{\prime }x_{+}^{j}+\alpha _{-}^{\prime }x_{-}^{j}){{a}_{2j-1}}$, $b^{\prime\prime} ={{\sum }_{j}}(\alpha _{+}^{{\prime\prime} }x_{+}^{-j}+\alpha _{-}^{{\prime\prime} }x_{-}^{-j}){{a}_{2j}}$. Here J is in general a function of $\Delta /w$, $\mu /w$ and N. Due to this coupling the two zero modes split in energy into a particle-hole symmetric pair of eigenenergies $E=\pm J$, and the ground state is no longer degenerate. The Majorana end modes can be combined into non-local Dirac fermions as $\tilde{c}=(b^{\prime} +{\rm i}b^{\prime\prime} )/2$ and ${{\tilde{c}}^{\dagger }}=(b^{\prime} -{\rm i}b^{\prime\prime} )/2$. The Hamiltonian can now be expressed as

$$\begin{eqnarray}&&{{H}_{f}}=J(2{{\tilde{c}}^{\dagger }}\tilde{c}-1).\end{eqnarray} \tag{ 7 }$$

The occupation number ${{\tilde{c}}^{\dagger }}\tilde{c}$ can either be zero or 1. Thus we see that the energies $\pm E$ come with corresponding eigenstates with opposite fermion parities. The sign of J decides which of these states is the ground state. The parity of the states in the bulk being fixed, the overall ground state parity is then decided by the lower one of the two split energy levels (which lie inside the bulk gap). The coupling J is a function of the chemical potential μ and oscillates, switching its sign at specific values of μ. Therefore the split energy levels cross zero at certain points in time if μ is varied linearly in time. This leads to oscillations in the overall parity of the ground state. Even though these split energy levels are exponentially smaller (for large system size N) than the energies in the bulk, the parity oscillations play a key role in the time evolution of the ground state in the quenching dynamics.

The energy splitting due to the coupling of the two Majorana end modes can be derived by evaluating the overlap between their associated wave functions. For a very long wire, the Majorana modes have zero energy; we calculate their exact wave functions by using the Heisenberg equations of motion $[H,{{a}_{n}}]=0$ and $[H,{{b}_{n}}]$, where a_n, b_n are the two Majorana operators at site n which were denoted by ${{a}_{2n-1}}$ and a_2n above. We obtain the difference equations for these operators as:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} (w+\Delta ){{b}_{n-1}}+\mu {{b}_{n}}+(w-\Delta ){{b}_{n+1}} & = & 0, \\ (w-\Delta ){{a}_{n-1}}+\mu {{a}_{n}}+(w+\Delta ){{a}_{n+1}} & = & 0, \\ \end{array}\end{eqnarray} \tag{ 8 }$$

for $2\leqslant n\leqslant N-1$.

These difference equations can be solved exactly by using z-transform methods and taking into account the form of the difference equations at the ends. The wave functions ${{\alpha }_{n}}$ and ${{\beta }_{n}}$ for the a_n and b_n modes on site n, respectively, are of the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\alpha }_{n}} & = & {{\alpha }_{0}}{{C}^{n}}\left[ {\rm cos} (\omega n)+\frac{1}{{\rm tan} \omega }{\rm sin} (\omega n) \right], \\ {{\beta }_{n}} & = & {{\beta }_{0}}{{C}^{-n}}[{\rm cos} (\omega n)+\frac{1}{{\rm tan} \omega }{\rm sin} (\omega n)], \\ \end{array}\end{eqnarray} \tag{ 9 }$$

where $C={{\left( \frac{w+\Delta }{w-\Delta } \right)}^{1/2}}$, and

$$\begin{eqnarray}&&\omega ={{{\rm tan} }^{-1}}\left( \frac{\sqrt{4{{w}^{2}}-4{{\Delta }^{2}}-{{\mu }^{2}}}}{\mu } \right).\end{eqnarray} \tag{ 10 }$$

Using these wave functions, we can now approximately calculate the energy splitting by assuming that the wire has a finite length N and computing the overlap of the wave functions ${{\alpha }_{{\rm L}}}$ at the left end and ${{\alpha }_{{\rm R}}}$ at the right end. This gives an expression of the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\alpha }_{{\rm L}}}{{\alpha }_{{\rm R}}} & = & \frac{{{\alpha }_{0{\rm L}}}{{\alpha }_{0{\rm R}}}{{a}^{N+1}}}{4{{{\rm sin} }^{3}}\omega }[2{\rm sin} (\omega N) \\ {} & {} & +N\{{\rm sin} [(N+2)\omega ]-{\rm sin} [(N+4)\omega ]\}]. \\ \end{array}\end{eqnarray} \tag{ 11 }$$

Figure 4(a) shows the oscillations in the splitting as calculated above. This is an approximate calculation because we have assumed the energies to be zero in the Heisenberg equations of motion for the Majorana operators and then calculated their overlap (which shifts the energies slightly away from zero). But one can see that it qualitatively agrees with the exact numerical calculation in figure 4(b). Figure 4(c) shows the oscillations of the ground state parity as a function of μ which is linearly varied with time. One sees an excellent correspondence with the zero crossings of the energy splitting in figure 4(b). Further, although these oscillations appear due to the degeneracy splitting in the topological phase, they do not exist in the entire topological region in the phase diagram. To understand this, we stress the fact that the key ingredient in getting these oscillations is the oscillatory component in the wave functions of the Majorana modes given in equations (9). One can see that the oscillatory functions ${\rm sin} (\omega n)$ and ${\rm cos} (\omega n)$ become hyperbolic if ω given by equation (10) becomes imaginary. The boundary at which this happens is given by the circle ${{\mu }^{2}}=4({{w}^{2}}-{{\Delta }^{2}})$. Beyond this circle the Majorana wave functions have only a decaying (but not oscillatory) component. This implies that there would not be any parity oscillations in this region. This can be clearly seen in figure 5(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We remark here that expressions for the energy splitting of the Majorana end modes are known for a continuum model [73, 74]. These results are consistent with the splitting being both oscillatory and decaying exponentially with increasing length.

Calculation of ground state fermion parity:to obtain a rigorous knowledge of the parity and its switching as a function of the chemical potential, as has been used in other Majorana wire contexts [75], we employ the measure introduced in Kitaev's well-known work [58]. Given a Hamiltonian of the form in equation (1), the transformation B, which reduces the Hamiltonian to the canonical form, can be represented as a conjugation by a parity preserving unitary operator if B has the form $B={{e}^{D}}$ i.e., if ${\rm det} (B)=1$. Otherwise B changes the parity. Therefore, the parity of the system is given by

$$\begin{eqnarray}&&P(H)={\rm sgn} [{\rm det} (B)].\end{eqnarray} \tag{ 12 }$$

In appendix B we illustrate this result with a simple problem of a two-site effective Hamiltonian. This illustration is particularly useful in our case since the effective Hamiltonian for only the coupled Majorana modes is in fact a two-site problem given by ${{H}_{f}}={\rm i}Jb^{\prime} b^{\prime\prime}$. Within the topological phase, the dynamics of only these end modes and their associated Dirac fermions determine the overall parity as we saw above.

In terms of the Majorana operators, the parity of a N-site system is given by

$$\begin{eqnarray}&&P\;=\;\mathop{\prod }\limits_{j=1}^{N}\;(-{\rm i}{{a}_{2j-1}}{{a}_{2}}j)\;=\;\mathop{\prod }\limits_{j=1}^{N}\;(1-2f_{j}^{\dagger }{{f}_{j}}).\end{eqnarray} \tag{ 13 }$$

We note that P is both Hermitian and unitary, and it commutes with the Hamiltonian in equation (2); since ${{P}^{2}}=I$, the eigenvalues of P must be ± 1. In the extreme case of $w=\Delta$ for the Majorana wire in equation (2), the terms of the form ${{a}_{2j}}{{a}_{2j+1}}$ with $j=1,\cdots ,N-1$ are equal to 1. Thus only the term ${{a}_{1}}{{a}_{2N}}$ remains, which is in fact the term in the effective Hamiltonian. Now, within the topological phase, small deformations of the parameters should not change this fact. Since the parity operator also commutes with the Hamiltonian, we can see that in the topological phase the end modes alone determine the parity of the ground state. The parity equivalence of ground state sectors has been studied in [76].

4.2. Exact expression for zero energy contours

While we obtained an approximate result for the tunneling amplitude that splits the degeneracy between Majorana end modes in the previous subsection, a formalism involving transfer matrices [34, 77] enables to track the exact points in the parameter space at which this energy changes sign, restoring the zero energy degeneracy and resulting in a parity switch. Previous work has presented similar derivations and results using the equivalent method of chiral decomposition [78].

For a mode with energy exactly equal to zero, equations (8) and (9) are applicable. We see that the a and b modes are decoupled; for definiteness, let us consider a zero energy mode involving the a_n's. Equation (8) shows that for $2\leqslant n\leqslant N-1$, the a_n's are related to each other by a transfer matrix M_a,

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \left( \begin{array}{ccccccccccccccc} {{a}_{n+1}} \\ {{a}_{n}} \\ \end{array} \right) & = & \;{{M}_{a}}\left( \begin{array}{ccccccccccccccc} {{a}_{n}} \\ {{a}_{n-1}} \\ \end{array} \right), \\ {\rm where}\;\;{{M}_{a}} & = & \left( \begin{array}{ccccccccccccccc} -\frac{\mu }{w+\Delta } & -\frac{w-\Delta }{w+\Delta } \\ 1 & 0 \\ \end{array} \right). \\ \end{array}\end{eqnarray} \tag{ 14 }$$

If ${{\lambda }_{1}}$ and ${{\lambda }_{2}}$ are the two eigenvalues of M_a, the general solution of equation (8) is ${{a}_{n}}={{c}_{1}}\lambda _{1}^{n}+{{c}_{2}}\lambda _{2}^{n}$. Next we note that the boundary equations become satisfied if we add fictitious sites with n = 0 and $n=N+1$ at the two ends of the system and demand that ${{a}_{0}}={{a}_{N+1}}=0$. This is possible if and only if

$$\begin{eqnarray}&&{{\left( \frac{{{\lambda }_{1}}}{{{\lambda }_{2}}} \right)}^{N+1}}\;=\;1.\end{eqnarray} \tag{ 15 }$$

We can show that equation (15) holds and that hence there is a zero energy mode if either

• (i)  
  the parameters lie insider the circular region ${{\mu }^{2}}+4{{\Delta }^{2}}\lt 4{{w}^{2}}$, and

  $$\begin{eqnarray}&&{{\left( \frac{\mu +{\rm i}\sqrt{4{{w}^{2}}-4{{\Delta }^{2}}-{{\mu }^{2}}}}{\mu -{\rm i}\sqrt{4{{w}^{2}}-4{{\Delta }^{2}}-{{\mu }^{2}}}} \right)}^{N+1}}\;=\;1,\end{eqnarray} \tag{ 16 }$$

  or
• (ii)  
  $\mu =0$ (which implies that ${{\lambda }_{1}}/{{\lambda }_{2}}=-1$ regardless of the relative values of w and Δ) and N is odd.

We get the same conditions if we look for a zero energy mode involving the b_n's.

In terms of ω defined in equations (10) and (16) is equivalent to saying that $(N+1)\omega$ is an integer multiple of π, namely, that ${\rm sin} [(N+1)\omega ]=0$. We see that this differs somewhat from the approximate condition that ${{\alpha }_{{\rm L}}}{{\alpha }_{{\rm R}}}$ given in equation (11) should be equal to zero.

The solutions of equation (16) are given by

$$\begin{eqnarray}&&4{{\Delta }^{2}}+{{\mu }^{2}}{{{\rm sec} }^{2}}\left( \frac{\pi p}{N+1} \right)\;=\;4{{w}^{2}},\end{eqnarray} \tag{ 17 }$$

where p is an integer equal to $1,2,\cdots ,N/2$ if N is even and $1,2,\cdots ,(N-1)/2$ if N is odd. In terms of the variables $\mu /w$ and $\Delta /w$, we observe that equation (17) defines a number of ellipses, which are labeled by the integer p; these are shown in figure 6(a) for N even and figure 6(b) for N odd. Note that all the ellipses pass through the two points given by $\mu =0$ and $\Delta =\pm w$. Figure 6(b) for N odd also contains a zero energy line lying at $\mu =0$ for all values of $\Delta /w$.

Zoom In Zoom Out Reset image size

Equation (17) can be understood in a simple way for the special case $\Delta =0$. Equation (17) then reduces to

$$\begin{eqnarray}&&\mu \;=\;\pm 2w\;{\rm cos} \left( \frac{\pi p}{N+1} \right).\end{eqnarray} \tag{ 18 }$$

We can understand this as follows. For $\Delta =0$, equation (1) describes a non-superconducting tight-binding model whose single particle spectrum is given, for an open chain with N sites, by ${{E}_{q}}=2w{\rm cos} [\pi q/(N+1)]-\mu$, where $q=1,2,\cdots ,N$. One of these energies vanishes whenever μ satisfies the condition given in equation (18), in particular, when p is equal to the smaller of the two integers q and $N+1-q$. This is where the ranges of p mentioned above, namely, $p=1,2,\cdots ,N/2$ for N even and $1,2,\cdots ,(N-1)/2$ for N odd come from. In addition, if N is odd (but not if N even), we have a zero energy state at $\mu =0$ corresponding to $q=(N+1)/2$.

Having found all the zero energy lines in the plane defined by $\mu /w$ and $\Delta /w$, we observe that the parity of the fermion number of the ground state flips sign whenever we cross one of these lines. As a check, this is again easy to see for the case $\Delta =0$. The number of energy levels which are occupied in the ground state changes by 1 and hence the fermion parity changes sign whenever one of the single particle energies E_q given above crosses zero. For $\Delta =0$ and very large negative values of μ, we can see that the ground state of equation (1) contains no fermions; hence, according to equation (13), the fermion parity is +1 for any value of N. For very large positive values of μ, the ground state of equation (1) is completely filled with N fermions; hence the fermion parity according to equation (13) is ${{(-1)}^{N}}$. We remark that the oscillations in the parity, which are related to the Kitaev's Pfaffian, map to spin–spin correlations in the transverse spin chain, and as with much of the literature on Majorana wires, these oscillations have been discussed in depth in the spin context [79].

Previous work involving the dynamics of quenching in the Majorana wire described above has focused on tuning through quantum critical points separating topological and trivial phases [31, 34, 61, 62]. There have been recent works on the effect of quenching on Majorana modes, signatures of Majorana modes in quenching dynamics and Kibble–Zurek scaling [80–84]. (As mentioned earlier, while the terms quench and ramp are sometimes used to distinguish between instantaneous change of parameter versus a time-dependent change at some given rate, respectively, here we use quench in a more general sense to encompass all such dynamic tuning. Our actual studies are restricted to the ramp case.) While our analysis also explores non-equilibrium dynamics within a particular topological phase, we use similar protocols for changing parameters of the Hamiltonian to tune from one phase to another.

Specifically, we consider a linear variation with time of the chemical potential of the system so as to go across the critical line $\mu =2w$,

$$\begin{eqnarray}&&\mu (t)=(2-{{\mu }_{i}})t/\tau +{{\mu }_{i}}.\end{eqnarray} \tag{ 19 }$$

Here ${{\mu }_{i}}$ is the initial chemical potential at t = 0 and $1/\tau$ is the quench rate. Due to the finite rate of variation of μ, the system cannot remain exactly in its ground state and will exhibit non-equilibrium behavior. Namely, excitations (or defects) will be produced in the ground state; this lead to an excess energy of the system and may also lead to a ground state which is in a different topological sector than the initial ground state. These effects can be characterized by the following quantities.

• Defect density: the number of defects produced in the ground state configuration, which is given by the sum over all the excitations.
• Adiabatic Fidelity $\mathcal{O}(t)$: this is the inner product of the instantaneous ground state $|{{\psi }_{{\rm ins}}}(t)\rangle$ of the time-dependent Hamiltonian with the time evolved initial ground state $|\Psi (t)\rangle$,
  $$\begin{eqnarray}&&\mathcal{O}(t)=|\langle \Psi (t)|{{\psi }_{{\rm ins}}}(t)\rangle |.\end{eqnarray} \tag{ 20 }$$
• Residual energy E_res: this is the energy in excess of the instantaneous ground state of the system. We will define this as the dimensionless quantity
  $$\begin{eqnarray}&&{{E}_{{\rm res}}}=[\langle \Psi (t)|H(t)|\Psi (t)\rangle -{{E}_{{\rm G}}}(t)]/|{{E}_{{\rm G}}}(t)|,\end{eqnarray} \tag{ 21 }$$
  where E_G is the energy of the ground state at time t.

Previous work: most of the analytical results for the above quantities obtained in earlier work are in the limit of very large system size or with periodic boundary conditions, where one can Fourier transform the Hamiltonian to momentum space. This in fact reduces the calculation to a well known problem of a transition between two states for each value of the momentum k in the Brillouin zone. This is the famous Landau–Zener–Majorana–Stueckelberg problem [85–88] which can, under a few assumptions, be solved exactly to obtain the probability of excitation from the ground state to the excited state. Using this probability, we can obtain expressions for the defect density, adiabatic fidelity and residual energy. In the limit of long time t, all these quantities have a universal power law scaling as a function of the quench rate which is related to the post-quench excitations. This is the well studied Kibble–Zurek scaling [1–42]. Given the quench rate $1/\tau$, the defect density and residual energy scale as $1/\sqrt{\tau }$ and the adiabatic fidelity scales as ${\rm exp} (c/\sqrt{\tau })$ in this one-dimensional Majorana wire.

For the Majorana wire with open boundary conditions, there have been investigations of the behavior of single particle states under a quench. While it has been found that the single particle bulk states still obey the Kibble–Zurek scaling for the defect density, the quench for an initial state with a Majorana end mode has been found to be non-universal and dependent on the topological features of the system [31]. The end states are not robust with respect to the quench and they delocalize to merge with the bulk states. This leads to a scaling of the defect density as ${{\tau }^{0}}$ (i.e., independent of the quenching rate), which is very different from the Kibble–Zurek scaling.

Another investigation for the open chain looks into a quantity called the Loschmidt echo, which is the survival probability of the Majorana end modes under a quench [61, 80]. Upon quenching across the critical point, the probability decays to extremely small values as the end modes merge with the bulk states when the system is near the critical point where the gap between the end and bulk states vanishes. But interestingly, when the system is quenched to exactly the critical point, although the probability of survival goes to zero initially, it revives itself completely at regular intervals of time. This is attributed to the nearly equal spacing of the low-lying energy levels in the bulk near the critical point; this spacing is of the order of $1/N$ while the scaling of the gap at the critical point is also of order $1/N$. This leads to oscillations in the survival probability with a time period which is proportional to the system size N. We will see below the remnant of this effect in the residual energy for an extremely slow quench.

In this work we will study the quenching dynamics in the above mentioned quantities for an open Majorana wire. We will mainly focus on the many-body states rather than the single particle states and the novel parity switching mechanism discussed in the previous section, which comes about in the topological phase due to the coupling between the Majorana end modes.

In comparison to the translationally invariant systems with periodic boundary conditions studied in previous work, a challenge encountered in these finite-sized systems with open boundary conditions is that one cannot exploit the momentum basis, which in previous works highly simplified the quench dynamical problem. Here,in principle, we are faced with the full-fledged ${{2}^{N}}$ -dimensional Hilbert space associated with fermions on a N-site lattice associated with the Fock space formed by fermion occupancy on each site. Here, we develop and present a dynamic many-body technique to reduce the problem to a numerically tractable form. Our technique hinges on two principles in calculating expectation values or overlaps between states in this time-dependent setting. The first is to use the Heisenberg picture so that the crux of the information on the time evolution is given by the relation between fermionic creation/annihilation operators at different times. The second is to invoke an analog of Wicks theorem for Majorana operators. The computation then reduces to dealing with time-dependent $2N\times 2N$ matrices, allowing us to embark on an exhaustive analysis of adiabatic fidelity and residual energies and to pinpoint attributes of parity blocked dynamics.

Let us start with a general time-dependent Hamiltonian which is quadratic in Majorana operators a_j $(j=1,2,\cdots ,2N)$,

$$\begin{eqnarray}&&H={\rm i}\mathop{\sum }\limits_{i,j=1}^{2N}{{a}_{i}}{{M}_{ij}}(t){{a}_{j}}.\end{eqnarray} \tag{ 22 }$$

Here M(t) is a real antisymmetric matrix; its elements will be functions of w, Δ and μ for the Majorana wire Hamiltonian in equation (2). This can be converted to the canonical form

$$\begin{eqnarray}&&H=4\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}(t)b_{j}^{\dagger }(t){{b}_{j}}(t),\end{eqnarray} \tag{ 23 }$$

up to a constant, by a time-dependent transformation B(t),

$$\begin{eqnarray}&&\bar{b}(t)=B(t)\bar{a}.\end{eqnarray} \tag{ 24 }$$

Here $\bar{b}(t)$ is a $(2N)$-component vector $\bar{b}={{({{b}_{1}},{{b}_{2}},\cdots ,{{b}_{N}},b_{1}^{\dagger },\cdots ,b_{N}^{\dagger })}^{T}}$ and so is $\bar{a}={{({{a}_{1}},{{a}_{2}},\cdots ,{{a}_{2N}})}^{T}}$. The $(2N)$-dimensional matrix B(t) comprises of the eigenvectors of H(t) and it belongs to the group $U(2N)$ with ${\rm det} (B)=\pm 1$. The eigenvalues of the Hamiltonian are $\pm {{\lambda }_{j}}$.

Adiabatic fidelity calculation: as defined in a previous section the adiabatic fidelity is given by $\mathcal{O}(t)=|\langle {{\psi }_{{\rm ins}}}(t)|\Psi (t)\rangle |$. The corresponding annihilation operators of these many-body states satisfy the relations ${{b}_{j}}(t)|{{\psi }_{{\rm ins}}}(t)\rangle =0$ and ${{\beta }_{j}}(t)|\Psi (t)\rangle =0$, where $|\Psi (t)\rangle =S(t,0)\Psi (0)$ and $\bar{\beta }(t)=B(0)S(t,0)\bar{a}$. Here $S(t,0)=\mathcal{T}{\rm exp} (-4\int _{0}^{t}M(t^{\prime} ){\rm d}t^{\prime} )$ is the evolution operator, with $\mathcal{T}$ denoting time ordering. The two sets of annihilation operators are related by

$$\begin{eqnarray}&&\bar{\beta }(t)=B(0)S(t,0){{[B(t)]}^{-1}}\bar{b}(t)=G(t)\bar{b}(t).\end{eqnarray} \tag{ 25 }$$

The key idea underlying the calculation in real space is to express the quantities of interest to us in terms of objects which can be calculated numerically in a simple way. Given the form of the initial Hamiltonian H(0) and the time-dependent H(t), we can see that the quantities B(0), B(t), $S(t,0)$ and G(t) can be easily computed. Given these and the annihilation operators for the ground states, the calculation of the adiabatic fidelity reduces to a computation of the determinant of an antisymmetric matrix A given by

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{A}_{jk}} & = & \langle {{\psi }_{{\rm ins}}}(t)|{{{\bar{\beta }}}_{j}}(t){{{\bar{\beta }}}_{k}}(t)|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\lt k, \\ {} & = & -\langle {{\psi }_{{\rm ins}}}(t)|{{{\bar{\beta }}}_{k}}(t){{{\bar{\beta }}}_{j}}(t)|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\gt k, \\ {} & = & 0\;\;{\rm for}\;\;j=k. \\ \end{array}\end{eqnarray*}$$

We now directly state an important result, deferring the detailed derivation to appendix. The adiabatic fidelity defined in equation (20) is given by

$$\begin{eqnarray}&&\mathcal{O}(t)=|{\rm det} (A){{|}^{1/4}}.\end{eqnarray} \tag{ 26 }$$

Given this relation and the Hamiltonian H(t), we can numerically calculate the adiabatic fidelity as a function of time for a system with open boundary conditions. This can naturally be extended to periodic/antiperiodic boundary conditions as well.

Residual energy calculation: another quantity of interest, the residual energy, defined in equation (21), measures the excess energy contained in the time-evolved quench dependent state compared to the instantaneous ground state energy. This quantity can also be calculated with the real space formalism developed in this section. Following the same strategy as for the adiabatic fidelity, the final expression can simply be expressed in terms of the matrix G(t),

$$\begin{eqnarray}&&{{E}_{{\rm res}}}=[4\mathop{\sum }\limits_{j,k}^{N}{{\lambda }_{j}}G_{N+j,k}^{-1}(t)G_{j,k+N}^{-1}(t)]/|{{E}_{{\rm G}}}(t)|.\end{eqnarray} \tag{ 27 }$$

The details of the derivation are given in appendix.

We now present the results that we obtain for an open Majorana wire with N sites with the Hamiltonian given in equation (2), where μ varies linearly in time as shown in equation (19). Given this specific form of H(t), we numerically calculate all the quantities B(0), B(t) in equation (24), $S(t,0)$, G(t) in equation (25) and finally the adiabatic fidelity $\mathcal{O}(t)$ in equation (26) and the residual energy E_Res in equation (27). In what follows, we provide an in-depth discussion of the novel phenomenon of topological and associated parity blocking as elucidated in section 2.

For periodic boundary conditions, the effect of the number of sites on the fermion parity of the ground state and the consequent phenomenon called topological blocking on the quenching dynamics of the ground state has been discussed in detail in [62]. For an open chain, we saw in section 4 that the Majorana end modes play an important role in determining the parity. For a fixed Δ and N, the parity changes sign as we sweep across the topological phase by varying μ. Here we will explicitly see this parity blocking effect in the evolution of the ground state within the topological phase. The initial parity of the system is set by the value of ${{\mu }_{i}}$ and the number of sites N. As we will see below, the choice of ${{\mu }_{i}}$ can have drastic consequences, especially for an odd number of sites.

7.1. Adiabatic fidelity and parity blocking

Figures 2 and 7 show the numerical results for the adiabatic fidelity $\mathcal{O}(t)$ along with the parity of the instantaneous ground state for an open chain with an even and odd number of sites, respectively. The case of the initial value ${{\mu }_{i}}=0$ for an odd number of sites will be discussed later.

Zoom In Zoom Out Reset image size

We can see from figures 2 and 7 that for both even and odd number of sites, the system starts in a particular fermion parity sector, and as it moves within the topological phase the instantaneous ground state switches parity regularly. On crossing the critical point it can either have opposite parity from the initial state or the same parity, depending on the initial parity sector. On the other hand, as we are dealing with parity conserved systems, the state which is time evolved from the initial ground state continues to have the same fermion parity. Thus the overlap of the time evolved state with the instantaneous ground state plunges to zero at times when the instantaneous parity becomes opposite to the initial parity. We call these parity oscillations, which occur for an open Majorana wire, as the parity blocking effect. The initial ground state is blocked from having any non-zero overlap with the instantaneous ground state for certain values of μ. Finally, on crossing the quantum critical point it becomes zero at all times if the parity is flipped from the initial ground state; this is also a manifestation of parity blocking. Hence, in figure 2, the case of an even number of sites, the parities of the initial and final ground states are the opposite and the system shows parity blocking for the entire topologically trivial phase, while in figure 7, the parities are matched and there is some residual overlap in the trivial phase.

Domain with no parity blocking: as we showed in section 4, the oscillations in the parity do not occur throughout the parameter space corresponding to the topological phase. Namely, there are no oscillations if ${{\mu }^{2}}\gt 4({{w}^{2}}-{{\Delta }^{2}})$. This implies that there ought to be no parity blocking in the adiabatic fidelity. We indeed see this in the numerical results shown in figure 8.

Zoom In Zoom Out Reset image size

Fermion parity degeneracy for an odd number of sites: for an open chain with an odd number of sites, and for states which belong to the odd fermion sector, the ${{\mu }_{i}}=0$ point is special in that it has two degenerate ground states. To see this, let us define an operator

$$\begin{eqnarray}&&C\;=\;{{{\rm i}}^{N(N-1)/2}}\;{{a}_{1}}{{a}_{4}}{{a}_{5}}{{a}_{8}}\cdots ,\end{eqnarray} \tag{ 28 }$$

where the last term on the right-hand side is given by ${{a}_{2N-1}}$ if N is odd and a_2N if N is even. We note that C is both Hermitian and unitary, so that ${{C}^{2}}=I$. Recalling that ${{f}_{n}}=(1/2)({{a}_{2n-1}}+{\rm i}{{a}_{2n}})$ and $f_{n}^{\dagger }=(1/2)({{a}_{2n-1}}-{\rm i}{{a}_{2n}})$, we find that C generates the particle-hole transformation

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} C{{f}_{n}}C & = & {{(-1)}^{n+N-1}}\;f_{n}^{\dagger }, \\ Cf_{n}^{\dagger }C & = & {{(-1)}^{n+N-1}}\;{{f}_{n}}. \\ \end{array}\end{eqnarray} \tag{ 29 }$$

This is a symmetry of the Hamiltonian in equation (1) if $\mu =0$. We now note that the parity and charge-conjugation operators P and C satisfy $PC={{(-1)}^{N}}CP$. Thus P and C anticommute if the number of sites N is odd. Since P and C both commute with H if $\mu =0$, every energy of a system with an odd number of sites will have a two-fold degeneracy with the two eigenstates having opposite fermion parities. (This can be shown as follows. If $|{{\psi }_{+}}\rangle$ is an eigenstate of H and P with eigenvalues E and +1 respectively, the relations $PC=-CP$ and HC = CH imply that $|{{\psi }_{-}}\rangle =C|{{\psi }_{+}}\rangle$ is an eigenstate of H and P with eigenvalues E and −1 respectively.)

Therefore, starting from ${{\mu }_{i}}=0$ means starting with parity states whose degeneracy is not split. The time evolution of an arbitrarily chosen ground state would therefore be that of a linear combination of both parity states. Even though the parity of the instantaneous ground state would keep switching as one sweeps through the topological phase, the overlap would be finite as the time evolved state will be in a superposition of both parities. But the value of the adiabatic fidelity will be smaller than in the parity blocked case as the amplitude is split between the two superposed states. Figure 9 shows the results for this unique case of quenching. We can clearly see in figure 9 that the splitting is zero at time t = 0. This change in the initial condition drastically affects the evolution of the ground state and we do not see complete parity blocking in this case.

Zoom In Zoom Out Reset image size

Comparison with a closed chain: the case of parity blocking in closed chains with periodic and antiperiodic boundary conditions has been studied in detail in a previous work [62]. While the previous analysis involved momentum modes, simplifying the problem to a set of two-level Landau–Zener systems, the real space formalism is easily extended here to compute all the quantities for a closed chain. The numerical results are shown in figure 10. Even though there is parity blocking in the case of an open chain, we may still expect the envelope of the adiabatic fidelity to be comparable with that of the periodic case. In figure 10, we see a good match between the envelopes in the two cases. For a closed chain, the final parity can flip from the initial state depending on the boundary conditions. As can be seen, the parity does not change for the periodic closed chain and therefore one still has a finite overlap. The match between the envelopes in the open and closed chain cases is very close. This suggests that overall the open chain would also respect the Kibble–Zurek behavior for the defect production and excitation density that was found in the closed chain case. The crucial difference between open and closed chains is the parity blocking and switching due to the coupling of the two Majorana end modes.

Zoom In Zoom Out Reset image size

7.2. Residual energy

The numerical results for the variation of the residual energy with time using the full many-body formulation in equation (27) is shown in figure 11 for both open and closed chains. The two cases have the same average behavior. Both show a rapid increase in the energy of the system above the instantaneous many-body ground state as they approach the critical point. This rapid rise is due to the system falling out of equilibrium upon approaching the critical point and thus losing adiabaticity. Far beyond the critical point, we find that the energy asymptotes to a fixed average value. The scaling analysis of this quantity for the transverse field Ising chain has been studied numerically in [89].

Zoom In Zoom Out Reset image size

Effect of Majorana end modes: the crucial difference between the open and closed chains in figure 11 is the presence of the abrupt jumps at small times in the case of the open chain. Reflecting the behavior of adiabatic fidelity, these jumps correspond to switching back and forth between the ground state and excited states due to parity blocking. Upon comparing the behavior of the residual energy with the energy splitting and parity switching of figure 4(b), we find that there is a complete match between the points at which the jumps in the residual energy take place and the points where the parity switching occurs.

Signatures of Loschmidt echoes: in addition to parity blocking, in our many-body system, we find evidence for the Majorana-mode related physics found in previous work on single particle dynamics in quenching [61]. The Loschmidt echo studied in this work calculates the probability for an initial Majorana mode to become a single-particle bulk excitation as a function of time as one sweeps through the critical point. If one varies $\mu (t)$ extremely slowly so as to be close to the adiabatic limit, the gap between the mid-gap end states and bulk states scales as $1/N$ on approaching the critical point since the dynamical critical exponent is equal to 1. The level spacing of the low lying bulk states also scales as $1/N$. Hence the Loschmidt echo turns out to be a periodic function with period N as one quenches to or across the critical point.

In our full-fledged many-body treatment, these echoes appears as 'chirps' of excitations whose occurrence has a period of N. These excitations contribute to the overall energy of the state. Note that this is true only if the quench rate is extremely slow. As shown in figure 12, one can see these oscillations in the numerical results at a low quench rate like $1/\tau =0.1$. The frequency of occurrence of these chirps in the excitations is indeed halved when the system size is doubled. These chirps also appear in the adiabatic fidelity.

Zoom In Zoom Out Reset image size

To summarize, open and closed chains broadly show similar average behaviors as expected for quench dynamics, However, both in the adiabatic fidelity and residual energy, distinct non-analytic features arise in the form of jumps only for an open chain, and these can be attributed directly to the presence of end modes and their associated fermion parity.

In conclusion, our study of non-equilibrium behavior in finite-sized Majorana wires demonstrates that the presence of topological order can dramatically alter quench dynamics. Previous work involving Majorana wires having periodic boundary conditions brought to light the notion of topological blocking in tuning between different quantum phases [62]. In contrast, we have seen here that the coupling between Majorana end modes and the associated ground state parity flips as a function of a tuning parameter gives rise to a more drastic manifestation of topological blocking due to a succession of switches between topological sectors within a single topological phase. As a result, some common measures studied in the quench dynamics literature, such as wave function overlaps between time-evolved states and instantaneous ground states (the adiabatic fidelity), and residual energies, show a series of non-analytic structures in the form of characteristic jumps which are not observed in standard Kibble–Zurek physics.

Our work has shown that there is a much richer texture in the phase diagram of the Kitaev chain or the Majorana wire than has been presented earlier. The circle which separates the regions of oscillating and purely decaying wave functions of the Majorana end modes exists in the thermodynamic limit. In addition, we have shown that a coupling between the end modes, due to a finite length of the chain, leads to further divisions within the circle in the form of ellipses, each division corresponding to a particular fermion parity of the ground state. Although the energy splitting between the Majorana modes goes to zero as we increase the chain length, the number of fermion parity switches increases linearly with the length. This has a dramatic consequence for the adiabatic fidelity under a quench, namely, parity blocking occurs more frequently as we increase the system size. (This is very different from conventional finite size effects which typically vanish in the thermodynamic limit.) We therefore see that parity blocking is not merely a finite size effect, but is a relevant manifestation of the physics of Majorana modes and their topological nature in any real system.

This study shows that quench dynamics serves well as a probe of topological order. While blocking features need not be unique to topological systems in that quantum invariants in other systems can possibly have similar effects, they are necessary conditions under appropriate circumstances (for instance, open boundary condition in the case studied here). Moreover, unlike in most other systems, such as ferromagnets having local order, we expect this blocking phenomenon to be robust against local perturbations. In the case of a finite-sized Majorana wire, the succession of parity switches associated with topological sectors is a crucial aspect of topological order; while studies of the static behavior have been extensive (see, for instance, [90–93]), here we have performed the nearly unexplored study of their effect on quench dynamics.

In fact, the issue of parity forms the basis of several discussions and proposals for Majorana wires, particularly in light of the potential experimental discovery of isolated Majorana end modes and their implications for topological braiding and quantum computing. Several schemes involve changing the on-site chemical potential at specific locations as a means of manipulating and dynamically moving the isolated end modes. A popular study regarding the end modes is the fractional Josephson effect (see, for instance, [58, 94–98]), which involves parity switches between two finite-sized Majorana wires connected to each other at their ends and their effect on Josephson physics (in principle, other zero energy end bound states could mediate such an effect). Our study here is highly relevant to these lines of investigation and it provides a dynamic quantum many-body formulation that goes beyond quasi-static approximations.

While our study primarily aims to understand the effects of parity switching on issues typically studied in the literature on quench dynamics, an experimental setup probing the predictions would be remarkable. While the arena of cold atomic gases is more ideally suited for measurements of residual energy and adiabatic fidelity, realizing topological order in these systems is still in its initial phases [99]. In the setting involving spin–orbit coupled wires, where the isolated end modes have potentially been observed, in principle, ground state parity switches can be observed by coupling the wire to another system. For instance, a possible read-out could involve tunnel-coupling to a quantum dot or STM tip (see, for example, [100, 101]). Further studies would involve pinpointing ways of measuring the behavior of adiabatic fidelity predicted here in such a setup.

Finally, the study of topological/parity blocking in quench dynamics presented here and the associated quantum many-body formulation offer wide scope for further exploration. Several aspects of this initial study require more detailed investigation, for instance, more involved studies of system size, and further connections with Kibble–Zurek scaling and single-particle physics, including anomalous scaling due to boundary effects and appearance of Loschmidt echoes. Oscillations have been found in the derivative of the Renyi entropy with respect to the chemical potential in [102] and it may be interesting to see if there are such effects related to the oscillations in the ground state parity. A host of open issues related to topological blocking in Majorana wires include constraints on thermalization imposed by topological order, effects of external potentials, such as quasiperiodic potentials and disorder, and higher dimensional analogs, such as the Kitaev honeycomb model.

We thank Bryan Clarke, Graham Kells and Jeffrey Teo for illuminating conversations. This work is supported by the National Science Foundation under the grants DMR 0644022-CAR (SH) and DMR-0906521 (VS), the U.S. Department of Energy, Division of Materials Sciences under Award No. DE-FG02-07ER46453 (SV), and by the DST, India, under Project No. SR/S2/JCB-44/2010 (DS). For their hospitality during the course of this work, DS gratefully acknowledges the Department of Physics, University of Illinois at Urbana-Champaign and SV the Aspen Center for Physics.

A.1. Calculation of adiabatic fidelity $\mathcal{O}(t)$

Let us start with a general Hamiltonian which is quadratic in Majorana operators a_j,

$$\begin{eqnarray}&&H={\rm i}\mathop{\sum }\limits_{i,j=1}^{2N}{{a}_{i}}{{M}_{ij}}(t){{a}_{j}}.\end{eqnarray} \tag{ A.1 }$$

The antisymmetric Hermitian matrix iM has real eigenvalues which come in pairs $\pm {{\lambda }_{j}}$ with corresponding eigenvectors x_j and $x_{j}^{*}$, which are orthonormal to each other. We define a set of linear combinations of Majorana operators in terms of these eigenvectors as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} b_{j}^{\dagger }(t) & = & \frac{1}{\sqrt{2}}\mathop{\sum }\limits_{i=1}^{2N}\;{{(x_{j}^{T})}_{i}}{{a}_{i}}, \\ {{b}_{j}}(t) & = & \frac{1}{\sqrt{2}}\mathop{\sum }\limits_{i=1}^{2N}\;{{(x_{j}^{\dagger })}_{i}}{{a}_{i}}. \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

In terms of a $(2N)$-component vector $\bar{b}={{({{b}_{1}},{{b}_{2}},\cdots ,{{b}_{N}},b_{1}^{\dagger },\cdots ,b_{N}^{\dagger })}^{T}}$, equation (A.2) can be expressed as a linear transformation

$$\begin{eqnarray}&&\bar{b}(t)=B(t)\bar{a}.\end{eqnarray} \tag{ A.3 }$$

The rows of the $(2N)$-dimensional matrix B are the eigenvectors x^T_j and $x_{j}^{\dagger }$, and B belongs to the unitary group $U(2N)$ with ${\rm det} (B)=\pm 1$. In terms of b_j, the Hamiltonian (A1) becomes, up to a constant

$$\begin{eqnarray}&&H=4\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}b_{j}^{\dagger }{{b}_{j}}.\end{eqnarray} \tag{ A.4 }$$

If $|\Psi (0)\rangle$ is the initial ground state of H at t = 0, then ${{b}_{j}}(0)|\Psi (0)\rangle =0$. The instantaneous ground state $|{{\psi }_{{\rm ins}}}(t)\rangle$ of H(t) is annihilated by b_j(t), namely, ${{b}_{j}}(t)|{{\psi }_{{\rm ins}}}(t)\rangle =0$.

We now want to find the adiabatic fidelity $\langle \Psi (t)|{{\psi }_{{\rm ins}}}(t)\rangle$. Let us examine the operators which annihilate the time evolved state $|\Psi (t)\rangle$ and some of their properties. Let us say ${{\beta }_{j}}(t)|\Psi (t)\rangle =0$. We want to put the information of the time evolution with H(t) into ${{\beta }_{j}}(t)$. The time evolution of the Majorana operators a_j is given by their Heisenberg equations of motion:

$$\begin{eqnarray}&&\frac{{\rm d}{{a}_{j}}(t)}{{\rm d}t}=-{\rm i}[H(t),{{a}_{j}}(t)]=-4\mathop{\sum }\limits_{k=1}^{2N}{{M}_{jk}}(t){{a}_{k}}(t).\end{eqnarray} \tag{ A.5 }$$

In terms of a $(2N)$-component vector $\bar{a}$, the solution of the equation ${\rm d}\bar{a}(t)/{\rm d}t=-4\;M(t)\bar{a}(t)$ is given by the evolution operator

$$\begin{eqnarray}&&\bar{a}(t)=S(t,0)\bar{a}(0).\end{eqnarray} \tag{ A.6 }$$

Here $S(t,0)=\mathcal{T}{\rm exp} (-4\int _{0}^{t}M(t^{\prime} ){\rm d}t^{\prime} )$ which can be calculated numerically for a given M(t). It can now be shown that $\bar{\beta }(t)=B(0)\bar{a}(t)=B(0)S(t,0)\bar{a}(0)$, where $\bar{\beta }(t)$ is the $(2N)$-component vector comprising of ${{\beta }_{j}}(t),\beta _{j}^{\dagger }(t)$ similar to $\bar{b}(t)$. Therefore the relation between $\bar{\beta }(t)$ and $\overline{b(t)}$ is given by

$$\begin{eqnarray}&&\bar{\beta }(t)=B(0)S(t,0){{[B(t)]}^{-1}}\bar{b}(t).\end{eqnarray} \tag{ A.7 }$$

The key idea underlying the calculation in real space is to express all the quantities of interest in terms of those which can be calculated numerically. We can see that given the form of the initial Hamiltonian H(0) and the time-dependent H(t), the quantities B(0), B(t) and $S(t,0)$ can be easily computed. Given these and the annihilation operators of the ground states, we will now derive the final expression for the adiabatic fidelity.

Consider the Fock space of ${{2}^{N}}$ states $|{{\phi }_{a}}\rangle$, $a=1,2,\cdots ,{{2}^{N}}$. Their fermionic occupation numbers are given by $\beta _{j}^{\dagger }(t){{\beta }_{j}}(t)=0$ or 1. For a $|\Psi (t)\rangle$ belonging to this Fock space, ${{\beta }_{j}}(t)|\Psi (t)\rangle =0$; hence $\beta _{j}^{\dagger }(t){{\beta }_{j}}(t)|\Psi (t)\rangle =0$ for all j. We now define the following operators L_j

$$\begin{eqnarray}&&{{L}_{j}}={{\beta }_{j}}(t)\;\;{\rm for}\;\;1\leqslant j\leqslant N,\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&\qquad =\beta _{2N+1-j}^{\dagger }\;\;{\rm for}\;\;N+1\leqslant j\leqslant 2N.\end{eqnarray} \tag{ A.9 }$$

In the Fock space, $|\Psi (t)\rangle$ is the only state which is not annihilated by the product ${{L}_{1}}{{L}_{2}}\;\cdots \;{{L}_{2N}}$. In fact,

$$\begin{eqnarray}&&{{L}_{1}}{{L}_{2}}\cdots {{L}_{2N}}|\Psi (t)\rangle =|\Psi (t)\rangle .\end{eqnarray} \tag{ A.10 }$$

We will make use of this fact to reduce the calculation in the ${{2}^{N}}$ dimensional Fock space to a matrix computation in $2N$ dimensions as follows. Consider the quantity

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \langle {{\psi }_{{\rm ins}}}(t)|{{L}_{1}}{{L}_{2}}\;\cdots \;{{L}_{2N}}|{{\psi }_{{\rm ins}}}(t)\rangle \\ \quad =\langle {{\psi }_{{\rm ins}}}(t)|{{L}_{1}}{{L}_{2}}\;\cdots \;{{L}_{2N}}\mathop{\sum }\limits_{a=1}^{{{2}^{N}}}|{{\phi }_{a}}\rangle \langle {{\phi }_{a}}|{{\psi }_{{\rm ins}}}(t)\rangle \\ \quad =\langle {{\psi }_{{\rm ins}}}(t)|{{L}_{1}}{{L}_{2}}\;\cdots \;{{L}_{2N}}|\Psi (t)\rangle \langle \Psi (t)|{{\psi }_{{\rm ins}}}(t)\rangle \\ \quad =|\langle {{\psi }_{{\rm ins}}}(t)|\Psi (t)\rangle {{|}^{2}}. \\ \end{array}\end{eqnarray} \tag{ A.11 }$$

This can be simplified further using a form of Wick's theorem given in [103], in terms of a $(2N)$-dimensional antisymmetric matrix A defined as

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{A}_{jk}} & = & \langle {{\psi }_{{\rm ins}}}(t)|{{L}_{j}}{{L}_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\lt k, \\ {} & = & -\langle {{\psi }_{{\rm ins}}}(t)|{{L}_{k}}{{L}_{j}}|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\gt k, \\ {} & = & 0\;\;{\rm for}\;\;j=k. \\ \end{array}\end{eqnarray*}$$

This matrix can be calculated using the relation between $\bar{\beta }(t)$ and $\bar{b}(t)$ and the fact that ${{b}_{j}}(t)|{{\psi }_{{\rm ins}}}\rangle (t)=0$. Finally we get

$$\begin{eqnarray*}&&\langle {{\psi }_{{\rm ins}}}(t)|{{L}_{1}}{{L}_{2}}\cdots {{L}_{2N}}|{{\psi }_{{\rm ins}}}(t)\rangle =Pf(A)=|\langle {{\psi }_{{\rm ins}}}(t)|\Psi (t)\rangle {{|}^{2}}.\end{eqnarray*}$$

Since the Pfaffian is given by $Pf(A)=\pm \sqrt{{\rm det} (A)}$, we see that the adiabatic fidelity is

$$\begin{eqnarray}&&|\langle {{\psi }_{{\rm ins}}}(t)|\Psi (t)\rangle |=|{\rm det} (A){{|}^{1/4}}.\end{eqnarray} \tag{ A.12 }$$

Now we need to calculate the matrix elements A_jk. Since A is antisymmetric, we need to calculate only the elements for $j\lt k$. These are given by

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {} & {} & {{A}_{jk}} \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|{{\beta }_{j}}{{\beta }_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\leqslant N,k\leqslant N, \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|{{\beta }_{j}}\beta _{2N+1-k}^{\dagger }|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\leqslant N,k\gt N, \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|\beta _{2N+1-j}^{\dagger }\beta _{2N+1-k}^{\dagger }|{{\psi }_{{\rm ins}}}(t)\rangle \;\;{\rm for}\;\;j\leqslant N,k\gt N. \\ \end{array}\end{eqnarray*}$$

We need to evaluate each of these terms We introduce $G(t)=B(0)S(t,0){{B}^{-1}}(t)$ so that $\beta (t)=G(t)b(t)$. Consider the first case $j\leqslant N,k\leqslant N$:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {} & {} & \langle {{\psi }_{{\rm ins}}}(t)|{{\beta }_{j}}{{\beta }_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle \\ {} & = & \mathop{\sum }\limits_{m,n}\langle {{\psi }_{{\rm ins}}}(t)|{{G}_{jm}}{{{\bar{b}}}_{m}}{{G}_{kn}}{{{\bar{b}}}_{n}}|{{\psi }_{{\rm ins}}}(t)\rangle \\ {} & = & \mathop{\sum }\limits_{m,n}{{G}_{jm}}{{D}_{mn}}{{G}_{kn}}. \\ \end{array}\end{eqnarray} \tag{ A.13 }$$

The matrix D is given by

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {} & {} & {{D}_{jk}}\;=\;\langle {{\psi }_{{\rm ins}}}(t)|{{{\bar{b}}}_{j}}{{{\bar{b}}}_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|{{b}_{j}}{{b}_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle =0\;\;{\rm for}\;\;j\leqslant N,k\leqslant N, \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|{{b}_{j}}b_{k-N}^{\dagger }|{{\psi }_{{\rm ins}}}(t)\rangle ={{\delta }_{j,k-N}}\;\;{\rm for}\;\;j\leqslant N,k\gt N, \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|b_{j-N}^{\dagger }{{b}_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle =0\;\;{\rm for}\;\;j\gt N,k\leqslant N, \\ {} & = & \langle {{\psi }_{{\rm ins}}}(t)|b_{j-N}^{\dagger }b_{k-N}^{\dagger }|{{\psi }_{{\rm ins}}}(t)\rangle =0\;\;{\rm for}\;\;j\gt N,k\gt N. \\ \end{array}\end{eqnarray*}$$

Using this fact, we get

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {} & {} & \langle {{\psi }_{{\rm ins}}}(t)|{{\beta }_{j}}{{\beta }_{k}}|{{\psi }_{{\rm ins}}}(t)\rangle \\ {} & = & \mathop{\sum }\limits_{m\leqslant N}{{G}_{jm}}(t){{G}_{k,m+N}}(t) \\ {} & = & ({\rm first}\;{\rm half}\;{\rm of}\;j{\rm th}\;{\rm row}\;{\rm of}\;{\rm G})\times {{({\rm second}\;{\rm half}\;{\rm of}\;k{\rm th}\;{\rm row}\;{\rm of}G)}^{T}}. \\ \end{array}\end{eqnarray} \tag{ A.14 }$$

Similarly the other elements of A_jk can be found and numerically evaluated as a function of time. Once we have the matrix A, the adiabatic fidelity is simply related to its determinant.

A.2. Parity in a two-site problem

For an open chain, we saw in section 4 that the overall parity of the ground state is decided by the fermion parity of the split energy states arising from the overlap of the Majorana end modes. We will consider here the effective Hamiltonian for such a system and illustrate how det(B) determines the parity of the system. The effective Hamiltonian for the coupled Majoranas is given by

$$\begin{eqnarray}&&{{H}_{f}}={\rm i}2\;J{{a}_{1}}{{a}_{2N}}.\end{eqnarray} \tag{ A.15 }$$

The eigenvalues are given by $\lambda =\pm J$ and the eigenvectors are x and its conjugate x*, where

$$\begin{eqnarray}&&x=\frac{1}{\sqrt{2}}\left( \frac{1}{-{\rm i}} \right).\end{eqnarray} \tag{ A.16 }$$

Using these eigenvectors we can construct the matrix B which transforms the Hamiltonian into the canonical form.

$$\begin{eqnarray*}B=\frac{1}{2}\left[ \begin{array}{ccccccccccccccc} 1 & {\rm i} \\ 1 & -{\rm i} \\ \end{array} \right]\end{eqnarray*}$$

Now suppose that J changes sign; then the eigenvalues are flipped and the rows in B are also flipped. The sign of the determinant of a matrix changes when two of its rows are interchanged. From equation (3) we see that as the value of J changes sign, the energy corresponding to the state with a particular fermion parity changes. The sign of det(B) precisely tracks this flip in the parity of the energy level which contributes to the ground state.

One can also see this from the calculation of the Pfaffian of the Hamiltonian for the two-site problem. The Hamiltonian in the Majorana basis is given by:

$$\begin{eqnarray*}\left( \begin{array}{ccccccccccccccc} 0 & -{\rm i}\mu /2 & 0 & {\rm i}(-w+\Delta )/2 \\ {\rm i}\mu /2 & 0 & {\rm i}(w+\Delta )/2 & 0 \\ 0 & -{\rm i}(w+\Delta )/2 & 0 & -{\rm i}\mu /2 \\ -{\rm i}(-w+\Delta )/2 & 0 & {\rm i}\mu /2 & 0 \\ \end{array} \right)\end{eqnarray*}$$

The Pfaffian of this Hamiltonian is given by

$$\begin{eqnarray}&&Pf(H)=\frac{{{\mu }^{2}}}{4}-\frac{{{w}^{2}}-{{\Delta }^{2}}}{4}.\end{eqnarray} \tag{ A.17 }$$

From this expression the Pfaffian changes its sign at $\frac{{{\mu }^{2}}}{4}=\frac{{{w}^{2}}-{{\Delta }^{2}}}{4}$. One can see that this is precisely the condition we have in equation (17) for N = 2. As shown in Kitaev's paper [58] we have the condition:

$$\begin{eqnarray}&&P(H)={\rm sgn} [Pf(H)]={\rm sgn} [{\rm det} (B)].\end{eqnarray} \tag{ A.18 }$$

A.3. Calculation of residual energy

The residual energy is defined as

$$\begin{eqnarray}&&{{E}_{{\rm res}}}=[\langle \Psi (t)|H(t)|\Psi (t)\rangle -{{E}_{{\rm G}}}(t)]/|{{E}_{{\rm G}}}(t)|,\end{eqnarray} \tag{ A.19 }$$

where E_G(t) is the instantaneous ground state energy. Let us calculate the first term in the expression, $\langle H(t)\rangle$. The matrix B(t) transforms the time-dependent Hamiltonian to the canonical form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} H(t) & = & 4\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}(t)b_{j}^{\dagger }(t){{b}_{j}}(t)-2\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}(t) \\ {} & = & 4\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}(t){{{\bar{b}}}_{N+j}}(t){{{\bar{b}}}_{j}}(t)-2\mathop{\sum }\limits_{j=1}^{N}{{\lambda }_{j}}(t). \\ \end{array}\end{eqnarray} \tag{ A.20 }$$

Therefore, $\langle H(t)\rangle =4\sum _{j\;\,=\;\,1}^{N}{{\lambda }_{j}}\langle \Psi (t)|{{\bar{b}}_{N\;\,+\;\,j}}(t){{\bar{b}}_{j}}(t)|\Psi (t)\rangle -2\sum _{j\;\,=\;\,1}^{N}{{\lambda }_{j}}(t)$. From equation (A.7) and using the definition $G(t)=B(0)S(t,0)B{{(t)}^{-1}}$ from the last section, we obtain

$$\begin{eqnarray}&&{{\bar{b}}_{i}}(t)=\mathop{\sum }\limits_{j}G_{ij}^{-1}(t){{\bar{\beta }}_{j}}(t).\end{eqnarray} \tag{ A.21 }$$

As before, let us try to reduce everything to quantities which can be numerically computed

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {} & {} & \langle H(t)\rangle +2\mathop{\sum }\limits_{j}^{N}{{\lambda }_{j}}(t) \\ {} & = & 4\mathop{\sum }\limits_{j}^{N}{{\lambda }_{j}}(t)\langle \Psi (t)|{{{\bar{b}}}_{N+j}}(t){{{\bar{b}}}_{j}}(t)|\Psi (t)\rangle \\ {} & = & 4\mathop{\sum }\limits_{j,k,l}^{N}{{\lambda }_{j}}(t)\langle \Psi (t)|G_{N+j,k}^{-1}(t){{{\bar{\beta }}}_{k}}(t)G_{j,l}^{-1}(t){{{\bar{\beta }}}_{{\rm L}}}(t)|\Psi (t)\rangle \\ {} & = & 4\mathop{\sum }\limits_{j,k,l}^{N}{{\lambda }_{j}}(t)G_{N+j,k}^{-1}(t)\langle \Psi (t)|{{{\bar{\beta }}}_{k}}(t){{{\bar{\beta }}}_{{\rm L}}}(t)|\Psi (t)\rangle G_{j,l}^{-1}(t). \\ \end{array}\end{eqnarray*}$$

Now using

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{{\bar{\beta }}}_{j}} & = & {{\beta }_{j}}\;\;{\rm for}\;\;j\leqslant N, \\ {{{\bar{\beta }}}_{j}} & = & \beta _{j}^{\dagger }\;\;{\rm for}\;\;N\lt j\leqslant 2N, \\ \end{array}\end{eqnarray*}$$

we have

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {} & {} & \langle \Psi (t)|{{{\bar{\beta }}}_{k}}(t){{{\bar{\beta }}}_{{\rm L}}}(t)|\Psi (t)\rangle \\ {} & = & \langle \Psi (t)|{{\beta }_{k}}(t){{\beta }_{l}}(t)|\Psi \rangle (t)=0\;\;{\rm for}\;\;k\leqslant N,l\leqslant N \\ {} & = & \langle \Psi (t)|{{\beta }_{k}}(t)\beta _{l-N}^{\dagger }(t)|\Psi (t)\rangle ={{\delta }_{k,l-N}}\;\;{\rm for}\;\;k\leqslant N,l\gt N, \\ {} & = & \langle \Psi (t)|\beta _{k-N}^{\dagger }(t){{\beta }_{{\rm L}}}(t)|\Psi (t)\rangle =0\;\;{\rm for}\;\;k\gt N,l\leqslant N, \\ {} & = & \langle \Psi (t)|\beta _{k-N}^{\dagger }(t)\beta _{l-N}^{\dagger }(t)|\Psi (t)\rangle =0\;\;{\rm for}\;\;k\leqslant N,l\leqslant N. \\ \end{array}\end{eqnarray*}$$

Hence the required expression simplifies to

$$\begin{eqnarray}&&\langle H(t)\rangle =4\mathop{\sum }\limits_{j,k}^{N}{{\lambda }_{j}}(t)G_{N+j,k}^{-1}(t)G_{j,k+N}^{-1}(t)-2\mathop{\sum }\limits_{j}^{N}{{\lambda }_{j}}(t).\end{eqnarray} \tag{ A.22 }$$

Finally the expression for the residual energy is given by

$$\begin{eqnarray}&&{{E}_{{\rm res}}}=[4\mathop{\sum }\limits_{j,k}^{N}{{\lambda }_{j}}(t)G_{N+j,k}^{-1}(t)G_{j,k+N}^{-1}(t)]/|{{E}_{{\rm G}}}(t)|.\end{eqnarray} \tag{ A.23 }$$

E_G is calculated simply by summing over all the negative eigenvalues of the numerically diagonalized Hamiltonian, namely, $-\sum _{j=1}^{N}{{\lambda }_{j}}(t)$, and the first term is calculated from G, which we already have to calculate numerically to find the adiabatic fidelity.