Topological excitations at time vortices in periodically driven systems

Periodically driven quantum many-body systems are known to feature a richer topological structure than their non-driven counterparts^{1,2,3,4,5,6,7,8,9,10,11,12}. This additional structure arises from the periodicity of the quasi-energy spectrum characterizing the Floquet states. For example, periodically driven systems support “anomalous” topological phases characterized by the presence of chiral edge states, despite all the Floquet bulk bands being topologically trivial³. Systems with particle-hole symmetry (such as Bogoliubov excitations in a superconductor) can support topologically protected Majorana excitations whose quasi-energy is either zero or π/T (known as 0 and π Majorana modes, respectively), where T is the driving period^2,13. These Majorana modes may appear at defects, such as domain walls, edges, or vortices in a superconductor. Systems combining 0 and π Majoranas open new possibilities for quantum information processing^14,15,16,17, not available in static systems. The emergence of π Majorana modes is closely related to discrete-time crystallinity^18,19,20.

Periodically driven systems may possess space-time symmetries that refine the classification of topological phases^{21,22,23,24,25}, as well as novel kinds of space-time topological defects^26,27. In two-dimensional periodically driven systems, it is natural to define a space-time topological defect, a “time vortex,” which is our focus in this work. A time vortex is a point in space around which the driving Hamiltonian’s phase lag changes as a function of position, returning to itself modulo an integer multiple of 2π upon encircling the vortex (Fig. 1a). We find that a time vortex can bind topological excitations in a system of fermions with particle-hole symmetry and no time-reversal symmetry (symmetry class D²⁸). Such a system may host chiral Majorana edge modes crossing the gaps at quasi-energy 0 and π/T². We show that a time vortex binds a Majorana π mode at its core if the number of chiral edge states through quasi-energy π/T is odd.

As a possible realization of this phenomenon, we consider a periodically driven version^29,30 of Kitaev’s spin model on the honeycomb lattice³¹. Time-reversal symmetry is broken by the chirality of the driving protocol. By time-reversal symmetry, we mean that there is an anti-unitary operator ({mathcal{T}}) such that the time-dependent Hamiltonian satisfies (H(t)={mathcal{T}}H(T-t){{mathcal{T}}}^{-1}), for an appropriate choice of the origin of time. As in Kitaev’s original model, the system can be described in terms of emergent fermions with particle-hole symmetry, coupled to a static background ({{mathbb{Z}}}_{2}) gauge field. For a certain parameter range, the model realizes an anomalous Floquet phase with chiral Majorana edge states traversing the gaps around quasienergy 0 and π/T; in this phase, a ({{mathbb{Z}}}_{2}) flux carries both Majorana 0 and π modes. We demonstrate explicitly that a time vortex in this system binds a π Majorana mode, but not a 0 mode. Remarkably, this implies that by combining ({{mathbb{Z}}}_{2}) fluxes and time vortices, one can realize any combination of localized 0 and π Majorana modes. The anomalous phase with or without a time vortex can be realized simply by periodically applying a sequence of Clifford gates to a 2D array of qubits, making it particularly suitable for implementation on near-term quantum hardware.

The Floquet time vortex

To define a time vortex, we consider a reference time-periodic tight-binding Hamiltonian H₀, with period T, defined on a 2D lattice:

where r, ({{boldsymbol{r}}}^{{prime} }) denotes the sites of the system. To introduce a time vortex into H₀, we define a modified Hamiltonian H_v with a spatially dependent time delay:

where (bar{{boldsymbol{r}}}=frac{1}{2}({boldsymbol{r}}+{{boldsymbol{r}}}^{{prime} })). To describe a time vortex at the origin, we take (theta (bar{{boldsymbol{r}}})) to be given by

where we measure positive angles in the counterclockwise direction from the vector (hat{{boldsymbol{x}}}). Thus, compared to H₀, H_v has a position-dependent delay in its time dependence, which smoothly accumulates a full period T when going in a closed loop around the center of the time vortex. This configuration realizes a topological defect that cannot be smoothly deformed to a trivial spatial dependence. Fig. 1a illustrates this space-time dependence.

Being a topological defect, it is natural to ask if the time vortex can carry topologically protected modes in its core. Here, we focus on systems in symmetry class D, such that the time-dependent Hamiltonian H₀(t) itself is particle-hole symmetric at each time t, satisfying ({{mathcal{P}}}^{-1}{H}_{0}(t){mathcal{P}}=-{H}_{0}(t)) and ({{mathcal{P}}}^{2}=1), where ({mathcal{P}}) is an anti-unitary operator. It follows that H_v(t) is also in class D. A Floquet system in class D is characterized by a set of integer-valued topological invariants ({{mathcal{W}}}_{varepsilon }) defined for each gap in the quasi-energy spectrum (ε denotes a quasi-energy within a gap). The invariant ({{mathcal{W}}}_{varepsilon }) is equal to the net number of chiral edge states that appear within the corresponding gap, weighted by their chirality³ (see Fig. 1b, c). Our main result relates the properties of a time vortex to the bulk invariants of the Floquet phase: we show that if ({{mathcal{W}}}_{pi /T}) is an odd integer, then a time vortex carries a topologically protected Majorana mode with quasi-energy π/T.

Semiclassical analysis of the π mode on the edge

We now formulate a semiclassical argument showing that inserting a time vortex through a system in a cylindrical geometry can toggle the existence of π modes at both edges on and off. To this end, we consider a strip of size L_x × L_y with periodic boundary conditions along the x direction [see Fig. 1b].

We start with a tight-binding Floquet model of fermions in class D, with translation invariance in the x direction (i.e., without a time vortex):

