Dynamical signatures of topological order in the driven-dissipative Kitaev chain

TL;DR

This paper introduces a method to identify dynamical signatures of topological order in open quantum systems, specifically in a driven-dissipative Kitaev chain, by analyzing bulk winding numbers and edge modes.

Contribution

It proposes a novel approach to classify topological phases in dissipative systems through dynamical signatures and links these to observable edge modes and currents.

Findings

Bulk winding numbers relate to Majorana edge modes.

Dissipation causes decay but preserves dynamical topological signatures.

Periodic driving reveals higher winding numbers and new edge modes.

Abstract

We investigate the effects of dissipation and driving on topological order in superconducting nanowires. Rather than studying the non-equilibrium steady state, we propose a method to classify and detect dynamical signatures of topological order in open quantum systems. Bulk winding numbers for the Lindblad generator $\hat{\mathcal{L}}$ of the dissipative Kitaev chain are found to be linked to the presence of Majorana edge master modes -- localized eigenmodes of $\hat{\mathcal{L}}$. Despite decaying in time, these modes provide dynamical fingerprints of the topological phases of the closed system, which are now separated by intermediate regions where winding numbers are ill-defined and the bulk-boundary correspondence breaks down. Combining these techniques with the Floquet formalism reveals higher winding numbers and different types of edge modes under periodic driving. Finally, we link the presence of edge modes to a steady state current.

Figures

Equations

$$H=-\sum _j\left(J{c}_j^{†}{c}_{j+1}+\gamma c_j^{†}{c}_{j+1}^{†}+\text{H.c.}\right)-\mu \sum _j{c}_j^{†}{c}_j,$$

$${ϵ}_k=\sqrt{4{\gamma }^2{\sin }^2(k)+(2J\cos (k)+\mu {)}^2}.$$

$$w_{2j-1}=c_j+c_j^{†},w_{2j}=i(c_j-c_j^{†}),$$

$$H=-\frac{i}{2}\sum _j\left((\gamma +J)w_{2j}{w}_{2j+1}+(\gamma -J)w_{2j-1}{w}_{2j+2}\right)-\frac{\mu }{2}\sum _j(1+i{w}_{2j}{w}_{2j-1}).$$