where a₀ is the lattice spacing. We assume that the Hamiltonian has a finite range ℓ, such that the matrix element (4) vanishes for (| {boldsymbol{r}}-{{boldsymbol{r}}}^{{prime} }|, >, ell).

For constructing the Floquet states in the presence of a time vortex, it will be useful to consider the Floquet states of this translation-invariant problem with twisted boundary conditions parameterized by a phase α. In accordance with Bloch’s theorem, we label the Floquet states of the system by a crystal momentum component k along (hat{{boldsymbol{x}}}) (around the cylinder) and band index n, along with the twist angle parameter α:

with the twisted boundary conditions

The corresponding quasi-energies ε_nk(α) follow from the application of Floquet’s theorem,

The twisted boundary conditions yield a discrete set of allowed crystal momentum values

with (m=0,1,ldots ,frac{{L}_{x}}{{a}_{0}}-1). While only α = 0, π are allowed for a class D Hamiltonian, here we consider general α as will be needed for the construction below.

We now consider the Floquet states in the presence of a time vortex threaded through the holes of the cylinder [as opposed to a point-like time vortex within the bulk of the system as described in Eq. (3)]. A time vortex can be threaded through the holes of the cylinder by setting (theta (bar{{boldsymbol{r}}})=2pi bar{x}/{L}_{x}) in Eq. (2). We will show that threading such a time vortex through a hole of the cylinder toggles the presence of the π mode in that hole on and off. In turn, we will argue that this behavior implies that a point-like time vortex within the bulk of the 2D system must carry a π mode.

In the section “Approximate Floquet states with a time vortex” we show that the Floquet eigenstate in the presence of the time vortex can be written as

for small ℓ/L_x, the Floquet states in the presence of the time vortex follow those of the defect-free system, with the local time lag of the Hamiltonian H_v imprinted directly onto χ_nk(r, t). Here, k satisfies Eq. (8). Note that in the presence of the time vortex, the Hamiltonian is not translationally invariant, but has a combined space-time translation symmetry. The label (k) of (chi_{nk}) is associated with the eigenvalue of the Floquet eigenstate under this combined space-time translation.

Next, we fix the value of α by imposing the physically desired periodic boundary conditions on χ_nk:

which is satisfied by setting α to ({alpha }_{n,m}^{star }), defined as the solution of the equation ({alpha }_{n,m}^{star }=-{varepsilon }_{nk}({alpha }_{n,m}^{star })T) with (k=(2pi m+{alpha }_{n,m}^{star })/{L}_{x}). This relation must be solved separately for each n and m. Therefore, comparing with Eq. (8), we find that to leading order in (frac{ell }{{L}_{x}}), the allowed values of k shift as

whereas the momentum in the model without the time vortex with periodic boundary conditions is discretized as (k=frac{2pi m}{{L}_{x}}). In particular, for states near quasi-energy π/T, the insertion of the time vortex effectively changes the boundary conditions from periodic to anti-periodic. In contrast, the boundary conditions for states near quasi-energy ε_nk = 0 do not change.

Assuming that there is a chiral mode through quasi-energy π/T, we can now show that for a finite system with circumference L_x, introducing a time vortex through the cylinder toggles between having a Majorana π mode at each edge and not having one. Linearizing the quasi-energy spectrum near ε = π/T gives ({varepsilon }_{nk}approx frac{pi }{T}+v(k-{k}_{0})), where particle-hole symmetry maps ε_k to −ε_−k, and hence constrains k₀ to be either 0 or π/a₀. Then, using Eq. (8) with α = 0 and Eq. (11) to describe systems with n_v = 0 or 1-time vortices, respectively, we find that the quasi-energy spectrum for a finite system, to order ℓ/L_x, is given by

Using Eq. (12), we can verify that if the system has a Majorana π mode at the edge (whose quasi-energy is exactly π/T) in the absence of a time vortex (n_v = 0), then the system with n_v = 1 does not have π mode, and vice versa. In particular, if k₀ = 0, the system with n_v = 0 has a π mode at the edges and the n_v = 1 system does not. If k₀ = π/a₀, the n_v = 0 system has a π mode if the number of sites around the cylinder, (L_x/a_0), is even, and does not have one if L_x/a₀ is odd. In the n_v = 1 case, the situation is reversed.

This argument implies that a point-like time vortex within the bulk must carry a π mode. To see this, consider a cylinder with a single time vortex piercing through its surface. Then, there must also be a time vortex through one of the holes of the cylinder but not the other. According to Eq. (12), adding a time vortex through a hole toggles between having a π mode or not at the corresponding edge. Therefore, inserting a time vortex through one of the holes of the cylinder and extracting it through the bulk changes the parity of the number of π modes on the edges from even to odd. Since the total number of Majorana π modes in the entire system must be even, the other Marjorana π mode has to be located in the bulk; this π mode is localized around the core of the bulk time vortex, which is the only point where the bulk gap around quasi-energy π/T can close.

Periodically driven Kitaev honeycomb model

As an application of these ideas, we use a periodically driven Kitaev spin model on the honeycomb lattice³⁰, which maps to a problem of free fermions in class D (Majorana fermions) coupled to a static background ({{mathbb{Z}}}_{2}) gauge field³¹. Introducing a time vortex enriches the types of topological excitations that are accessible in this system. In addition, the spin model and the associated time vortices can be realized on near-term digital quantum computers.

The model consists of spin-(frac{1}{2}) degrees of freedom on the vertices of a honeycomb lattice with time-dependent anisotropic exchange interactions. The Hamiltonian is

Here, I = (i, j, s) denotes the sites of the honeycomb lattice, where i ∈ {1, …, N_x} and j ∈ {1, …, N_y} label a unit cell, and s = A, B labels the two sublattices. Each site hosts a spin σ_I. The bonds of the lattice are partitioned into three sets a = {x, y, z} according to the three types of bond orientations in the honeycomb lattice. We take ({{mathcal{J}}}_{{bf{I}},{bf{J}}}^{a = x,y,z}) to be periodically time-dependent with a period T [see Fig. 2a].

This model can be elegantly solved with arbitrary coupling strengths by fermionization as is done in Ref. ³¹. We introduce a set of four Majorana operators ({{b}_{{bf{I}}}^{x},{b}_{{bf{I}}}^{y},{b}_{{bf{I}}}^{z},{c}_{{bf{I}}}}) for each site I, and a constraint ({b}_{{bf{I}}}^{x}{b}_{{bf{I}}}^{y}{b}_{{bf{I}}}^{z}{c}_{{bf{I}}}=1), which results in a two-level system for each site. We map the spins to fermions via ({sigma }_{{bf{I}}}^{a}=i{b}_{{bf{I}}}^{a}{c}_{{bf{I}}}) and identify a conserved ({{mathbb{Z}}}_{2}) gauge field defined on the bonds of the lattice,

This leads to the Hamiltonian

The variables u_I,J are not gauge-invariant: applying the gauge operator ({b}_{{bf{I}}}^{x}{b}_{{bf{I}}}^{y}{b}_{{bf{I}}}^{z}{c}_{{bf{I}}}) changes the sign of u_I,J on the three bonds touching vertex I. However, the ({{mathbb{Z}}}_{2}) flux W_p given by the product of the gauge fields around a hexagonal plaquette p is gauge-invariant. Note also that the free fermion Hamiltonian (15) is purely imaginary, and is hence particle-hole symmetric with ({mathcal{P}}={mathbb{I}}).

Using the fact that the operators {u_I,J} are conserved by the Hamiltonian, we fix a gauge and replace them with their eigenvalues, u_I,J = ±1. With this substitution, the Hamiltonian (15) describes free fermions, {c_I}.

We start from a translation-invariant model given by ({{mathcal{J}}}_{{bf{I}},{bf{J}}}^{a}(t)={tilde{{mathcal{J}}}}^{a}(t)), where we choose each ({tilde{{mathcal{J}}}}^{a}(t)) to be a train of rectangular pulses of width Δt and amplitude ({{mathcal{J}}}_{0}^{a}), starting at time lT + t_a and ending at lT + t_a + Δt, for (lin {mathbb{Z}}). We take t_x = 0, ({t}_{y}=frac{1}{3}T), ({t}_{z}=frac{2}{3}T). The Kitaev honeycomb model is illustrated in Fig. 2.

In order to insert a time vortex at r₀, we follow the prescription in the section “The Floquet time vortex” by adding a delay as a function of the position of the midpoint of each bond ({bar{{boldsymbol{r}}}}_{{bf{I}},{bf{J}}}=frac{1}{2}left({boldsymbol{r}}({bf{I}})+{boldsymbol{r}}({bf{J}})right)) such that

In a time-reversal symmetry-broken version of the static Kitaev model, a phase with a non-zero Chern number of the c_I fermions can occur in which ({{mathbb{Z}}}_{2}) fluxes (where W_p = −1) bind Majorana zero modes. In contrast, Ref. ³⁰ showed that in the anomalous Floquet phase, realized in the vicinity of ({{mathcal{J}}}_{0}^{a}Delta t=pi /4), such ({{mathbb{Z}}}_{2}) fluxes carry Majorana modes at both 0 and π/T quasi-energies. Section “Phase diagram of the Floquet Kitaev honeycomb model” presents the phase diagram of the Floquet Kitaev honeycomb model.

The limit of short pulses with no overlaps, Δt < T/3, at “resonance,” ({{mathcal{J}}}_{0}^{a}Delta t=pi /4), turns out to be very simple to analyze. Each pulse serves to swap the Majorana operators connected by the bond on which it acts. After the action on the bond 〈I, J〉 with I in the A sublattice and J in the B sublattice, c_I → c_J and c_J → −c_I. In the spatially uniform model, the order of the pulses dictates that a full period of time-evolution effectively swaps a pair of Majoranas located on opposite corners of each plaquette in the bulk, c_{(i, j, A)} → −c_{(i+1, j−1, B)} and c_{(i+1, j−1, B)} → c_{(i, j, A)}. The result is that Floquet eigenmodes are split across pairs of sites as (left[{c}_{(i,j,A)}pm i{c}_{(i+1,j-1,B)}right]/sqrt{2}), with quasi-energies (pm frac{pi }{2T}). Creating a ({{mathbb{Z}}}_{2}) flux on a plaquette adds a minus to the exchange of the two fermions, altering the Floquet eigenmodes to real combinations of the form (left[{c}_{(i,j,A)}pm {c}_{(i+1,j-1,B)}right]/sqrt{2}), with quasi-energies (varepsilon =0,frac{pi }{T}). In terms of the spin degrees of freedom, at the resonant limit each pulse simply implements a Clifford gate on neighboring spins (U=exp (ifrac{pi }{4}{sigma }_{{bf{I}}}^{a}{sigma }_{{bf{J}}}^{a})), where we note that an operator of the form (exp (ifrac{pi }{4}P)) is a Clifford unitary for any (multi-qubit) Pauli operator P.

Time vortex