$$\nu =\frac{i}{\pi }{\int }_{-\pi }^{\pi }⟨{\psi }_k|\frac{\partial }{\partial k}|{\psi }_k⟩\mathrm{d}k\mathrm{mod}2=\{\begin{array}{ll}0,& |\mu |>2J\\ 1,& |\mu |<2J\end{array},$$

$$\hat{\mathcal{L}}\rho =\frac{\mathrm{d}\rho }{\mathrm{d}t}=-i[H,\rho ]+\sum _m\left(L_m\rho L_m^{†}-\frac{1}{2}{L}_m^{†}{L}_m\rho -\frac{1}{2}\rho L_m^{†}{L}_m\right),$$

$$\hat{\mathcal{L}}=\frac{1}{2}\sum _{i,j}\left(\begin{array}{cc}{\hat{c}}_i^{†}& {\hat{c}}_i\end{array}\right){\mathbf{A}}_{ij}\left(\begin{array}{c}{\hat{c}}_j\\ {\hat{c}}_j^{†}\end{array}\right)-A_0\widehat{\text{mathds}1},$$

$${\hat{c}}_j^{†}\rho =\frac{1}{2}\left(w_j\rho +(\hat{P}\rho )w_j\right),{\hat{c}}_j\rho =\frac{1}{2}\left(w_j\rho -(\hat{P}\rho )w_j\right),$$

$$\mathbf{A}=\left(\begin{array}{cc}-{\mathbf{X}}^{†}& -i\mathbf{Y}\\ 0& \mathbf{X}\end{array}\right),$$

$$\mathbf{X}\equiv -4i\mathbf{H}+2\left(\mathbf{M}+{\mathbf{M}}^T\right),\mathbf{Y}\equiv -4i\left(\mathbf{M}-{\mathbf{M}}^T\right),$$

$$H$$

$$L_m$$

$${\mathbf{X}}^{†}\mathbf{C}+\mathbf{CX}=i\mathbf{Y},C_{ij}\equiv \mathrm{Tr}(w_i{w}_j{\rho }_{ss})-{\delta }_{i,j}.$$

$$\hat{\mathcal{L}}=-\sum _j{\beta }_j{\hat{b}}_j^{\prime }{\hat{b}}_j.$$

$$H=\sum _{j=1}^N\sum _{\ell =-N/2}^{N/2}\left(\begin{array}{cc}{w}_{2j-1}& w_{2j}\end{array}\right){\mathbf{h}}_{\ell }\left(\begin{array}{c}{w}_{2(j+\ell )-1}\\ w_{2(j+\ell )}\end{array}\right),L_j=\sum _{\ell =-N/2}^{N/2}{\underline{l}}_{\ell }\cdot \left(\begin{array}{c}{w}_{2(j+\ell )-1}\\ w_{2(j+\ell )}\end{array}\right),$$

$$\mathbf{h}(k)=\sum _{\ell =-N/2}^{N/2}{e}^{-ik\ell }{\mathbf{h}}_{\ell },\mathbf{m}(k)=\underline{l}(k)\otimes {\underline{l}}^{\ast }(k),\underline{l}(k)=\sum _{\ell =-N/2}^{N/2}{e}^{-ik\ell }{\underline{l}}_{\ell },$$

$$\mathbf{h}(k)=\frac{1}{2}\left(\begin{array}{cc}0& i(J\cos k+\frac{\mu }{2})-\gamma \sin k\\ -i(J\cos k+\frac{\mu }{2})-\gamma \sin k& 0\end{array}\right).$$

$$\mathbf{m}(k)=\frac{g}{4}\left(\begin{array}{cc}(1+\Delta {)}^2& i(1-{\Delta }^2)\\ -i(1-{\Delta }^2)& (1-\Delta {)}^2\end{array}\right).$$

$$\mathbf{x}(k)=$$

$$\mathbf{y}(k)=$$

$$\beta (k)=g(1+{\Delta }^2)\pm \sqrt{4g^2{\Delta }^2-4{\gamma }^2{\sin }^{2}k-(2J\cos k+\mu {)}^2}.$$

$$\mathcal{S}(\mathbf{x}(k)-a_0\text{mathds}1)\mathcal{S}=-(\mathbf{x}(k)-a_0\text{mathds}1),$$

$$\mathbf{x}(k){\underline{u}}_i(k)={\beta }_i(k){\underline{u}}_i(k),{\mathbf{x}}^{†}(k){\underline{v}}_i(k)={\beta }_i^{\ast }(k){\underline{v}}_i(k),$$

$${\underline{v}}_i^{\ast }(k)\cdot {\underline{u}}_j(k)={\delta }_{ij}.$$

$$\displaystyle\mathbf{x}(k)\mathcal{S}\underline{u}_{i}(k)=\big{(}2a_{0}-\beta_{i}^{*}(k)\big{)}\mathcal{S}\underline{u}_{i}(k)\,.$$

$${\nu }_i$$

$$\mathcal{S}{\underline{u}}_i(k)={\sigma }_x{\underline{u}}_i^{\ast }(k)=e^{i{\phi }_i(k)}{\underline{u}}_i(k),\mathcal{S}{\underline{v}}_i(k)={\sigma }_x{\underline{v}}_i^{\ast }(k)=e^{i{\phi }_i(k)}{\underline{v}}_i(k),$$

$${\nu }_i$$

$${\hat{b}}_i^{\prime }={\underline{v}}_i\cdot {\underline{\hat{c}}}^{†},{\hat{b}}_i={\underline{u}}_i\cdot \underline{\hat{c}}-\mathbf{C}{\underline{u}}_i\cdot {\underline{\hat{c}}}^{†}.$$

$$⟨\mathcal{O}(t)⟩=\mathrm{Tr}(d_0^{†}{d}_0\rho (t)),d_0^{†}\equiv \frac{1}{2}(w_1+i{w}_{2N}).$$

$$⟨\mathcal{O}(t)⟩=e^{-2g(1+{\Delta }^2)t}.$$

$$r_L=\frac{⟨c_1^{†}{c}_1⟩-⟨c_{N/2}^{†}{c}_{N/2}⟩}{⟨c_{N/2}^{†}{c}_{N/2}⟩},r_R=\frac{⟨c_N^{†}{c}_N⟩-⟨c_{N/2}^{†}{c}_{N/2}⟩}{⟨c_{N/2}^{†}{c}_{N/2}⟩}.$$

Full text

Dynamical signatures of topological order in the driven-dissipative Kitaev chain

M. T. van Caspel1, S. E. Tapias Arze1, I. Pérez Castillo2,3*

1 Institute of Physics and Delta Institute for Theoretical Physics,

University of Amsterdam, Amsterdam, The Netherlands

2 Department of Quantum Physics and Photonics, Institute of Physics,

National Autonomous University of Mexico, Mexico

3 London Mathematical Laboratory, 8 Margravine Gardens,

W6 8RH London, United Kingdom

• isaacpc@fisica.unam.mx

Abstract

**We investigate the effects of dissipation and driving on topological order in superconducting nanowires. Rather than studying the non-equilibrium steady state, we propose a method to classify and detect dynamical signatures of topological order in open quantum systems. Bulk winding numbers for the Lindblad generator $\hat{\mathcal{L}}$ of the dissipative Kitaev chain are found to be linked to the presence of Majorana edge master modes – localized eigenmodes of $\hat{\mathcal{L}}$. Despite decaying in time, these modes provide dynamical fingerprints of the topological phases of the closed system, which are now separated by intermediate regions where winding numbers are ill-defined and the bulk-boundary correspondence breaks down. Combining these techniques with the Floquet formalism reveals higher winding numbers and different types of edge modes under periodic driving. Finally, we link the presence of edge modes to a steady state current. **

Contents

1. 1 Introduction

2. 2 The Kitaev chain with bulk dissipation

3. 2.1 Free Lindbladian time evolution

4. 3 Topological order and edge modes

5. 3.1 Topological properties of $\mathbf{x}(k)$

6. 3.2 Bulk-boundary correspondence and its breakdown

7. 3.3 Physical interpretation of edge master modes

8. 3.4 Steady state expectation values

9. 4 Higher winding numbers in the driven-dissipative Kitaev chain

10. 4.1 Topology and the Floquet formalism

11. 4.2 Lindblad-Floquet theory

12. 4.3 Numerical results

13. 5 Conclusions and outlook

14. A Geometric interpretation of the winding number

15. B Computing the winding number

16. C Topology of X vs. A

17. D Edge modes of $\mathbf{X}$

18. E Time evolution of observables

19. F Off-diagonal Floquet blocks and their importance

1 Introduction

Since its discovery, topological order is a subject that has attracted significant interest in condensed matter physics, with a number of research directions showing recent exciting activity [1, 2]. One of those is in the field of periodically driven systems, where research into Floquet topological insulators continues to reveal surprising properties [3, 4, 5]. Although these systems are inherently out-of-equilibrium, many defining features, such as long-range order and the bulk-boundary correspondence, remain present. Furthermore, entirely new physics emerges due to the periodic time evolution that results into additional topological invariants and phases.

Floquet theory has proven to be very effective in describing periodically driven quantum systems, which are used to engineer effective couplings leading to new states of matter [6]. However, for closed systems, particularly interacting ones, these periodic protocols often predict an unbounded heating, where they are driven to infinite temperature states at long times. More realistic models should therefore add some particle losses – inevitable in an experimental setting – which compete with the driving and lead the system to a Floquet steady state at finite temperature. Combining the Floquet formalism with the theory of open quantum systems results in an elegant framework for such driven-dissipative models [7, 8]. Unfortunately, this poses a major challenge for the study of topological order, typically defined for the system’s ground state(s), as the non-unitary time evolution produces mixed quantum states. The topological classification of mixed states is an ongoing and contentious effort, with a number of open problems and conflicting findings [9, 10, 11, 12].

Rather than focusing on topological properties of the steady state density matrix, the present work analyzes the full time evolution of an open quantum system. We draw inspiration from the study of topology in non-Hermitian Hamiltonians, which also lack a ground state. A flurry of recent work in this field shows that topological order can be inferred from their complex spectrum [13, 14, 15, 16, 17, 18, 19]. By writing the time evolution of a non-interacting open system in terms of non-Hermitian matrices, we apply some of these techniques and investigate which aspects of the topological order remain when a system is coupled to a bath in two cases: firstly when the Hamiltonian part of the evolution is time independent, and secondly when it is subject to periodic driving.

As a case study, we consider the Kitaev chain, a spinless approximation of a superconducting nanowire with a topological phase that exhibits unpaired Majorana edge states [20]. The natural description in terms of Majorana fermions is well-suited for our treatment of the dissipation, which exploits the properties of Clifford algebras [21, 22].

The structure of this paper is as follows. In Section 2, we describe the model in detail and summarize its main equilibrium properties. In Section 3, we turn to our main goal of describing what happens to the Majorana edge modes and the bulk-boundary correspondence in the dissipative Kitaev chain. For the type of dissipation that we consider – a single, identical channel coupled to each site, with a parameter $\Delta$ interpolating between loss and gain – we find that the topological order of the Kitaev chain is preserved in a dynamical sense: edge modes now decay exponentially over time, but their existence is still guaranteed by symmetry. In addition, we find that the different topological phases are now separated by an intermediate region, where the band gap is closed and exceptional points appear [13]. This is the result of a spontaneously broken anti-unitary symmetry, similar to $\mathcal{PT}$ symmetry in non-Hermitian systems [23].

Finally, in Section 4, periodic driving is added to the system. In this case, the behavior of the system becomes even more complex as a consequence of Floquet resonances [24]. For example, one can now have a seemingly unlimited number of edge modes, which are split into two types [25, 26]. Moreover, in some regions of the parameter space, the edge modes can repel each other and break symmetries in a complex way. Instead of trying to describe and understand the full array of phenomena that can be seen in the driven-dissipative Kitaev chain, our goal is to demonstrate the promise of this approach and the richness of this kind of systems. We conclude in Section 5 with a number of open questions and possible future applications of our work.

2 The Kitaev chain with bulk dissipation

Originally studied in this context in [20], the Kitaev chain is a mean-field model of a 1D $p$-wave superconductor with spinless fermions, given by the Hamiltonian

[TABLE]

where $J$ is the hopping amplitude, $\gamma$ the $p$-wave pairing and $\mu$ the chemical potential. Under periodic boundary conditions (PBC), the system has dispersion relation

[TABLE]

The band gap closes at the critical point $\mu=\pm 2J$ and Kitaev showed that this corresponds to a topological phase transition. The easiest way of seeing this is by mapping to Majorana fermions:

[TABLE]

resulting into the following Hamiltonian

[TABLE]

The appearance of topological edge modes in the non-trivial phase, when considering open boundary conditions (OBC) with $N$ sites, becomes clear at the point $J=\gamma>0,~{}\mu=0$. The Majorana fermions then form uncoupled dimers across neighboring sites, exactly like the limiting case of the prototypical Su-Schrieffer-Heger (SSH) model. At the edges, the modes $w_{1}$ and $w_{2N}$ disappear from the Hamiltonian entirely and therefore commute with $H$. Together they can be interpreted as a single delocalized fermionic mode with zero energy, which is protected from perturbations by symmetry as long as the gap does not close, resulting in a two-fold degeneracy of the ground state.

The symmetry of the system is crucial for topological order in 1D[27]. In the case of the Kitaev chain there is a combination of particle-hole and spatial inversion symmetry, which is often referred to as chiral symmetry. Spinless fermionic systems of this symmetry class are characterized by a $\mathbb{Z}_{2}$ topological invariant. Computing this so-called winding number from the system parameters can be done in a variety of ways. One possible way of doing so is the Zak phase, defined as the Berry phase over the full Brillouin zone:

[TABLE]

where $|\psi_{k}\rangle$ is the eigenstate of $H$. A gauge transformation $|\psi_{k}\rangle\rightarrow e^{if(k)}|\psi_{k}\rangle$ for any (real, continuous) function $f(k)$ will contribute a multiple of $2\pi$ to the integral, and thus $\nu$ is only defined modulo 2. In the next section, we will discuss how the Zak phase is generalized to open systems.

2.1 Free Lindbladian time evolution

The time evolution of a quantum system in contact with a Markovian bath is described by the Lindblad master equation [28, 29]:

[TABLE]

where $\hat{\mathcal{L}}$ is the Liouvillian superoperator, $H$ is the system’s Hamiltonian and the Lindblad operators $L_{m}$ encode the effect of the bath. If $H$ is quadratic and all $L_{m}$ are linear in fermionic operators, then the system is non-interacting and $\hat{\mathcal{L}}$ can be written in the following bilinear form:

[TABLE]

where $\hat{c}_{j}$ and $\hat{c}^{\dagger}_{j}$ are fermionic superoperators, satisfying $\{\hat{c}_{i},\hat{c}^{\dagger}_{j}\}=\delta_{i,j}$ and acting on a Fock space of operators [21]. In particular, their action on the density operator is:

[TABLE]

where $\hat{P}=e^{i\pi\sum\hat{c}_{j}^{\dagger}\hat{c}_{j}}$ is the parity superoperator. One can show that $\hat{P}$ and $\hat{\mathcal{L}}$ commute, which implies that the Fock space is a direct sum of even and odd subspaces. The bilinear form (7) needs to be defined separately for each of these subspaces. For physical states, we can restrict ourselves to the even subspace, in which the $4N\times 4N$ structure matrix $\mathbf{A}$ takes the following block-triangular form [22]:

[TABLE]

where the Hamiltonian matrix $\mathbf{H}$ and the bath matrix $\mathbf{M}$ have dimensions $2N\times 2N$ and are defined by:

[TABLE]

Here $w_{j}$ are Majorana fermions with anticommutator $\{w_{i},w_{j}\}=2\delta_{i,j}$, as defined by eq. (3). The overall shift $A_{0}=\frac{1}{2}\operatorname{Tr}\mathbf{X}$ in eq. (7) is required for trace conservation.

This formalism has a number of advantages. First, due to the block-triangular form of the structure matrix $\mathbf{A}$, the spectrum of $\hat{\mathcal{L}}$ is completely determined by the real matrix $\mathbf{X}$, specifically its $2N$ eigenvalues $\beta_{i}$. These are known as rapidities and have $\operatorname{Re}(\beta_{i})\geq 0$. If $\operatorname{Re}(\beta_{i})>0$ for all $i$, then there is a unique non-equilibrium steady state $\rho_{ss}$ such that $\hat{\mathcal{L}}\rho_{ss}=0.$ Another advantage is that two-point functions in $\rho_{ss}$ are readily computed by solving the continuous Lyapunov equation [22] for the covariance matrix $\mathbf{C}$:

[TABLE]

Finally, by computing the (generalized) eigenvectors of $\mathbf{A}$, we can construct normal master modes $\hat{b}_{j}$ and $\hat{b}^{\prime}_{j}$, which act as free excitations on the operator Fock space. In this picture, the steady state is the vacuum for the master modes and, assuming $\mathbf{X}$ is diagonalizable, the Liouvillian takes the diagonal form

[TABLE]

If $\mathbf{X}$ is a defective matrix, it is possible to write $\hat{\mathcal{L}}$ in a Jordan normal form with one or more nontrivial Jordan blocks [22]. The eigenmodes of the Liouvillian can then be constructed by acting with a combination of $\hat{b}^{\prime}_{i}$ superoperators onto the steady state. Note that these modes are not density matrices, because they are traceless. To remain in the even sector of the operator Fock space, the number of master modes should be even.

To describe translationally invariant open systems we attach an identical, local bath to each site, represented by the Lindblad operator $L_{j}$. Then, the Hamiltonian and the Lindblad operators can be rewritten as follows:

[TABLE]

where we have two Majorana modes per site $j$ and we consider periodic boundary conditions (PBC) $w_{j+2N}=w_{j}$. Now, one can perform a Fourier transform on $\ell$, to find the $2\times 2$ matrices

[TABLE]

as functions of the quasi-momentum $k\in[-\pi,\pi)$. From these, the matrices $\mathbf{x}(k)$ and $\mathbf{y}(k)$ can be constructed according to eq. (10), which implies in turn that eq. (13) becomes a $2\times 2$ matrix equation for the covariance matrix $\mathbf{c}(k)$ in Fourier space [30].

Applying eqs. (15) and (16) to the Kitaev chain Hamiltonian (4) yields:

[TABLE]

For the bath, we choose general single-site Lindblad operators $L_{j}=\sqrt{g}(c_{j}+\Delta c_{j}^{\dagger})$ where the real parameter $\Delta$ interpolates between pure gain ($\Delta\rightarrow\infty$) and loss ($\Delta\rightarrow 0$). The bath matrix then becomes:

[TABLE]

The Lindbladian dynamics is therefore governed by the matrices:

[TABLE]

where the $\sigma_{\alpha}$ are Pauli matrices. The corresponding spectrum of rapidities for $\mathbf{x}(k)$ is:

[TABLE]

3 Topological order and edge modes

Previous studies on the topology of open quantum systems have largely focused on the steady state(s) [31, 32, 33, 10, 11], in analogy to topological order in the ground state of closed quantum systems. For mixed steady states, however, this is a complicated issue, since the conventional understanding of topological order does not necessarily hold for density matrices, resulting in a poor understanding of finite-temperature topological insulators [9].

For our study of the open Kitaev chain, the Lindblad operators $L_{j}$ break all symmetries of the Hamiltonian [34], resulting in a unique steady state independent of initial conditions. This implies that the ground state degeneracy of the closed system, caused by the zero-energy edge modes, does not survive the long-time limit [35]. However, we will see that the full dynamics of the system retains some of its original symmetries and, following studies of topology for non-Hermitian Hamiltonians [13, 15] and driven systems[3], we can extract some topological properties by studying the spectrum and eigenmodes of the Liouvillian operator $\hat{\mathcal{L}}$.

3.1 Topological properties of $\mathbf{x}(k)$

Since the spectrum of $\hat{\mathcal{L}}$ is completely determined by the matrix $\mathbf{x}(k)$, it is worth exploring its symmetry properties. Actually, one can show that the matrices $\mathbf{A}$ and $\mathbf{X}$ share the same topological invariants (see Appendix C). For sake of clarity, let us start by considering the case of pure loss $\Delta=0$. Here, $\mathbf{x}(k)$ and $\mathbf{h}(k)$ only differ by a constant shift and multiplicative factor, neither of which affects its eigenvectors. This implies that the winding number, given by eq. (5), of the closed system is well-defined and unchanged. One can similarly show that the same holds for the case of pure gain and for separate loss and gain channels on each site. The authors of [36] found a similar result for a dissipative SSH model, in which the topological band structure of the Liouvillian was identical to that of the Hamiltonian.

For $\Delta\neq 0$, even though the chiral symmetry of $\mathbf{x}(k)$ is broken, there remains an anti-unitary symmetry (AUS) [37, 38], up to a constant shift:

[TABLE]

where the operator $\mathcal{S}=K\sigma_{x}$ corresponds to a complex conjugation $K$ and a swapping between the even and odd Majorana modes, while the constant shift $a_{0}$ is related to the trace of $\mathbf{x}(k)$, that is $a_{0}=A_{0}/N=\frac{1}{2}\operatorname{Tr}\mathbf{x}(k)=2g(1+\Delta^{2})$. Now, consider the biorthogonal left and right eigenvectors of $\mathbf{x}(k)$ [39]:

[TABLE]

The symmetry (21) automatically implies that

[TABLE]

Now there are two options: either $\underline{u}_{i}(k)$ is also an eigenvector of $\mathcal{S}$, which implies that $\operatorname{Re}\beta_{i}(k)=a_{0}$ and $\operatorname{Im}\beta_{i}(k)\neq 0$; or $\underline{u}_{i}(k)$ and $\mathcal{S}\underline{u}_{i}(k)$ are the two different eigenvectors of $\mathbf{x}(k)$, with eigenvalues $\beta_{i}(k)$ and $2a_{0}-\beta_{i}^{*}(k)$, respectively.

We say that the AUS is unbroken if the first case holds for all $k$ in the Brillouin zone. Otherwise, the symmetry is said to be spontaneously broken. Examples of the rapidity spectrum with unbroken and broken AUS can be seen in panels (a) and (b) of Figure 1. If $\operatorname{Re}\beta_{i}(k)=a_{0}$ for a subset of the Brillouin zone, as is the case in panel (b), then this subset is delimited by a pair of so-called exceptional points (EPs), marked with a red diamond in Figure 1. At these exceptional points the matrix $\mathbf{x}(k)$ becomes non-diagonalizable. Although the present analysis is reminiscent of $\mathcal{PT}$-symmetric Hamiltonians [40, 23], where the eigenvalues are real when the symmetry is unbroken, in our case the conjugation $K$ cannot be seen as time reversal as we are considering a non-unitary time evolution. Thus we refrain from labeling the AUS as a $\mathcal{PT}$ symmetry, although they are mathematically equivalent.

For the open Kitaev chain111In the present study, we assume $|\gamma|>\sqrt{\frac{g\Delta}{J}}J$ and $g\Delta<J$. The latter is warranted since we typically require weak dissipation for the Lindblad formalism to apply. If the former is not satisfied, then the gap can close at values of $k$ other than [math] or $\pi$, which increases the size of the AUS-broken regions. If $\gamma<g\Delta$, the topologically nontrivial phase will disappear completely. See Appendix A for more details., the matrix $\mathbf{x}(k)$ has broken AUS whenever $2J-2g\Delta\leq|\mu|\leq 2J+2g\Delta$. This implies that between the two original phases of the Kitaev chain, an intermediate region of width $4g\Delta$ appears in which the band gap of $\mathbf{x}(k)$ is closed, as the two bands touch at the exceptional points. This is consistent with the findings of [13].

Without a chiral symmetry, it is not clear whether the $\mathbb{Z}_{2}$ topological invariant that characterizes the topology of the closed Kitaev chain is well defined, i.e. a quantized winding number. To investigate this, we generalize the winding number (5) to the biorthogonal basis (22), yielding:

[TABLE]

where $\nu_{i}$ is the winding number associated to the band $i$ of $\mathbf{x}(k)$. In general, these numbers will not be real or quantized when chiral symmetry is broken. However, in the regime with unbroken AUS, one can show that their real part is quantized. To see this, note that the unbroken symmetry implies that:

[TABLE]

where the eigenvalues must lie on the unit circle, since $\mathcal{S}$ is antiunitary. Using eq. (24), we can show that the winding number satisfies the following relation:

[TABLE]

If we demand the eigenvalue $e^{i\varphi_{i}(k)}$ to be single-valued within the Brillouin zone, then the last term must be an even integer. Therefore we must have $\nu_{i}=-\nu_{i}^{*}+2n$, or equivalently $\operatorname{Re}(\nu_{i})=n$ for $n\in\mathbb{Z}$. Requiring gauge invariance again restricts this to $\mathbb{Z}_{2}$. From here on we will simply refer to the quantized real part as the winding number. Whether the imaginary part also has a physical interpretation in this context is not clear. We also note that $\sum_{i}\nu_{i}\mod 2=0$, such that it is enough to study the winding number of only one of the two rapidity bands.

The resulting phase diagram of the open Kitaev chain can be found in Figure 1. As we can see in panel (d), we have now three phases: A trivial topological phase (denoted in the figure as $I$), a non-trivial topological phase (indicated as $II$) and, finally, an intermediate phase (indicated as $III$). In the non-trivial topological phase the AUS is unbroken and the real part of the winding number (for further details on its derivation see Appendix B) has quantized value $1$, while in phase $I$ its value is zero. Let us now explore the effect of the bath on the bulk-boundary correspondence; whether a nonzero winding number still implies the existence of topologically protected edge modes, and what those modes might look like.

3.2 Bulk-boundary correspondence and its breakdown

While the bulk-boundary correspondence seems to generally fail for non-Hermitian systems [41], this is not the case for $\mathcal{PT}$-symmetric Hamiltonians [42]. As we have argued before, the open Kitaev chain has an AUS that is mathematically equivalent to a $\mathcal{PT}$ symmetry, and therefore we expect bulk-boundary correspondence to hold in our case as well. That is to say, under OBC we expect two eigenvectors of $\mathbf{X}$ to be localized at the edges when $\nu=1$, but none when $\nu=0$. We have investigated this in two complementary ways. Firstly, we have analyzed the modes close to the gap in the continuum limit and studied the effect of domain walls which give rise to a localized edge mode (details of this approach can be found in Appendix D). Secondly, we have studied numerically the spectrum and eigenmodes of the $2N\times 2N$ matrix $\mathbf{X}$ with OBC. The resulting spectrum is shown in Figure 2.

Both of these analyses indicate the presence of two exponentially localized edge modes when $|\mu|<2J$, i.e. in the nontrivial phase of the closed Kitaev chain. This can be clearly seen in Figure 2, where the red points mark the rapidities of the edge modes at $\beta=g(1\pm\Delta)^{2}$, consistent with the derivation in Appendix D. The corresponding eigenvectors with localization length $\xi=2\gamma/(2J-|\mu|)$ can be seen in the insets of that figure for one value of $\mu$. Comparing these results to the winding number of Figure 1(d), we conclude that the bulk-boundary correspondence holds as long the AUS is unbroken. In that case $\nu$ is quantized and accurately reflects the presence or absence of edge modes. In the intermediate region $2J-2g\Delta<|\mu|<2J+2g\Delta$ where the AUS is broken, the bulk-boundary correspondence breaks down. Even though the band gap is closed, the topologically protected edge modes remain as long as $|\mu|<2J$. Only when $\mu$ reaches $\pm 2J$, the real part of the bulk bands connects to the edge modes and they disappear, as shown in Figure 2. Interestingly, the winding number does not even show a discontinuity at this point.

3.3 Physical interpretation of edge master modes

An important question remains: what is the physical meaning of these edge modes? How could they be observed experimentally? In the closed Kitaev chain, Majorana edge modes can in principle be observed by scanning tunneling microscopy [43]. The localized eigenvectors of $\mathbf{X}$, on the other hand, are not physical excitations on top of a ground state, and thus one cannot expect to detect them in the same way. Instead, they correspond to boundary master modes, as defined in eq. (14). In terms of the left and right eigenvectors of $\mathbf{X}$, $\underline{v}_{i}$ and $\underline{u}_{i}$ respectively, and the covariance matrix $\mathbf{C}$, these master modes $\hat{b}^{\prime}$ and $\hat{b}$ are given by (see Appendix C for details)

[TABLE]

A localized eigenvector then implies that the corresponding $\hat{b}^{\prime}$ is localized too, while localization of $\hat{b}$ also requires that there are no long-range correlations in the steady state.

Now, define $\hat{b}^{\prime}_{L}$ and $\hat{b}^{\prime}_{R}$ to be the two edge master modes, corresponding to the localized eigenvectors. In the limiting case $\mu=0$ and $J=\gamma$, where the Majorana edge modes decouple completely from the bulk, we can analyze the system’s behavior exactly. The edge master modes then take the simple form: $\hat{b}^{\prime}_{L}=\hat{c}^{\dagger}_{1},~{}~{}~{}\hat{b}^{\prime}_{R}=\hat{c}^{\dagger}_{2N}$. All the other master modes are restricted to the bulk, since they must be orthogonal to the edge modes. To understand the physical consequences of this type of localization, let us consider the expectation value of the edge occupation number in this limit:

[TABLE]

Taking as initial condition the (pure) ground state with the edge mode occupied, i.e. $\langle\mathcal{O}(0)\rangle=1$, we can then turn on the dissipation and study the time evolution. Since the edges are decoupled, only the decay mode $\hat{b}^{\prime}_{L}\hat{b}^{\prime}_{R}\rho_{ss}$ contributes to this expectation value, with a decay rate $2g(1+\Delta^{2})$ given by its eigenvalue. In the steady state, we find $\lim_{t\rightarrow\infty}\langle\mathcal{O}(t)\rangle=0$, as can be computed using eq. (13). Therefore the time evolution is simply given by:

[TABLE]

In other words, the Majorana zero modes of the Kitaev chain are still present in the system, but acquire an exponential decay. These results are in agreement with those found in [35] in a similar context.

If we look at the evolution of local observables, such as the single-site fermionic occupation number, then the fact that the edge master modes decay at different rates becomes apparent. This is shown in Figure 3, where we plot the time evolution of the local density at the edges (dashed red line for the left edge, and solid blue line for the right one) and at the center of the chain (green dash-dot line) (details of this derivation are described in Appendix E). The occupation at the right edge of the system approaches its steady state value more slowly than the occupation number at the left edge. The difference between these decay rates is controlled by the parameter $\Delta$, responsible for the split between the real parts of the edge mode rapidities, and it vanishes for $\Delta=0$. In the middle of the chain, neither of the edge modes contributes to the single-site occupation and the decay rate is halfway between that of the two edges. Since the decay modes in the bulk all have complex rapidities, there is also an oscillating part to the expectation values.

Note that, since the edge modes have rapidities $\beta=g(1\pm\Delta)^{2}$, something interesting happens at $\Delta=\pm 1$: one of the rapidities vanishes. This means that one of the pair of edge modes does not decay and the steady state becomes degenerate. Intuitively, this makes sense: the Lindblad operators $L_{j}$ are now proportional to one of the Majorana operators $w_{2j}$ or $w_{2j-1}$, and will consequently leave one of the Majorana edge modes isolated from the bulk. However, this degenerate steady state $\rho_{ss}^{\prime}$ is restricted to the odd sector of the operator Fock space and, as such, it can never be reached from a physical state in the even sector, implying that is not a true degeneracy. One way around this might be to consider a trivial domain ($|\mu|>2J$) within the bulk of a nontrivial Kitaev chain with OBC. With two pairs of edge modes at the four different domain walls, it is possible to construct a steady state degeneracy in the even sector, where two Majorana edge modes decay and two remain.

3.4 Steady state expectation values

As we have already argued (see Figure 3), the presence of edge modes affects the rates at which some observables approach their steady state expectation values. A logical next step would be to determine whether the steady state expectation values themselves also reflect some topological order. Although the Majorana edge modes decay under influence of the bath, they might leave some signature in the steady state covariance matrix. By numerically solving the Lyapunov equation (13) with OBC, we can extract the expectation values of various observables. In particular, it is possible to study the expectation value of the occupation number at the edges, as a function of $\mu$, to see whether the phase transition could be detected through these observables. In order to compare the occupation at the edges relative to the bulk, we define the edge occupation ratios:

[TABLE]

We expect this quantity to be nonzero even in the absence of topological edge modes, due to the boundary conditions. However, when plotted as a function of the chemical potential, these ratios have a clear peak in the topologically nontrivial phase. This can be seen in the panels (a) and (b) of Figure 4, where we show the numerical results for two different values of the dissipation strength $g$. In the limit $g\rightarrow 0$, there are sharp kinks at the critical points $\mu=\pm 2J$ and the left and right edge show identical behavior. As the dissipation strength is increased, the phase transition is smoothened out and, for $\Delta\neq 0$, there is an increased asymmetry between the left and right edge, as can be verified in panel (b).

This asymmetry can be seen as an unequal charge build-up on the edges of the system. Therefore one might expect that, if we consider PBC instead of OBC, there would be a current flowing in the steady state. One way to define the electronic current is:

[TABLE]

The steady state expectation value of this operator can be seen in panels (c) and (d) of Figure 4, again for the same two values of $g$. As with the edge occupation ratio, there are sharp kinks at the transitions between the topologically trivial and the nontrivial regimes when $g\rightarrow 0$, and stronger dissipation smoothens out this transition. The expectation value of the current operator also increases with the strength of the dissipation. The current reverses direction if we change the sign of $\Delta$ and vanishes at $\Delta=0$.

We stress that these observables are not order parameters: in the open system, there is no sharp transition or discontinuity at the critical value of $\mu$. However, in the limit $g\rightarrow 0$, kinks appear and the derivatives become discontinuous. Furthermore, we have verified the presence of a steady state current under PBC, which appears to susceptible to topological order. This may seem counter-intuitive at first glance, but note that the $p$-wave pairing terms in the Hamiltonian (1) are not symmetric under spatial inversion, as they pick up a minus sign. It is therefore natural that the steady state does not have this symmetry.

4 Higher winding numbers in the driven-dissipative Kitaev chain

Driven closed quantum systems with chiral symmetry have been shown to accommodate multiple pairs of edge modes, of which one can distinguish two different types [25, 3]. Topological phases are in this case characterized by a $\mathbb{Z}\times\mathbb{Z}$ invariant. A relevant question is then what happens to these properties when a bath is coupled to a driven system. This has been somewhat explored in the $XY$ Heisenberg spin chain [8], which maps to the Kitaev chain by a Jordan-Wigner transformation [44]. These studies were further extended to identify the existence of Floquet Majorana edge modes [45].

We will now explore what remains of the structure previously unveiled in the open Kitaev chain, once driving is added. In particular we will look into greater detail at the non-equilibrium properties of these Floquet Majorana edge modes and their classification into two types as described in [25, 26]. We show that the combined effects of the driving and the bath induce regions of higher winding number which are separated by intermediate phases of broken AUS.

4.1 Topology and the Floquet formalism

As it is well known, a closed periodically driven system, with $H(t+T)=H(t)$, can be studied using Floquet theory [46]. In this case, the stroboscopic time evolution is generated by the Floquet Hamiltonian $H_{F}(t_{0})$, defined as:

[TABLE]

Similarly to what happens in a periodic lattice, where quasi-momenta can be restricted to the Brillouin zone, the eigenvalues $\epsilon_{F}$ of $H_{F}(t_{0})$ are only defined up to multiples of $2\pi/T$, resulting in a Floquet Brillouin zone with $\epsilon_{F}\in[-\pi/T,\pi/T)$. For this reason these eigenvalues are referred to as quasi-energies. Note that $H_{F}(t_{0})$ depends on the initial time $t_{0}\in[0,T)$ for which the Floquet propagator (32) is defined, but the quasi-energy spectrum does not [46].

The periodicity of quasi-energies in a driven two-band system with chiral symmetry results in the presence of two band gaps, one at $\epsilon_{F}=0$ and another one at $\epsilon_{F}=\pi/T$. In each band, topologically protected edge modes may arise, characterised by their own winding number. It has been shown that the driven Kitaev chain has a $\mathbb{Z}\times\mathbb{Z}$ topological invariant [25, 3] – an expected result, since the periodic driving introduces an additional $\mathcal{S}_{1}$ manifold.

In order to see this, we need to determine whether the Floquet Hamiltonian $H_{F}(t_{0})$ inherits the chiral symmetry of the instantaneous Hamiltonian $H(t)$. For sake of clarity we will consider periodic quenches of the chemical potential $\mu$, alternating between the values $\mu_{1}$ during a time $t_{1}$, and $\mu_{2}$ during a time $t_{2}$, with $T=t_{1}+t_{2}$ being the driving period. It is easy to check that the Floquet Hamiltonian in this scenario has chiral symmetry only at two particular points $t^{\prime}$ and $t^{\prime\prime}$ within each period, precisely where the driving protocol is time-reversal invariant - in our case, at the center of the two time plateaus. Let us denote the corresponding Floquet Hamiltonians as $H_{F^{\prime}}=H_{F}(t^{\prime})$ and $H_{F^{\prime\prime}}=H_{F}(t^{\prime\prime})$, with

[TABLE]

These two effective Hamiltonians share the same spectrum, but can have different bulk winding numbers $\nu^{\prime}$ and $\nu^{\prime\prime}$, since they are related by the $k$-dependent unitary transformation $G=e^{-iH_{2}t_{2}/2}e^{-iH_{1}t_{1}/2}$. Using this transformation, one is able to derive that:

[TABLE]

where $|\psi_{F^{\prime}}(k)\rangle$ are eigenstates of $H_{F^{\prime}}$. One can show that this $\mathbb{Z}$ topological invariant counts the number of Floquet Majorana edge modes with quasi-energy $\epsilon_{F}=\pi/T$ under open boundary conditions [25, 47]. Note that, unlike the winding number of undriven systems given by eq. (5), $\nu_{\pi}$ does not have to be defined modulo 2 in order to be gauge invariant, which implies that more than two Floquet Majorana edges modes can accumulate at the same band gap. Similarly, the other $\mathbb{Z}$ invariant, which is given by the combination $\nu_{0}\equiv(\nu^{\prime}+\nu^{\prime\prime})/2$, counts the edge modes with $\epsilon_{F}=0$. Therefore, by driving the system, one can get multiple orthogonal edge modes at each end of the chain in contrast to a single one in the undriven case.

4.2 Lindblad-Floquet theory

To see how the two new topological invariants are affected by the presence of a bath, we first need to generalize Floquet theory to encompass the Lindbladian time evolution. For a time-dependent Liouvillian, the time-evolution superoperator over one period is $\hat{\mathcal{U}}(t_{0}+T,t_{0})=\hat{T}\exp\left(\int_{t_{0}}^{t_{0}+T}\,\mathrm{d}\tau\hat{\mathcal{L}}(\tau)\right)$. Whether this can be written in terms of a Floquet Liouvillian $\hat{\mathcal{L}}_{F}$ depends on the details of the system [48], but in our case it does not pose a problem and therefore we write

[TABLE]

where $\hat{\mathcal{L}}_{F}(t_{0})$ can be expressed, as done in eq. (7), in terms of a effective structure matrix $\bm{A}_{F}$. Next, to preserve the antiunitary symmetry (21) of the dissipative Kitaev chain we follow the procedure described in eq. (33). We pick as the starting time $t_{0}$ one of the two time-reversal invariant points $t^{\prime}$ and $t^{\prime\prime}$ in the two-step driving protocol and write:

[TABLE]

The computation of the off-diagonal block $\mathbf{Q}^{\prime}$ is somewhat complicated and shown in Appendix F. Similar equations define the matrices $\mathbf{A}_{F^{\prime\prime}}$, $\mathbf{X}_{F^{\prime\prime}}$, and $\mathbf{Q^{\prime\prime}}$ for the Floquet Liouvillian $\hat{\mathcal{L}}_{F^{\prime\prime}}$ at time $t^{\prime\prime}$. Our main focus will be on the properties of $\mathbf{X}_{F^{\prime}}$ and $\mathbf{X}_{F^{\prime\prime}}$ as their eigenvectors determine the topological properties of the driven-dissipative system. Their eigenvalues, which we refer to as quasi-rapidities, are defined up to a multiple of $2\pi i/T$ and can therefore be restricted to the first Floquet Brillouin zone: $-\pi/T<\operatorname{Im}(\beta)\leq\pi/T$. For PBC, translational invariance allows us to decompose the matrices $\mathbf{X}_{F}$ in Fourier components of $2\times 2$ matrices $\bm{x}_{F}(x)$. It turns out that if the Fourier components $\mathbf{x}_{1}(k)$ and $\mathbf{x}_{2}(k)$ satisfy the AUS condition (21), so do the matrices $\mathbf{x}_{F^{\prime}}(k)$ and $\mathbf{x}_{F^{\prime\prime}}(k)$, due to the symmetric composition of the exponentials in eq. (37). Therefore, from the result we discussed for the undriven case, we can expect regions in parameter space where the AUS is unbroken and the winding numbers $\nu_{0}$ and $\nu_{\pi}$ are quantized, separated by intermediate regions containing exceptional points.

As we will see, in this case, the phase diagram turns out to be much richer: aside from the known phases of the undriven Kitaev chain described in Sec. 2, there are new transitions induced by Floquet resonances. Such resonant regimes are characterized by long-range order [8], nonlocal Floquet Hamiltonians [24] and the possibility of high winding numbers [45].

To obtain an expression for the spectrum of $\mathbf{x}_{F^{\prime}}(k)$ and $\mathbf{x}_{F^{\prime\prime}}(k)$, as well as the topological winding numbers $\nu_{0}$ and $\nu_{\pi}$, we proceed as follows. Using Euler’s formula we write

[TABLE]

with $\underline{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z}$), and where $\underline{m}$, $\underline{n}$ and $\underline{p}$ are complex valued vectors. Once again we have $a_{0}=g(1+\Delta^{2})$. After some algebraic manipulation, one can show that:

[TABLE]

The quasi-rapidity spectrum of $\mathbf{x}_{F^{\prime}}(k)$ is then given by $\beta_{F}(k)=a_{0}\pm i\cos^{-1}(p_{0})/T$, and its winding number is the same as that of $\underline{p}(k)\cdot\underline{\sigma}$. Moreover, the AUS guarantees that $p_{z}$ is real while $p_{x}$ and $p_{y}$ are purely imaginary. Starting from eq. (24) and following the derivation in Appendix B, we finally obtain an expression for the winding number:

[TABLE]

The same calculation can be done for $\mathbf{x}_{F^{\prime\prime}}(k)$ at the other time-reversal invariant point to obtain the winding number $\nu^{\prime\prime}$ . By adding and subtracting the resulting winding numbers, we finally find the two topological invariants $\nu_{0}$ and $\nu_{\pi}$ of the driven-dissipative system.

4.3 Numerical results

In order to build some intuition for the problem, let us first consider the infinite-frequency limit ($T\rightarrow 0$), for which the stroboscopic time evolution is governed by the average Liouvillian [7]. In this case, that is simply the dissipative Kitaev chain with $\mu=(\mu_{1}t_{1}+\mu_{2}t_{2})/T$. As the driving period is increased, the Floquet Brillouin zone shrinks. Eventually, the quasi-rapidity bands reach the edges of the Floquet Brillouin zone, at which point the gap at $\pi/T$ closes. As we further increase $T$, we first encounter a finite intermediate region with broken AUS, after which the gap reopens and the first Floquet resonance occurs. This phase is characterized by avoided crossings of the quasi-rapidity bands and nonlocal hopping in the effective generator of (stroboscopic) time evolution [24]. According to the bulk-boundary correspondence, if we were to use OBC, an edge mode would be present inside the reopened gap. For large enough $T$, the resonances become more frequent and start to overlap, i.e. several resonances can be present at once. This is the mechanism that leads to multiple edge modes per gap.

The behavior described above can be seen in the bottom panel of Figure 5, which shows the imaginary part of the quasi-rapidities as a function of the period. As the Floquet resonances accumulate, the number of edge modes in each gap can be clearly seen from the real part of the spectrum, shown in the middle panel of the same figure. Away from the AUS-broken regions, these numbers neatly coincide with the two bulk winding numbers $\nu_{0}$ and $\nu_{\pi}$, computed numerically from eq. (41) and shown in the top panel of Figure 5. The real part of the quasi-rapidities, which arises purely as a consequence of non-Hermiticity, is more than a visual aid to distinguish the different edge modes: in the absence of chiral symmetry, it provides the mechanism that prevents two edge modes from recombining and scattering, even when they occupy the same edge and the same gap. As long as the real parts of the quasi-rapidities are different, the corresponding modes cannot hybridize.

However, the hybridization of two edge modes can occur when they have the exact same quasi-rapidity. At that point, the imaginary parts of the pairs of edge modes can repel each other, resulting in two pairs with the same real part. These hybridized edge modes have an imaginary part that is neither [math] nor $\pi/T$ and are not counted by the winding numbers, although they are still localized at the edges. Whether they are also topologically protected in some way is an interesting open question. If they are, then this can be seen as further violation of the bulk-boundary correspondence. These hybridized edge modes can be observed for certain values of $T$ as additional black lines in the real part of the spectrum in the middle panel of Figure 5. They are present around $T{\approx}3.7$ (enlarged in an inset), as well as for $T>4.5$. A similar phenomenon seems to happen around $T\approx 4.2$, where $\nu_{\pi}$ jumps from $-2$ to [math] while the bulk band gaps do not close. In the real part (also enlarged in an inset), we see two pairs of edge modes (blue lines) coming together and disappearing.