Now we consider inserting a time vortex at the center of a plaquette in a system with no ({{mathbb{Z}}}_{2}) fluxes. As described above, the time vortex adds delays to the time dependence of the couplings of each bond relative to the uniform system. The simplicity of the analysis above for the uniform case with short pulses was enabled by the fact that the exchange pulses on neighboring bonds never acted simultaneously (allowing disconnected pairs of sites to evolve within each pulse). This simplicity can be retained for the analysis of a situation with a time vortex by taking the pulses to be sufficiently short, with Δt < T/6, such that the pulses acting on the bonds of the plaquette hosting the time vortex core are sufficiently separated in time that no overlaps are created. We find it furthermore convenient to take Δt ≲ T/N_x, Δt ≲ T/N_y, so that at the initial time t = 0, none of the pulses anywhere on the lattice are active. In this “instantaneous pulse” limit, the effect of the time vortex is only a rearrangement of the order of the pulses acting on the bonds of the lattice. The order of the pulses acting on the bonds and the action of the full period evolution operator on the Majoranas are illustrated in Fig. 3. Below, we will demonstrate with numerical simulations that the qualitative results we obtain from this simple limit persist under more generic driving conditions, including the case Δt > T/3 where neighboring pulses overlap.

It is easy to see that in a region far away from the vortex core, the order of pulses acting on the bonds will locally be the same as in the absence of the vortex. However, depending on the orientation of the region relative to the vortex core, the order may change by a cyclical permutation to become y → z → x or z → x → y. In any case, plaquettes in such regions will host localized Floquet modes with quasi-energies (pm frac{pi }{2T}). Thus the bulk spectrum remains undisturbed by the insertion of the time vortex. The effect of the time vortex at the vortex core is also easy to understand. The order of the pulses on the bonds around the plaquette hosting the time vortex core becomes a cycle going around the plaquette in the clockwise direction, starting from the z pulse on the left bond of the plaquette (note that all six pulses around the core act within one driving period T). Therefore, the Majorana at the lower left corner of the plaquette c_{(i, j, B)} will hop entirely around the plaquette in a single period. To determine the quasi-energy of the corresponding localized mode, we note that each pulse implements a transformation c_A → c_B, c_B → −c_A. Following the arrows in Fig. 2a, three factors of (−1) are incurred over the entire period. Consequently, the overall phase is π, and we find that this eigenstate is a Majorana π mode bound to the time vortex.

π modes in a finite-size simulation

After confirming the existence of the π mode bound to the time vortex in the resonant, instantaneous pulse limit, we now consider the generic situation where the pulses do not satisfy the resonance condition and pulses on neighboring bonds overlap in time. We demonstrate that the π mode survives and remains localized around the vortex core. The phase diagram of the Floquet Kitaev honeycomb model highlighting the region in parameter space in which a time vortex binds a π mode is shown in section “Phase diagram of the Floquet Kitaev honeycomb model”.

To treat the generic situation, we numerically perform time-evolution of the Majorana operators c_I in the Heisenberg picture. We then calculate the spatially resolved time-averaged density of states^12,32

Here ({widetilde{G}}_{{bf{I}},{bf{J}}}^{R}(omega )) is the Fourier transform of the time-averaged single-particle retarded Green’s function given by

where the Green’s function itself is

Here, since the c_I fermions are free, the expectation value is independent of the state it is taken in, as long as there are no ({{mathbb{Z}}}_{2}) fluxes in the system.

Figure 4a shows the time-averaged density of states summed over the sites of the central plaquette (which hosts the time vortex), in a system of size 16 × 16 unit cells. The parameters of the simulation are given in the figure caption. In this generic regime, the flat bulk bands at quasi-energy ±π/(2T) discussed above for the resonant limit acquire a finite width, as expected. Importantly, the gaps around quasi-energy 0 and π/T remain finite. Consequently, the π Majorana mode at the vortex core remains robust, as can be seen by the sharp peak at ω = π/T. We also find an additional (unprotected) bound state at the vortex core, near the top of the upper band and the bottom of the lower band. Figure 4b shows the spatially resolved time-averaged density of states, integrated within the range ω ∈ [0.99π/T, 1.01π/T], in a system of size 8 × 8 unit cells. The spectral function has support near the center, where the time vortex is located (the plaquette indicated by p in the figure), and at the edge of the system. The π mode at the vortex core hybridizes with the π mode at the edge, but this hybridization is expected to decrease exponentially with the system size.

Topological invariant

Topological defects in periodically driven non-interacting systems of different symmetry classes have been classified in Ref. ²⁷. In particular, it was found that a zero-dimensional (point-like) defect in a two-dimensional periodically driven system is characterized by a ({{mathbb{Z}}}_{2}times {{mathbb{Z}}}_{2}) invariant. The two ({{mathbb{Z}}}_{2}) invariants, ν₀ and ν_π, are equal to the number of Majorana 0 and π modes bound to the center of the defect, respectively. Our goal is to compute ν₀ and ν_π for the time vortex; in particular, we show that in the anomalous phase of the driven Kitaev honeycomb model of section “Periodically driven Kitaev honeycomb model”, ν₀ = 0 and ν_π = 1.

We consider the Hamiltonian acting on a region far from the time vortex core, at an angle θ around the core [as defined in Eq. (3)]. The Hamiltonian is a slowly varying function of the position (bar{{boldsymbol{r}}}), so by treating (bar{{boldsymbol{r}}}) as a constant in a given region we can define the Bloch Hamiltonian H₀(k, t; θ)³³, obtained by taking θ in the right-hand side of Eq. (2) to be constant and Fourier transforming with respect to ({boldsymbol{r}}-{{boldsymbol{r}}}^{{prime} }). As defined, H₀(k, t; θ) is periodic both in t and in θ. A construction of the ({{mathbb{Z}}}_{2}times {{mathbb{Z}}}_{2}) invariant for H₀(k, t; θ) was given in Ref. ²⁷; here, we give an alternative formulation of the invariant, which is more directly related to the number of Majorana modes at the defect.

To this end, we construct a family of time-dependent Hamiltonians parameterized by λ ∈ [0, 1], which smoothly interpolates between the vortex-free model, obtained for λ = 0 (with trivial topological invariants) and the model with the vortex for λ = 1:

This family of Hamiltonians is not periodic in θ for λ ∉ {0, 1}. The non-periodicity has a simple consequence in real space: the position-dependent Hamiltonian around the time vortex changes abruptly along the ray θ = 0 unless λ = 0 or 1 (Fig. 5a). One-dimensional states that are bound to this ray may close the Floquet gaps around quasi-energy 0 and π/T at some values of λ. Such a gap closing and reopening, corresponding to a topological phase transition along the θ = 0 ray, is necessary in order to obtain a Majorana mode at the core of the defect (the edge of the θ = 0 ray) for λ = 1.

Thus, we have reduced the problem to determine the topological properties of a one-dimensional system along the θ = 0 ray in real space. Far away from the vortex core, the θ = 0 ray can be viewed as a domain wall between a y > 0 region where the phase of the driving Hamiltonian is 0 and a y < 0 region where the phase is 2πλ (see Fig. 5a). The system with such a domain wall is translationally invariant along the θ = 0 ray, and hence k_x is conserved.

As is often the case for systems in class D, to probe for topological gap closings and re-openings in the θ = 0 ray as a function of λ, it suffices to consider the two particle-hole symmetric momenta, k_x = 0 and π/a₀³⁴. This can be done by studying ({H}_{0}({k}_{x}in {0,frac{pi }{{a}_{0}}},{k}_{y},t;2pi lambda )), which describes the region below the domain wall (y < 0). For fixed values of k_x, this Hamiltonian describes a family of one-dimensional (along k_y) Floquet Hamiltonians parametrized by λ. Each such family of 1D Hamiltonians is itself characterized by a pair of ({{mathbb{Z}}}_{2}) invariants²⁷, one for each of the particle-hole symmetric quasi-energies, ε = 0 and π/T.

Physically, ({H}_{0}({k}_{x}in {0,frac{pi }{{a}_{0}}},{k}_{y},t;2pi lambda )) defines a Floquet version of an “adiabatic fermion parity pump”^33,35. We denote the resulting ({{mathbb{Z}}}_{2}) invariants as ν(ε, k_x), corresponding to the fermion parity pumped by the system labeled by ({k}_{x}in {0,frac{pi }{{a}_{0}}}) at quasi-energy (varepsilon in {0,frac{pi }{T}}) as the parameter λ changes adiabatically from 0 to 1. We identify the topological invariants of the time vortex as ({nu }_{0}=nu (0,0)+nu (0,frac{pi }{{a}_{0}}),(mathrm{mod},,2)) and ({nu }_{pi }=nu (frac{pi }{T},0)+nu (frac{pi }{T},frac{pi }{{a}_{0}}),(mathrm{mod},,2)). If ν(ε, k_x) = 1, then the gap at the end of the system has to close at quasi-energy ε for some value of λ. This, in turn, implies a topological gap closing along the ray at the same k_x at this value of λ (i.e., a level crossing of two states related by particle-hole symmetry). Viewed as an adiabatic process, such a level crossing changes the fermion parity corresponding to the Dirac fermion formed by the two Majorana states that cross, hence it is referred to as a fermion parity pump.

The invariants ν(ε, k_x) can be determined as follows. We examine the evolution operator from time 0 to t,

in the three-dimensional space of ({k}_{y}in [0,frac{pi }{{a}_{0}}]), λ ∈ [0, 1), and t ∈ [0, T), for fixed k_x = 0 or (frac{pi }{{a}_{0}}) (Fig. 5b). We denote the phases of the eigenvalues of U (known as the phase bands⁴) by ϕ_n(k, t; 2πλ), where n is the band index. Along the plane t = 0, we have (U={mathbb{I}}). At t = T, U(k, T; 2πλ) is equal to the Floquet operator, whose spectrum is assumed to be gapped around ε = 0 and π/T.

For every plane of fixed time t > 0 for which the spectrum of (U({k}_{x}in {0,frac{pi }{{a}_{0}}},{k}_{y},t;2pi lambda )) is gapped (i.e., the phase bands do not cross 0 and π), a ({{mathbb{Z}}}_{2}) invariant can be defined as follows. We consider the effective Hamiltonian

for fixed ({k}_{x}in {0,frac{pi }{{a}_{0}}}), as a function of k_y and λ (since H_eff(k, t; 2πλ) is assumed to be gapped at (varepsilon =frac{pi }{T}) for this value of t, we can use the usual convention for the branch cuts of (log (z)), with a discontinuity along Re(z) < 0, ({rm{Im}}(z)=0)). This can be viewed as a time-independent 1D Hamiltonian in class D that depends periodically on a parameter λ; the ({{mathbb{Z}}}_{2}) invariant is equal to the fermion parity pumping over a single cycle of changing λ adiabatically^33,35. For an explicit prescription for the computation of this invariant, see section “Computation of the topological invariants”.

The sum of the invariants at quasi-energies 0 and π/T, (nu (0,{k}_{x})+nu (pi /T,{k}_{x}),(mathrm{mod},,2)), is simply equal to the fermion parity pumping invariant of H_eff (k, T; 2πλ) as a function of k_y and λ for fixed ({k}_{x}in {0,frac{pi }{{a}_{0}}}). Disentangling ν(0, k_x) and (nu (frac{pi }{T},{k}_{x})) is more subtle, and depends on the evolution operator at intermediate times, 0 < t < T. To this end, we note that ν(π/T, k_x) can only change via a topological gap closing of H_eff(k, t, 2πλ) at ε = π/T for some intermediate t. We therefore track the topological singularities where the gap of H_eff(k, t; 2πλ) closes at quasi-energy π/T⁴, and compute the ({{mathbb{Z}}}_{2}) fermion parity pumping invariants of H_eff(k, t; 2πλ) on constant t planes before and after the singularities, illustrated in Fig. 5b. For each such pair of planes, we define an index ζ_i equal to the difference (mod 2) of the ({{mathbb{Z}}}_{2}) invariants of the two planes, that captures the change in the invariant across the singularity. The sum (mod 2) of ζ_i’s across all (varepsilon =frac{pi }{T}) topological singularities for 0 < t < T gives (nu (frac{pi }{T},{k}_{x})).