In order to better understand the behavior of the winding numbers, we show in Figure 6 a density plot of $\nu_{0}$ (left panel) and $\nu_{\pi}$ (right panel) in the $(T,\mu_{1})$ plane. Regions with broken AUS are shown as hatched areas. We can see that, for $T$ slightly larger than 4, one of the winding numbers changes by $2$ without an broken-AUS intermediate region. This supports the notion that there may be other mechanisms of AUS-breaking responsible for the noise around $\operatorname{Re}(\beta)=a_{0}$ seen in Figure 5 for $T>3$ which only appears under OBC. A better understanding of these phenomena will be left for future work.

Aside from the anomalous transitions described above, the phase diagram of the driven-dissipative Kitaev chain (as shown in Figure 6) has the same rib structure of Floquet resonances that appears in previous work [8, 45, 24]. The total number of edge modes, given by $|\nu_{0}|+|\nu_{\pi}|$, corresponds to that found by [45]. Although the splitting into $\nu_{0}$ and $\nu_{\pi}$ is clear from the viewpoint of topological classification, it is less obvious what the physical difference is between the two types of edge modes. It is not clear if these modes can be qualitatively distinguished by an appropriately chosen physical observable.

So far, for the driven-dissipative system, we have analyzed the quasi-rapidity spectrum and its link to topological phases. Presently, we will consider the impact of the transitions between these phases on the expectation values of observables in the Floquet steady state. These values can be computed by solving the following discrete-time Lyapunov equation (see Appendix F for further details):