In practice, the topological singularities may occur on higher-dimensional regions (lines, surfaces, or volumes) in the (k_y, t, λ) space, rather than at isolated points (as assumed, for simplicity, in Fig. 5b). In that case, there may not be planar surfaces that separate the ε = 0 and (frac{pi }{T}) singularities. To handle this situation, U(k, t; 2πλ) can in principle be deformed continuously without closing the gap at t = T such that there are planes of constant t that separate all the topological singularities. Alternatively, one can find non-planar surfaces that separate the singularities, on which the invariants of H_eff are computed; these surfaces are required to be periodic in λ for the invariant to be well-defined.

Using this procedure, we find the invariants of the time vortex in the anomalous phase of the Kitaev honeycomb model to be ν₀ = 0, ν_π = 1, as expected. The details of the computation are provided in the section “Computation of the topological invariants”. The topological invariants of the two types of topological defects in this phase, the ({{mathbb{Z}}}_{2}) flux and the time vortex, are summarized in Table 1.

To summarize, in this work, we have considered topological defects in periodically driven systems that utilize their intrinsic symmetry under translation in time. A spatially uniform shift in the origin of time does not affect the Floquet spectrum and changes the Floquet states by a unitary transformation. In contrast, a space-dependent time delay of the drive can have physical consequences. Intuitively, a time vortex, where this time delay winds by a full period, has a simple effect on the Floquet states far away from the center of the defect: when going around the defect, the wavefunction acquires a phase that is proportional to the quasi-energy of the Floquet state times T. This phase can be thought of as arising from the additional driving period “experienced” by the wavefunction upon encircling the defect. As we have shown, in class D, the change of phase implies that a time vortex binds a Majorana π mode. More generally, a time vortex defect with n_v-fold winding binds a Majorana π mode if n_v is odd.

We note that the appearance of topological π modes at time vortices is not unique to symmetry class D. Similar phenomena are expected to occur in any symmetry class whose topological classification in d = 2 is non-trivial (such that topological edge states exist), and that includes particle-hole symmetry (such that π modes can be topologically protected). In such cases, if the gap at quasi-energy π/T is non-trivial, a time vortex may bind π modes. In addition to class D, these conditions are met also in classes DIII and C of the Altland-Zirnbauer classification of topological insulators and superconductors^4,27.

In terms of practical utility, this construction may be useful for topological quantum computation as it allows manipulating logical qubits encoded in zero and π modes independently. Using both types of Majorana modes leads to a more efficient encoding rate while still protecting the logical qubits from local noise, as long as the noise varies slowly in the time since the logical operations involve pairs of Majorana operators separated in space or in quasienergy.

The time vortex is an externally imposed topological point defect and can be toggled and manipulated by the driving Hamiltonian. In contrast, the ({{mathbb{Z}}}_{2}) flux is an intrinsic excitation that depends on the quantum state rather than on the driving Hamiltonian, perhaps making it more challenging to control. In fact, as demonstrated in the section Periodically driven Kitaev honeycomb model—Time vortex, the time vortex can be implemented in a very simple spin model requiring only two-body local Clifford gates, prepared in the ({{mathbb{Z}}}_{2}) flux-free sector (topologically equivalent to the ground state of the toric code). The time vortex can be toggled on and off simply by rearranging the order of the periodic gate sequence. This makes the time vortex particularly suitable for implementation on quantum devices which allow for a high accuracy application of Clifford gates.

Approximate Floquet states with a time vortex

In this section, we explicitly verify that the ansatz for the Floquet states on the cylinder in the presence of the time vortex in Eq. (9) satisfies the Schrödinger equation to leading order in (frac{ell }{{L}_{x}}). By substitution we find

where we have used the definition of H_v(t) in Eq. (2), with (theta (bar{{boldsymbol{r}}})=2pi bar{x}/{L}_{x}) as discussed above Eq. (9), and suppressed the spatial indices r and ({{boldsymbol{r}}}^{{prime} }) for notational simplicity, invoking a matrix notation, e.g., ({H}_{{rm{v}}}(t){chi }_{nk}(t)equiv {sum }_{{{boldsymbol{r}}}^{{prime} }}{H}_{{rm{v}}}({boldsymbol{r}},{{boldsymbol{r}}}^{{prime} },t){chi }_{nk}({{boldsymbol{r}}}^{{prime} },t)).

Our aim is to show that the right-hand side of Eq. (23) is of the order of (frac{ell }{{L}_{x}}). Performing a position-dependent shift of the time variable in Eq. (23), (tmapsto t+frac{x}{{L}_{x}}T), we find:

To obtain the second line, we expanded H₀ around (t-frac{{x}^{{prime} }-x}{{L}_{x}}T), and used (| {x}^{{prime} }-x|, <, ell), where ℓ is the range of the hopping in H₀. In obtaining the final expression, we used the Schrödinger equation for ψ_nk(t) to eliminate the first term in the second line. (Notice that for any fixed, continuous driving profile, the first-order term in the Taylor expansion can be made dominant over the higher-order terms by choosing the circumference of the cylinder, L_x, large enough.)

In addition to solving the Schrödinger equation, using Eq. (7) and Eq. (9), we find that χ_nk(t) satisfies the Floquet boundary condition in time: ({chi }_{nk}(t+T)={e}^{-i{varepsilon }_{nk}(alpha )T}{chi }_{nk}(t)), where α is determined as described below Eq. (10). This shows that χ_nk(t) is a Floquet eigenstate with quasi-energy ε_nk(α).

We note that, while χ_nk(r, t) is not a Bloch state, it is an eigenstate of a combined space-time translation:

which is a symmetry of the Hamiltonian in the presence of a time vortex.

Phase diagram of the Floquet Kitaev honeycomb model

In this section, we map out the phase diagram of the Floquet Kitaev honeycomb model [Eq. (15)] in the flux-free sector. This is done by searching the space of parameters (pulse duration Δt and the pulse integral ({{mathcal{J}}}_{0}^{a}Delta t)) for gap closings at quasienergy 0 and π/T, as done in Ref. ³⁰. The topological invariants ({{mathcal{W}}}_{0}), ({{mathcal{W}}}_{pi /T}) are then determined by computing the spectrum on a cylindrical geometry, and counting the number of chiral edge states at each quasi-energy gap, weighted by their chirality. The phase diagram is shown in Fig. 6, where the gap closings are indicated by black lines. The resonant point analyzed in the main text is located at ({{mathcal{J}}}_{0}^{a}Delta t=pi /4) with non-overlapping pulses (Δt < 1/3) (the red circle belongs to this region).

In the shaded region in Fig. 6 (that includes the resonant point) we find that ({{mathcal{W}}}_{pi /T}=1,mathrm{mod},,2), except on the phase boundaries where ({{mathcal{W}}}_{pi /T}) is not defined. Therefore, within this region, a time vortex defect binds a Majorana π mode, as argued in the main text. In particular, the point examined away from the resonant limit in Fig. 1c and Fig. 4 is within the same phase as the resonant point.

Kitaev honeycomb model in momentum space

The ({{mathbb{Z}}}_{2}) flux-free sector W_p = 1 can be described by choosing a gauge in which u_I,J = 1 for I in the A sublattice and J in the B sublattice. Here, taking spatially uniform couplings (i.e., in the absence of a time vortex), translation symmetry allows us to write the Hamiltonian (15) in momentum space as

where c_k,s is the Fourier transform of c_I,

Here, for convenience, we take the components k_x, k_y as the projections of the momentum k on the (non-orthogonal) primitive vectors of the honeycomb lattice. H_k(t) is given by

where the Pauli matrices τ^a act on the sublattice degree of freedom. We resolve the time-evolution operator in momentum space as U(k, t) and find that particle-hole symmetry acts as

The time-evolution operator can be expressed in terms of a spectral decomposition

where ϕ_n(k, t) are called the phase bands⁴, and are defined modulo 2π. At t = T, ε_n(k) = ϕ_n(k, T)/T is referred to as the quasienergy, and unlike the bands of a static system, these are periodic with period (frac{2pi }{T}). Particle-hole symmetry maps ϕ_n(k, t) to − ϕ_n (−k, t), such that there are two particle-hole symmetric quasi-energies, ε = 0 and (frac{pi }{T}), at which protected Majorana modes can appear.

Computation of the topological invariants

This section provides the details of the numerical calculation yielding the topological invariants of the time vortex defect in the driven Kitaev honeycomb model. We specialize to the resonant limit discussed in the section “Periodically driven Kitaev honeycomb model—Time vortex” and follow the procedure outlined in the section “Periodically driven Kitaev honeycomb model—Topological invariant”. We set the duration of the pulses during which ({tilde{{mathcal{J}}}}^{a}) are non-zero to be Δt = T/3. Since the invariant only requires knowledge of the time-evolution of a translationally invariant system, the pulses corresponding to different couplings do not overlap.

We begin by numerically mapping out the locations of the “topological singularities” in the three-dimensional space of k_y, λ, t for both k_x = 0 and (frac{pi }{{a}_{0}}) by exponentiating the Hamiltonian (28) for each pulse separately and taking products of these unitaries. For k_x = 0, we find that the ϕ = 0 singularities are located on a one-dimensional curve lying in the ({k}_{y}=frac{pi }{{a}_{0}}) plane, originating at t = 0 and extending up to (t=frac{2}{3}T). On this plane, ({H}_{{boldsymbol{k}}}=left({tilde{{mathcal{J}}}}^{x}+{tilde{{mathcal{J}}}}^{y}-{tilde{{mathcal{J}}}}^{z}right){tau }^{y}). Therefore, time-evolution on this plane corresponds to rotation around the τ^y axis of the Bloch sphere in the positive direction during the action of the x and y pulses, and in the negative direction during the action of the z pulse. The topological singularities at ϕ = 0 correspond to the cancellation of the z pulse with either the preceding y pulse or the successive x pulse, or to a cancellation of complementary parts of x and y pulses with the z pulse acting in between. Singularities at ϕ = π are found on the plane (t=frac{2}{3}T), both on the line at k_y = 0 and on the line at λ = 0. At k_x = k_y = 0 all pulses act identically, so for any λ, two thirds of a driving period leads to the accumulation of two pulses. At λ = 0, the singularity is due to the accumulation of an x and a y pulse which are equivalent at k_x = 0 for all k_y.

At ({k}_{x}=frac{pi }{{a}_{0}}), the ϕ = 0 singularities are found both on 1D curves and on a 2D plane. These are all smoothly connected to the plane t = 0 and extend up to (t=frac{2}{3}T). On this plane there is a cancellation of the x and y pulses. Since during this part of the period ({tilde{{mathcal{J}}}}^{z}=0), the plane is parallel to the k_y axis. The line segment singularities appear on the planes ({k}_{y}=0,frac{pi }{{a}_{0}}) and correspond to the cancellation of the z pulse with the y or x pulse, respectively. The singularities at ϕ = π are again confined to the (t=frac{2}{3}T) plane, but here they appear at two isolated points at ({k}_{y}=0,lambda =frac{1}{3}), and ({k}_{y}=frac{pi }{{a}_{0}},lambda =frac{2}{3}), and are due to the accumulation of z and x or y and z pulses, respectively.