[TABLE]

Here $\mathbf{C}_{F^{\prime}}$ is the stroboscopic covariance matrix in the infinite time limit at the first time-reversal invariant point, i.e. $t_{0}=t^{\prime}$. From this covariance matrix, we can derive stroboscopic expectation values of all other observables. Once again we consider the steady state current $\mathcal{J}$ under PBC, studied in section 3.4 for the undriven system. The results are shown in Figure 7 for two values of the system-bath coupling $g$. As we can appreciate, with weak dissipation (left panel), the rib-like resonance structure is clearly visible and there appears to be a direct correspondence between the steady state current and the total winding number $|\nu_{0}|+|\nu_{\pi}|$. On the other hand, the stronger the dissipation the smoother the transitions become, so that for a strong enough dissipation this rib-like structure is barely visible (right panel).

5 Conclusions and outlook

By using tools from the study of non-Hermitian Hamiltonians, we have analyzed the topological properties of a driven-dissipative extension to the Kitaev chain. Dissipation is generated by an identical Markovian bath coupled locally to each site of the system, as described by a Lindblad master equation. By introducing a parameter $\Delta$ to the bath that interpolates between gain and loss, we can study a wide range of effects. For example, intermediate values of $\Delta$ break the chiral symmetry that protects topological order in the Kitaev chain, but create a new antiunitary symmetry (AUS) reminiscent of $\mathcal{PT}$ symmetry in non-Hermitian Hamiltonians. While this AUS leaves intact the different topological phases of the Kitaev chain, it introduces intermediate regions between the phases, where the AUS is broken and the winding number is no longer quantized. In these regions, the bulk-boundary correspondence breaks down.

The addition of periodic driving makes for some incredibly rich physics and we have only explored the tip of the iceberg. The phase diagram of the driven-dissipative system mostly follows the regions defined by Floquet resonances, which can overlap to produce higher winding numbers and multiple pairs of edge modes. Having multiple modes localized on each edge can lead to interesting phenomena such as hybridization. Precisely under which conditions these hybridized edge modes exist and what their relation is to the bulk topological order are still open problems. Like in the undriven open system, topological phases are separated by intermediate AUS-broken regions. However, there also seem to be other mechanisms of AUS-breaking, connected to the driving and the OBC, that are not yet understood.

The topological edge modes that we have studied, both in the driven and undriven case, are transient properties of the dissipative system. The presence of the bath leads the system to a unique state at long times and bestows decay rates onto the edge modes. Interestingly, these decay rates can be different for the two edges of the chain, resulting in a left-right asymmetry that is reflected in the dynamics and steady state expectation values of various observables. Numerical results show some signatures of the topological order in the steady state, but an analytical connection has eluded us so far. Various ways to define winding numbers for density matrices have been proposed, for example using the Uhlmann phase [49], the interferometric phase suggested in [50], and more recently the Ensemble Geometric Phase [12]. While we stress that these mixed-state winding numbers are very different from the ones we have calculated for the band Liouvillian, it would be intriguing to apply them to the steady state of the dissipative Kitaev chain and compare this to our results.

Another promising angle for linking Floquet topological order to observables is borrowed from the study of spin chains and known as stroboscopic spin textures [45, 51, 26] or nonlocal string order [52, 53], although it might not work with broken chiral symmetry. If applicable to our system, these techniques could potentially distinguish the two different winding number $\nu_{0}$ and $\nu_{\pi}$ in the driven system.

Finally, as potential experimental realizations of Majorana zero modes have been proposed as information processing devices to implement unitary gate operations in topological quantum computation [54], the generation of higher number Majorana modes using periodic systems opens the door to further explore in which way these new localized modes could optimally be used in quantum computation together to its robustness against the presence of external noise. Moreover, a thorough understanding of the decay times of Majorana modes under influence of the environment may prove important for future developments in this field.

Acknowledgments

This work is part of the Delta-ITP consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW) and of the research program of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organization for Scientific Research (NWO). The authors want to thank V. Gritsev, E. Ilievski, N. Robinson and J. van Wezel for guidance and interesting discussions.

Appendix A Geometric interpretation of the winding number

A $2\times 2$ matrix with antiunitary symmetry given by (21) can only have the following form:

[TABLE]

where $a_{0}$, $b(k)$, $c(k)$ and $d(k)$ are real for all $k$. The matrix becomes defective when $b^{2}+c^{2}=d^{2}$ and this equation forms the surface of a double cone in the space spanned by $b$, $c$ and $d$. Inside of the cone, the antiunitary symmetry is spontaneously broken.

As $k$ traverses the Brillouin zone, $\mathbf{x}(k)$ defines a loop in this three-dimensional space, characterized by the system parameters such as $J$ and $\mu$. If the loop intersects the cone, then it does so at a pair of so-called exceptional points. Those are the points where the band gap closes. We can now define a topological invariant as the winding number of $\mathbf{x}(k)$ around this defective cone, interpreting it as an obstruction in $\mathbb{R}^{3}$. If the loop winds around the cone, it cannot be contracted without closing the gap (see Figure 8).