Next, for k_x = 0 and ({k}_{x}=frac{pi }{{a}_{0}}) separately, we construct surfaces that are periodic in λ and separate the two types of singularities. Since both types of singularities exist on the plane (t=frac{2}{3}T), this cannot be accomplished by planes of constant t. Instead, we design suitable piece-wise planar surfaces composed of edge-sharing polygons, as shown along with the singularities in Fig. 7.

Then, as explained in the main text, we compute a ({{mathbb{Z}}}_{2}) topological invariant on each one of these surfaces, and on the t = T plane. On each surface, we define the effective gapped Hamiltonian H_eff(k, t; 2πλ) [Eq. (22)]. We describe the procedure for computing the invariant on a plane of constant t; the generalization to any surface characterized by a function t = f(k_y, λ) which is periodic in λ is straightforward.

The ({{mathbb{Z}}}_{2}) invariant characterizing the surface is the fermion parity pump invariant³³ of the effective one-dimensional static Hamiltonian defined by H_eff as a function of k_y, when the parameter λ changes adiabatically from 0 to 1. To compute it, we follow the procedure described in Ref. ³⁵. We first diagonalize the Hamiltonian H_eff for each k_y and λ, denoting the diagonalizing matrix by S(k_y, λ), where each row of S is an eigenvector of H_eff. Importantly, it is possible to find a continuous gauge choice for the rows of S(k_y, λ) in the domain ({k}_{y}in [0,frac{pi }{{a}_{0}}]), λ ∈ [0, 1), because S(k_y, λ) is not required to be periodic in k_y within this domain (which spans only part of the Brillouin zone).

For k_x and k_y both equal to 0 or (frac{pi }{{a}_{0}}), H_eff(k, t; 2πλ) is particle-hole symmetric, and can be chosen to be purely imaginary, ({H}_{{rm{eff}}}=-{H}_{{rm{eff}}}^{* }) [as is manifestly the case for the driven Kitaev model, Eq. (28)]. Then, the eigenvectors of H_eff come in pairs, ({| uleft.rightrangle ,| {u}^{* }left.rightrangle }), with eigenvalues {ε, − ε}, where (| {u}^{* }left.rightrangle) is the complex conjugate of (| uleft.rightrangle). We define a matrix (O({k}_{y}in {0,frac{pi }{{a}_{0}}},lambda )) whose rows are formed of the real and imaginary parts of the rows of (S({k}_{y}in {0,frac{pi }{{a}_{0}}},lambda )): the corresponding rows of O are ({frac{1}{sqrt{2}}(| uleft.rightrangle +| {u}^{* }left.rightrangle ),frac{1}{isqrt{2}}(| uleft.rightrangle -| {u}^{* }left.rightrangle )}). By the hermiticity of H_eff, (| uleft.rightrangle) and (| {u}^{* }left.rightrangle) are orthogonal, and hence O is an orthogonal matrix. Since ({pi }_{1}[SO(N)]={{mathbb{Z}}}_{2}), two ({{mathbb{Z}}}_{2}) invariants can be defined per ({k}_{x}in {0,frac{pi }{{a}_{0}}}), corresponding to O(k_x, k_y = 0, λ) and (O({k}_{x},{k}_{y}=frac{pi }{{a}_{0}},lambda )).

To compute the invariants in practice, we note that as explained in Ref. ³⁵, the ({{mathbb{Z}}}_{2}) invariants are equal to the parity of the number of times the phases of the eigenvalues of (O({k}_{x},{k}_{y}in {0,frac{pi }{{a}_{0}}},lambda )) cross π modulo 2π as λ varies from 0 to 1. The invariants corresponding to k_y = 0 and (frac{pi }{{a}_{0}}) depend on the gauge choice, but their sum (mod 2) is gauge-invariant, because of the requirement that S(k_x, k_y, λ) is continuous in the entire domain ({k}_{y}in [0,frac{pi }{{a}_{0}}]), and is equal to the ({{mathbb{Z}}}_{2}) invariant characterizing H_eff on the surface t = f(k_y, λ). The values of the invariants on the surfaces separating the phase singularities and on t = T are shown in Table 2.

Finally, as discussed in the main text, we identify (nu (varepsilon =0,{k}_{x})+nu (varepsilon =frac{pi }{T},{k}_{x}),(mathrm{mod},,2)) as the two ({{mathbb{Z}}}_{2}) invariants computed on the t = T planes. Adding the two invariants for k_x = 0 and ({k}_{x}=frac{pi }{{a}_{0}}) gives ({nu }_{0}+{nu }_{pi },(mathrm{mod},,2)). To disentangle ν₀ and ν_π, we subtract (mod 2) the invariant on the t = T plane from the invariant of the surface that separates the ϕ = π phase singularities from the ϕ = 0 phase singularities. For the Kitaev model example, these surfaces are shown in Fig. 7 for k_x = 0 and (frac{pi }{{a}_{0}}). As explained in the main text, the result is the change in the ({{mathbb{Z}}}_{2}) invariant that corresponds to the ϕ = π gap, since the invariant can only change by a phase singularity (gap closing). The change in the invariant across the phase singularities gives ν_π, because the ϕ = π invariant is trivial on the plane t = 0. We thus have both ({nu }_{0}+{nu }_{pi },(mathrm{mod},,2)) and ν_π, from which we can deduce also ν₀. Using this procedure for the time vortex in the resonantly driven Kitaev model, we see from Table 2 that ν₀ = 0, ν_π = 1, as expected.