It is also clear from this geometric picture how $\gamma$ affects the transitions between the different phases. The loop defined by $\mathbf{x}(k)$ is an ellipse with semi-axes of length $2J$ and $2\gamma$. The radius of curvature at its vertex ($k=0$ or $\pi$) is $2\gamma^{2}/J$. If this is larger than $2g\Delta$, i.e. the radius of the cone, then the cone will touch the ellipse at its vertex when $|\mu|=2J\pm 2g\Delta$. In the main text, we therefore require $|\gamma|>\sqrt{Jg\Delta}$. If $|\gamma|$ is smaller, the ellipse becomes too eccentric and will touch the cone for lower values of $|\mu|$. In that case one can have four exceptional points, as the ellipse intersects the cone in four different places. Furthermore, it is obvious that the ellipse cannot wind around the cone at all when $|\gamma|<g\Delta$, in which case the topologically nontrivial phase disappears entirely.

Appendix B Computing the winding number

In order to computer the winding number of $\mathbf{x}(k)$, we define:

[TABLE]

The left and right eigenvectors of $\mathbf{x}(k)$ can then be written as:

[TABLE]

such that $\underline{v}_{i}^{*}\cdot\underline{u}_{j}=\delta_{ij}$. The two winding numbers (24) associated with these bands are:

[TABLE]

Note that the first term $\oint\!\,\mathrm{d}\theta$ is simply the winding number of $q(k)$ around the origin in the complex plane, which corresponds to the winding number of the closed Kitaev chain. The second term only contributes to the real part of $\nu$ in the regime where $\mathcal{PT}$ symmetry is broken. We can further simplify:

[TABLE]

The numerical results of the integral are shown in Figure 1(d).

Appendix C Topology of X vs. A

The matrix $\mathbf{X}$ does not give all the necessary information about the full Liouvillian time evolution, despite determining the entire spectrum. The decay modes and the steady state also depend on the off-diagonal block $\mathbf{Y}$. To build a full description of the time evolution, we also need to know the eigenvectors of $\mathbf{A}$, which form the so-called normal master modes, superoperators that map one decay mode into another. These eigenvectors can be build up from the eigenvectors of $\mathbf{X}$ and the steady state covariance matrix $\mathbf{C}$, which in turn depends on $\mathbf{Y}$. To see this, we can write:

[TABLE]

Therefore the left and right eigenvectors of $\mathbf{A}$, denoted by $\underline{\phi}_{i}^{\pm}$ and $\underline{\psi}_{i}^{\pm}$ respectively, can be written as:

[TABLE]

where $\underline{v}_{i}$ and $\underline{u}_{i}$ are the left and right eigenvectors of $\mathbf{X}$, as given in (22). The right eigenvectors $\underline{\psi}_{i}^{\pm}$ provide the normal master modes of eq. (27). All of the above can be written in Fourier space, as functions of the quasimomentum $k$, and this allows us to compute the winding numbers for the four bands of the matrix $\mathbf{A}(k)$:

[TABLE]

where $\nu_{i}$ are the winding numbers associated with the bands of $\mathbf{x}(k)$, as given by eq. (24). Integration by parts shows that the two results are complex conjugates. This provides a thorough justification to only study the topological band structure of $\mathbf{x}(k)$.

Appendix D Edge modes of $\mathbf{X}$

One way to prove the presence of localized edge modes is to consider the linearized Dirac form of $\mathbf{x}(k)$ around $k=0$:

[TABLE]

A Fourier transform to real space yields

[TABLE]

where we have defined the Dirac mass $m\equiv-2J-\mu$. Now we consider $m$ to vary in space, such that $m(r)=m_{1}\Theta(-r)+m_{2}\Theta(r)$, with $\Theta(r)$ the Heaviside step function. An eigenvector $\underline{u}(r)$ with eigenvalue $\beta$ must satisfy:

[TABLE]

with $\delta\beta\equiv g(1+\Delta^{2})-\beta$. The solution of this differential equation is of the form $\underline{u}(r)=\exp\left[\int_{0}^{r}\!\,\mathrm{d}r^{\prime}\mathbf{B}(r^{\prime})\right]\underline{u}(0)$. A localized solution is only possible when the eigenvalues of $\mathbf{B}(r)$ have a real part that changes sign at $r=0$, such that $\underline{u}(r)$ falls off exponentially on either side of the domain wall. This is the case when $\delta\beta=\pm 2g\Delta$ and $m_{1}$ and $m_{2}$ have different signs. Since changing the sign of $m$ puts us in a different topological phase, this proves that an edge mode appears on the boundary between two phases222The same analysis can be done around $k=\pi$, in order to find the phase transition at $\mu=+2J$.. The eigenvalues corresponding to these localized solutions are:

[TABLE]

where the sign depends on the orientation of the domain wall.

Finally, note that we can model the vacuum as a fermionic chain with $\mu\rightarrow-\infty$, placing it firmly in the topologically trivial regime (as expected). We can therefore use this same analysis to describe edge modes of the Kitaev chain with open boundary conditions. As long as $|\mu|<2J$, each edge forms a domain wall and allows for a localized edge mode.

Appendix E Time evolution of observables

The time evolution of the covariance matrix $\mathbf{C}(t)$ can be expressed in terms of the eigenvalues and eigenvectors of $\mathbf{X}$. The covariance matrix satisfies the following differential matrix equation [8]:

[TABLE]

Assuming $\mathbf{X}$ is diagonalizable and that the steady state is unique, we find as a solution:

[TABLE]

where $\mathbf{C}_{ss}$ is the covariance matrix in the steady state, found by solving eq. (13) (with $\mathbf{C}_{ss}=\mathbf{C}$), and $\mathbf{C}_{0}$ is the covariance matrix for the initial state. In the limit $t\rightarrow 0$, we recover $\mathbf{C}(0)=\mathbf{C}_{0}$ by identifying two resolutions of the identity in our biorthogonal basis. If the initial state is a pure energy eigenstate (e.g. the ground state), we can compute $\mathbf{C}_{0}$ from the eigenvectors of the $2N\times 2N$ matrix $\mathbf{H}$:

[TABLE]

with $\mathbf{H}\underline{\psi}_{n}=\epsilon_{n}\underline{\psi}_{n}$ and $\epsilon_{n}\leq\epsilon_{n+1}$. For the ground state, we set $n_{f}=N$ as the Fermi level. In the topologically nontrivial phase, the Kitaev chain has a degenerate ground state with $\epsilon_{N}=\epsilon_{N+1}=0$. We can then choose $n_{f}=N-1$ or $n_{f}=N+1$, depending on whether the edge zero-mode is filled or not.

Appendix F Off-diagonal Floquet blocks and their importance

In order to computer the off-diagonal block $\mathbf{Q}$ in eq. (36), we first write:

[TABLE]

where $i=1,2$ and $\mathbf{Z}_{i}$ must satisfy the Lyapunov equation

[TABLE]

as is easily shown by requiring that $[\mathbf{A}_{i},\exp(\mathbf{A}_{i}t_{i})]=0$. A formal solution can be obtained iteratively, but this is not very illuminating.

Now, by simple matrix multiplication, $\mathbf{Q}$ takes the form:

[TABLE]

While this expression seems intractable, it can be evaluated numerically. The importance of this matrix becomes clear when we compute the stroboscopic time evolution of the covariance matrix. We start with eq. (58), which also holds when $\mathbf{X}$ is time-dependent. This equation can be rewritten using the structure matrix $\mathbf{A}(t)$ as

[TABLE]

which in turn is satisfied by the following form:

[TABLE]

We can now find a solution for the Floquet steady state by requiring that $\mathbf{D}$ is unchanged after one period:

[TABLE]

where eq. (36) was used to obtain a discrete-time Lyapunov equation for the Floquet steady state covariance matrix $\mathbf{C}_{F}$ [8].