Probing Majorana Modes via Local Spin Dynamics

We investigate Majorana modes in a quantum spin chain with bond-dependent exchange interactions by studying its dynamics. Specifically, we consider two-time correlations for the Kitaev-Heisenberg (KH) Hamiltonian close to the so-called Kitaev critical point. Here, the model coincides with a phase boundary of two uncoupled instances of Kitaev's model for p-wave superconductors, together supporting a degenerate ground state characterized by multiple Majorana modes. In this regime, the real-time dynamics of local spins reveal a set of strong zero modes, corresponding to a set of protruding frequencies in the two-time correlation function. We derive perturbative interactions that map the KH spin chain onto the topological regime of Kitaev's fermionic model, thus opening up a bulk gap whilst retaining almost degenerate modes in the mesoscopic regime, i.e., for finite system sizes. This showcases the emergence of Majorana modes in a chain of effective dimers. Here, the binding energy within each unit cell competes with the inter-dimer coupling to generate a finite size energy gap, in analogy with local energy terms in the transverse-field Ising model. These modes give rise to long coherence times of local spins located at the system edges. By breaking the local symmetry in each dimer, one can also observe a second class of Majorana modes in terms of a beating frequency in the two-time correlations function of the edge spin. Furthermore, we develop a scenario for realizing these model predictions in ion-trap quantum simulators with collective addressing of the ions.


INTRODUCTION
Topological modes are ubiquitous in many-body (MB) models, but their experimental detection and control in naturally occurring quantum systems can be challenging [1]. A prominent example is the Majorana fermion (MF), a non-Abelian anyon with non-trivial exchange statistics [2], which has been studied for a wide range of MB systems [3]. The perhaps simplest manifestation of an MF was proposed by Kitaev, who introduced a toy model for a fermionic quantum wire in the form of a one-dimensional (1D) p-wave paired superconductor [4]. The microscopic origin of this model was worked out for d 5 transition metals [5], and it is up until this day an important tool in the active pursuit of controllable MFs [6]. Kitaev's fermionic model is intimately connected to the Ising [7] and Kitaev-Heisenberg (KH) spin models. While the Ising model is ubiquitous and studied extensively in many contexts the KH model, with a potential realization in 4d 5 ruthenium trichloride α-RuCl3 [8,9], is less commonplace. The KH model attains frustration due to bond-dependent exchange couplings, and it may support long-range magnetic order [10,11] and quantum spin-liquid states (SQLs) [9,12]. Specifically in 1D the prospect for SQLs and topological modes has also been investigated [13,14].
For large-scale calculations of SQLs and topological modes, the use of quantum simulation with engineered lattice Hamiltonians in cold atom systems is a viable * Johannes.Bjerlin@matfys.lth.se pathway [15,16]. However, these systems typically require extremely low temperatures. A favorable alternative is given by simulators based on trapped ions [17][18][19]. Such setups are versatile and can function at comparatively high temperatures [20]. Recent successful examples of ion-trap simulations include a dynamical phase transition for a 53-qubits system [21], as well as quasiparticle dynamics in an Ising spin chain [22]. While quantum simulation promises remarkable speed up in the characterization of complex systems [23], so far most implementations have focused on well-studied static properties for which other highly effective numerical and analytical tools are available. Quantum simulation of dynamical features is hence especially compelling [24,25], as classical calculations are typically very costly [26,27]. While there have been some advancements in classical computation of time-dependent observables [28], there is to this date no general method to efficiently simulate the dynamics of large and strongly correlated systems.
An interesting theme, adjacent to quantum simulation and condensed matter physics, concerns the territory of few-to-many body physics. Here, recent advances in computational and experimental techniques (particularly within ultracold atomic gases [40,41]), has sparked experimental studies of, e.g., few-body magnetism without a lattice in one dimension [42], the formation of a Fermi sea [43] and, more recently, a few-body analogue of a quantum phase transition in two dimensions [44,45].
The advancements mentioned above highlight the growing interest in the controlled simulation of mesoscopic systems and number-conserving models with exotic features, which can shed light on the origins of quantum MB phases. Already, several studies have been conducted on MFs and topological phases in number conserving lattice models, motivated by the quest for a topological quantum computer [46]. This includes numerical studies of the topological features themselves, using density matrix renormalization group techniques [6,47,48], as well as studies focusing on the microscopic origins and possible realizations of the models in which they arise [21,47,[49][50][51]. This also extends to studies of dynamical observables [24,37,50], and in a recent preprint some of the few-body aspects of Majorana quasiparticles were laid out [52], underlining the promise of quantum simulation of few-body physics as a way to study complex MB phenomena using a bottom-up approach.
Here we focus on few-body phases of an interacting 1D quantum spin model (SM) that emulates Majorana edge modes (MEMs), investigating its dynamical features in the few-to-many body limit. The term "emulate" refers to the fact that the MEMs are topologically non-trivial only in the fermionic representation of the model [7]. Interestingly, the dynamical features of the SM still manifest a large discrepancy between bulk and edge. We begin by presenting an appropriate form of the Kitaev-Heisenberg Hamiltonian [6], using two parameters to tune the system between different phases around one of its critical points. We briefly discuss the various relevant phases in the static regime before investigating their individual dynamical signatures in local spin observables, focusing on MEMs. We use two-time correlation (TTC) functions to probe the Majorana bulk gap as well as the interaction-induced energy splitting between edge modes. This detection protocol elucidates the few-to-many-body development of MEMs without the requirement of deterministic preparation of any particular quantum state. Finally, we discuss a possible experimental realization of these findings in an ion-trap setup.
The static properties of the Hamiltonian are studied by means of exact diagonalization, using the full basis set ofŜ z i eigenstates. We utilize a sparse representation of the Hamiltonian and obtain the low-lying eigenvectors using the open-source library Eigen [53] developed for c++. For determining dynamical features, we numerically solve the time-dependent Schrödinger equation, using the fourth-order Runga-Kutta method for temporal discretization. Here, the sparse matrix-vector multiplication can be easily parallelized and distributed over multiple cores. Using this setup, we can currently treat systems of up to chain lengths L ∼ 20 on a single standard machine.

II. A TUNABLE MODEL FOR MAJORANA EDGE MODES
We focus on the 1D KH model describing an even number L of spin-1/2 subsystems interacting via nearest neighbor (NN) couplings. The unit cells consist of two spins, where the interaction inside the unit cell is different from the interaction between neighboring unit cells. The interaction between the spins is described by the Hamiltonian Here j is the unit cell index, and i is the spin index. S x i corresponds to a local operator of spin i, describing spin along the x-axis. We further set = 1, so that S D i = 1/2 · σ D i in terms of Pauli matrices σ D i . Initially neglecting the last term,V δ , this quantum MB Hamiltonian, with tunable parameters J and K, can be considered an inhomogeneous Heisenberg XY -model with exchange terms and additional sign-alternating double spin-flip interactions [6,54]. Similar models have been studied in the context of quantum phase transitions, criticality and magnetic long-range order [14,55,56]. We parametrize it in terms of a polar parameter, θ, governing the relative strength and signs of the interactions according to and we use √ K 2 + J 2 = 1 as the unit of energy throughout.
We begin by studying the phase diagram of the system withV δ = 0 close to the so-called Kitaev points, located at θ KP = ±π/2 → K = ±1, J = 0 and θ KP = 5π/4 ± π/2 → K = ±1/ √ 2, J = ∓1/ √ 2. Applying the Jordan-Wigner transformation, we can find the corresponding fermionic model (see Supplemental materials VII). The fermionic Hamiltonian can be directly decomposed into two separate systems, A and B, , located at (δ = 0, θ = π/2). At this point, the Hamiltonian describes a fermionic p-wave superconductor that supports a degenerate ground state with multiple Majorana modes. The hierarchy of degenerate multiplets at the KP is directly observed in the spectrum in subplot (a). Adding a termV δ with strength δ > 0 drives the system into the Majorana edge mode (MEM) phase, giving rise to a bulk gap, dividing all energy states into two sectors. Perturbing around the KP with the Y-bond interaction termV inter δ in Eq. 5 realizes a dimer Majorana edge mode phase, with each level attaining a 2 · 2 L degeneracy for L → ∞, whereas a non-vanishing gap is attained for finite L, as seen in subplot (b). For θ = π/2 − dθ and δ ∼ 0, shown in subplot (a), the system becomes a Tomonaga Luttinger liquid (TLL), which persists for moderate perturbation strengths |δ|. For θ = π/2 + dθ and δ ∼ 0, the system is in a spiral XY phase, also persisting for moderate perturbations. For θ ∼ π/2 and δ < 0, the system is in a gapped (G) phase with no MEMs. The MEM phases show a distinctly different behavior than the other phases in terms of the dynamical development of local spins. of length L/2, each corresponding to one instance of Kitaev's model for a p-wave paired superconductor [4] at the boundary point between the trivial and topological phase (see Supplemental materials VII). Exactly at the Kitaev points, only one of the subsystems A or B contributes energy in the Hamiltonian, so the full system acquires one free spin per unit cell, leading to groundstate degeneracies 2 L/2 and 2 L/2 -1 for open and closed chains, respectively [6]. For the open chain this amounts to L/2 Majorana operators, which are entirely absent from the Hamiltonian, so that the entire spectrum exhibits the same degeneracies as found in the groundstate. This global degeneracy is a stronger condition than what is usually required for general topological order [57]. The system here hosts multiple bulk Majorana modes distributed all across the chain, with a hierarchy of multiply degenerate states. Specifically, the highly degenerate groundstate multiplet is separated from the excited states by a gap, a necessary condition for the presence of non-Abelian quasiparticles [2,3]. Throughout the text, a globally N -fold degenerate spectrum means that each level in the spectrum is at least N -fold degenerate, but additional degeneracies may be present.
In this work, we focus specifically on the realization of MEMs around the Kitaev point θ KP = π/2. In Kitaev's original model, the MEM phase supports topologically protected modes at the edges [4], which correspond to a spontaneously broken spin-reflection symmetry when mapped to the Ising spin model [7]. We will nevertheless use the term MEM also in the spin picture.
To achieve the MEM phase in our setup we must invoke the additional termV δ into the Hamiltonian. Starting at the Kitaev point θ K = π/2 → K = 1 we map the system ontô with fermionic creation(annihilation) operatorsd † j (d j ).
Here K 1 = K 2 = 1, and j is the unit cell index (see Supplemental materials VII). Comparing to the Kitaev model [4], this gives the boundary point of the p-wave paired superconductor. Therefore, for topological modes the relative size of the first term must be decreased, so that |K 1 | < K 2 . We may thus either decrease |K 1 | or increase K 2 to enter the topological regime. We first consider the (local) energy term proportional to K 1 and map this back to the spin picture (see Supplemental materials VIII), revealing the appropriate perturbation term,V with the MEM phase occuring for |K 1 + δ| < K 2 . This term corresponds to interactions within a unit cell of two spins. We can compare this situation to the equivalence of the transverse-field Ising model and the Kitaev model [7,58], where the local fermionic energy term maps onto the local energy of a single spin in a magnetic field. For our case, each term in Eq. 4 instead represents the local energy of the unit cell dimer j. Precisely at the Kitaev point, where K 1 = K 2 = 1, the dimer energy equals that of the inter-dimer bond, and the system remains gapped for L → ∞. Here, the global degeneracy is that of L/2 dimers with one free spin each, giving 2 L/2 states. For |K 1 | < K 2 the inter-dimer bonds instead dominate, and an additional global two-fold symmetry arises for L → ∞, corresponding to zero-energy Majorana modes, giving a global degeneracy of 2 · 2 L/2 . The degeneracy is perfect in the limit of infinite chains, whereas the finite size gap between the two degenerate multiplets scales with e -L/2 . Aside from additional degeneracies, the energy spectrum of this system coincides perfectly with that of a transverse field Ising model with L/2 spins,Ĥ I = K 2Ŝ y iŜ y i+1 + K 1Ŝ x i . The Hamiltonian 1 is thus very similar to the transverse field Ising model but differs in its dynamical properties due to the additional degeneracies.
As noted above, we can also enter the MEM phase by increasing the relative size of the terms scaling with K 2 in the Hamiltonian 3, givinĝ with the unit cell index j and the MEM phase occurring for |K 1 | < K 2 + δ. This term corresponds to interactions between two unit cells. We now invoke a third option for the perturbing interaction,V δ , corresponding to a fully connected Ising term, where i is the spin site index. This perturbation does not map the fermionic Hamiltonian onto a Kitaev model, but we nevertheless see the emergence of an MEM phase for |K 1 | < K 2 + δ. We will see that this perturbation simultaneously creates MEMs and breaks local symmetries within each dimer, giving rise to a beating pattern in the time-dependent edge spin correlation functions. In conclusion, we use δ as an effective parameter that controls the onset of MEMs, using either of the perturbations in Eq. 5 or Eq. 6. The two different perturbations are used to highlight two different effects in dynamical simulations of the MEM regime. The spectrum due to the inter-dimer perturbation Eq. 5 is depicted in Fig. 1.

III. PHASE DIAGRAM AND STATIC PROPERTIES
Let us now briefly discuss the four phases in the phase diagram shown in Fig. 1, spanned by the parameters θ and δ in the vicinity of the critical Kitaev point at θ = π/2 and δ = 0. The characterization of these phases will be helpful when discussing the dynamical features of local spins in the later sections.
• Majorana Edge Mode (MEM) phase (θ = θ KP = π/2, δ > 0): the bulk energy spectrum in this regime is gapped, with two zero-energy edge modes in the thermodynamic limit. However, in finite systems, their energies remain small but finite, yielding a finite-size gap that vanishes exponentially with increasing system size. For the perturbing term,V inter δ , each level has a global 2 L/2 -fold degeneracy due to free parameters in the Hamiltonian, so that for L → ∞ the spectrum becomes 2 · 2 L/2 -fold degenerate. We call this the dimer MEM phase. ForV Ising δ the spectrum becomes globally two-fold degenerate for L → ∞. We call this the Ising MEM phase.
Further characterization and discussion of the static properties of these phases are presented in Supplemental materials IX.

IV. DYNAMICAL FOOTPRINTS OF THE MAJORANA MODES
Time-dependent observables are a powerful tool for the analysis of physical systems beyond their groundstate phases [25,30,31,34,36]. In particular, local measurements,Ŝ D i (t), of a spin i along D ∈ {x, y, z} are intuitive and experimentally accessible probes that can be used to highlight the emergence of MEMs [36,59].
We first consider the spin operator, The eigenstates ofĜ D are denoted | ± Φ D n = |s D 1 s D 2 s D 3 ... , where ± denotes positive or negative parity, respectively. These eigenstates will serve as the initial states for the dynamical simulations, where we numerically evolve each state in time under the Hamiltonian operator 1 and study the dynamical evolution of local spins. We also note that this spin-operator flips all spins along the axes perpendicular to D, i.e., along x and y for G z .
For quantitative measures, we consider the mean autocorrelation function, where the sum over N produces the average over a randomly sampled set of N initial states | ± Φ D n ∈ {| ↓↓↓ ... , ..., | ↓↓↑ ... } in the basis of spin D. We also consider the (discrete) Fourier transformed evolution func- , again taking the average over a large set of initial states, We further calculate variances to highlight which features are largely independent of the particular input states we choose. 1 By sampling over multiple initial states and taking the average, we specifically access robust features of the system in the sense that an experimental setup would not rely on repeated and deterministic preparation of any specific initial state. Measurements can instead be performed with mixed states for those spins which are not directly probed, which is especially relevant for detection of strong zero modes [36,37]. To simplify the computations we, however, consider pure initial states for the individual runs and take the average afterward, i.e. we essentially perform a Monte-Carlo sampling of a completely mixed density matrix. Fig. 2 shows mean autocorrelations for a set of randomly sampled states developing in time under two different Hamiltonians. As will be discussed in the following section, there are several robust features in the mean autocorrelations of edge spins (like the constant spin-y projection in each plot), even though they essentially represent time-development of mixed states.

A. Zero modes of the Hamiltonian
A Hamiltonian that supports MEMs can be represented in terms of Majorana operators in such a way that some of them drop out of the Hamiltonian in the infinite (L → ∞) system limit [4], giving rise to global symmetries and corresponding degeneracies in the entire energy spectrum. In line with the procedures in Refs. [4,35,36], we elucidate the dynamical properties of the finite-size system by first deriving the zero modes, which correspond to the Majorana modes in the corresponding fermionic model. These modes are constructed to approximately commute with the Hamiltonian, with corrections ∼ e −L , and are associated with long-time coherent features in the dynamical evolution.
We first consider the Hamiltonian 1 at the Kitaev point θ = π/2, using the perturbationV inter δ , so that For convenience, we have rescaled the Hamiltonian in the last line, so we end up witĥ where we recognize that δ > 0 yields the dimer MEM phase.
We initially aim to construct a zero mode corresponding to a local spin along y at the edge of the spin chain. This means that we simply put for the zeroth-order approximation of the mode operator Ψ which we offset by introducing first and second order terms Ψ (1) where the factor of four is absorbed by an emerging op-eratorŜ y 3Ŝ y 3 = 1/4. The first-order term actually commutes with the full Hamiltonian, so we can freely choose the constant M 1 . Ψ (2) A does not commute withŴ , which again can be offset by additional terms of higher order. Continuing in this fashion, we get which we can show commutes with the Hamiltonian up to an exponentially small factor. In this expression, we have changed the representation to Pauli spin operators, S D i → σ D i /2, and separated the odd and even orders for convenience. Each term in the sums now corresponds to a Majorana fermion [4]. For the second sum we have one free choice for the constant M i per term, meaning each unit cell adds a degree of freedom for the zero mode. This is consistent with the global 2 L/2 -fold degeneracy in the spectrum. For simplicity, we have chosen M j = M j-1 .
Each term in the sum of Eq. 11 anticommutes with all the others, so that and we proceed to choose the normalization so that Ψ 2 A = 1 for |M| < 1 and δ > 0. It is now clear the each of the zero modes constructed in Eq. 11 commutes with the Hamiltonian H, now in matrix representation, up to an exponentially small term, where G z is the spin-flip operator in Eq. 7. We also note that the zero mode anticommutes with the spin-flip operator so that {G z , Ψ A } = 0. This means that G z toggles between different eigenstates of Ψ A and vice versa. Together with normalizability and the vanishing commutator ε rem , these properties of the zero mode constitute the necessary conditions for long-time edge spin coherence [36]. Alternatively, zero modes can be constructed similarly by starting at the other edge,Ŝ y L . We also find that forŜ y 2 andŜ z 1 , which also commute withĤ 0 , two corresponding zero modes, Ψ B and Ψ C can be found (see Supplemental materials X). These modes differ in their construction, since they are each derived from different starting points. However, the different zero modes Ψ A , Ψ B , and Ψ C each span the same operator space, so that the implicated global degeneracy of the energy spectrum for L → ∞ is still 2 · 2 L/2 . In summary, we can construct normalizable zero modes Ψ A ,Ψ B , and Ψ C for the Hamiltonian 10, which includes the perturba-tionV inter δ with δ > 0. Interestingly, the zero modes Ψ A and Ψ B also commute, up to an exponentially small factor, with the HamiltonianĤ which instead uses the perturbing Ising termV Ising δ . This is true for δ > 0, providing we set all constants in front of odd orders to zero, so that M A = 0 and M B = 0. This means that for each of the operators,Ŝ y 1 ,Ŝ y 2 , one may use the same construction of zero modes for Hamiltonian H Ising δ as forĤ inter δ . ForĤ Ising δ the free parameters in the zero modes are however removed, so the global degeneracy in the spectrum becomes only two-fold in the limit L → ∞. The commutations between the finite-size zero modes and the HamiltoniansĤ Ising δ andĤ inter δ are identical, so the finite-size gap between zero modes ∆ L are also the same for the different perturbations. We will see that the long-time coherence depends crucially on this gap, and we therefore expect some identical long-time features for both Hamiltonians.

B. Long time dynamics of edge spins
We now proceed to summarize the impact of zero modes on the long-time dynamics of our system. For details, we refer the reader to the Supplemental materials XI, and for a comprehensive theoretical background to Ref. [35,36]. We evaluate the autocorrelation function Γ D 1 (t) for an eigenstate |S D , with corresponding eigenvalue s D 1 , of the edge spin operatorŜ D 1 along direction D. In practice we will use states | ± Φ D n which are eigenstates of G D , but here we consider a general state |S D . We get for the autocorrelation where n|, m| are eigenstates of some HamiltonianĤ with corresponding zero mode Ψ A . Since Ψ A (almost) commutes with the Hamiltonian we may divide all energy states into two sectors denoted by positive or negative sign, corresponding positive or negative eigenvalues of Ψ A so thatĤ We can now re-write the autocorrelation function with new indicies For long times t and large system size L, terms with n A = m A add up incoherently while terms with n A = m A add up coherently, i.e. terms with n A = m A get a random phase so that we can ignore them. The double sum may then be approximated for long times by with terms T D 1 -T D 4 relating to the time-independent matrix elements between |n A ± .
This shows that for an infinite system, the long-time spin oscillations are stable. For finite systems this is no longer the case, and the oscillations will eventually decay. The coherence time, i.e. the time during which the spin autocorrelation function remains stable, either displaying a finite value or a persistent oscillation, is generally set by the commutation between the Hamiltonian and the zero mode, which vanishes with L → ∞ [36]. Interestingly, if the finite size gaps ∆ L between semi-degenerate states in a systems' spectrum are all identical, spin autocorrelations which first appear to decay will have a revival time of order 1/∆ L ∝ 1/ ε rem . If one other hand the gaps are different we will only have partial revivals. For integrable systems, like the Ising model, these revival times may be directly calculated [60]. This is also true for the Hamiltonian in Eq. 10, which will be apparent from the dynamical simulations.
1. Coherence for σ y 1 ,σ y 2 and σ z 1 So far we have not specified the direction D in which we aim to measure the spin, and we proceed to study the effect of two particular choices ofŜ D 1 based on the timeindependent terms in Eq. 17. Since we specifically use the Pauli-spin representation of operators in this chapter, we will use σ D i to represent a spin operator at site i along D. We begin with the HamiltonianĤ inter δ and σ y 1 with corresponding zero mode Ψ A . For T y 2 and T y 3 we exploit the spin-flip operator in Eq. 7 which anticommutes with Ψ A so that G z |n A ± = |n A ∓ and we get Now we note the anticommutation {σ y 1 , G z } = 0, leading to T y 2 = T y 3 = 0. For T y 1 and T y 4 we instead employ Eq. 16 For T y 4 we may use Eq. 11 directly, giving {Ψ A , σ y 1 } = N e + C, where C represents (exponentially) small corrections. For This gives the long-time limit of the autocorrelation which only depends on the initial state and how much overlap it has with each sector of eigenstates for Ψ A . The exact form of the corrections C depend specifically on the model [36], but they are always exponentially decreasing with L/2. We can use an identical derivation for σ y 2 by making the substitutions Fig. 2 a shows the simulated dynamical development of spins forĤ inter δ , confirming that the mean autocorrelation for σ y 1 ,σ y 2 and σ z 1 is long-lived compared to other spins. This is directly explained by the fact that we can construct corresponding zero modes, as shown in the previous section. This is not true for the other operators shown in the plot, where the autocorrelation vanishes for long times.

C. Beating patterns for edge spins
In Fig. 2 we study how the autocorrelation function compares for the HamiltoniansĤ inter δ andĤ Ising δ . The results for σ y 1 and σ y 2 are the same for the different Hamiltonians, whereas σ z 1 is strikingly different. Curiously, long time coherence is still present forĤ Ising δ , but with an oscillating factor which we find is independent of system size. The coherence time of the oscillation is however set by system size, as for σ y 1 and σ y 2 . We can relate this result directly to the zero modes. We derive in the Supplemental materials XII that We see that the expression by symmetry is, except for the time-dependent factor, identical to the autocorrelation Γ y 1 (t), but here for an initial state |S z . This precession of the edge spin σ z 1 is hence given by an oscillation, with frequency δ, and an envelope function given by the coherence time of σ y 1 . Comparing the evolution of different spin-components in Fig. 3 a we see that σ x 1 has a spin precession which does not decohere, in contrast to σ z 1 . We note that the x-component of the edge spin commutes withĤ KH (θ = π/2), and using the same technique as for σ z 1 one may again calculate the spin precession of σ x 1 from the pertur-bationV Ising δ . This gives a time dependence ∝ (cos 2 δt − sin 2 δt) without any decoherence, owing to the fact that σ x 1 does not couple the different zero mode eigenstates. We find that the spin precession frequency of σ x 1 and σ z 1 is independent of the chosen initial state, since it explicitly depends on global gaps ∆ δ in the spectra. These correspond to the gap within the edge unit cell, given by the perturbation δŜ y i S y i+1 . Since σ x 1 toggles between levels within the unit cell split by δ, but not between zero mode eigenstates, it only corresponds to the oscillation frequency from the gap ∆ δ . The operator σ z 1 on the other hand, toggles between both zero modes and unit cell energy levels. The decoherence time now relates directly to the toggling between different zero modes, corresponding to gaps ∆ L given by the commutation of the zero modes and the Hamiltonian. This correspondence is shown in Fig. 4, where the autocorrelation is evaluated for longer times. Away from the Kitaev point θ KP the global gaps are no longer present, so the oscillations decay quickly in both the TLL and XY phase. are no global degeneracies in the Hamiltonian, and oscillations at the edges simply decay in the same manner as for spins in the bulk, as seen from Fig. 5. Simulations show that the decoherence becomes more profound, meaning that the system is fully decohered for a longer time without signs of revival, as the systems grow in size. We stress that the edges are still interacting with the bulk in this MEM phase, which becomes apparent from the fact that the coherence times increase with L (not shown here). This means that we cannot think of the edges and bulk as two completely separate systems described in terms of a tensor product between them. However, observables like Γ x 1 (t) may still have trivial behavior if they fully commute with the Hamiltonian, as seen in Fig. 3.
From the above results, we conclude that the autocorrelation function can reveal a clear signature for the existence of zero modes in the Hamiltonian, providing the correct spins are measured.
The oscillations are only visible for the edges and not in the bulk since the bulk spins generally decohere fast [35]. A handy explanation is offered in the limit δ 1 where we notice that applying a local spin-z at the edge adds the energy cost of breaking (or creating) exactly one antiferromagnetic bond: The operatorσ z 1 , which acts as a spin-flip operator on a local spin-y state at the edge, will necessarily break or create exactly one bond, connecting states separated in energy by the gap ∆ δ = δ. These gaps are present throughout the spectrum and result in a large set of coherent terms in the autocorrelation function 15, giving a significant contribution for long times. Forσ z i , acting on a bulk state L > i > 1, a spin flip is instead associated with either simultaneous creation and destruction of one bond (or simultaneous creation or destruction of two bonds) and will thus con- The characteristic frequency ∝ ∆ δ does not change with L. The revivals of the oscillations correspond to a beat frequency f beat ∝ ∆L which is, however, directly dependent on the zero mode gap, ∆L, which in turn depends on L. The expected beat frequencies calculated directly from the spectral gaps are T L=8 beat = 259 and T L=10 beat = 673, which is in good agreement with the observed dynamics.
nect states with energy differences smaller than the gap. This gives rise to a set of low frequencies contributing to the temporal evolution of the bulk spin, causing an effective decay of the oscillations. In the limit of small δ, this discrepancy between bulk and edge is reduced since ∆ δ is then of the same order as other, small, gaps in the spectrum. From the dynamical simulations, we find that the number of dominating frequencies for the first bulk spin corresponds to N peaks = L − 2 which means that for successively larger systems, the oscillations in real time will be washed out for this site. The high relative variance for the bulk spins, see Fig. 3 furthermore shows that the exact dynamics here depend heavily on the input state. Directly at the Kitaev point, the number of dominating frequencies corresponds to the system size L (N peaks = L/2). The frequencies are independent of the site index and input state, whereas their relative weights depend on the site index. This explains why no coherent oscillations are seen directly at the Kitaev point.
In the limit of long-time spin precession for the MEM phase we observe decay and revival of the oscillation at the edge, resulting in a beating pattern. The beating pattern, with period T L beat , observed in Fig. 4 is a finite-size effect directly related to the revival of decohering spins, which was noted for the pure Ising model in Refs. [36,60], where the revival time scaled with system size. The envelope function in Eq. 22 is given by the same function as for the spin σ y 1 , which has some important consequences. We know that the decoherence and revival times of the autocorrelation function here are related to the correction term C [36], so these properties are therefore identical for σ y 1 and σ z 1 . The revival time scales with the few-body gap ∼ 1/∆ L , as evident in Fig. 4. As noted before, ∆ L is identical forĤ Ising δ andĤ inter δ . We have shown that the long-time properties of these two systems can be mapped onto those of an effective Ising model with two-site unit cells. The complete revival and resulting beating pattern in the autocorrelation plots of Fig. 4, scaling with T L beat ∝ 1/∆ L , is therefore not surprising.

V. QUANTUM SIMULATION WITH TRAPPED IONS
For a possible quantum simulation of the MEM phase, we here sketch a setup with ions trapped utilizing radio frequency (RF) fields [61][62][63]. Experiments with such systems typically realize effective Ising, XY, or XYZ spin-spin interactions, and may be supplemented with global transverse fields terms [64,65]. These setups have been used in a large number of studies, for example simulation of quantum magnets [66] entanglement propagation [22], and variations on the quantum Ising spin chains [21,67,68], along with more general quantum computing implementations [69,70]. Whereas most quantum simulation experiments with trapped ions use a linear configuration we will consider a geometry where the trapping frequencies perpendicular to the tap axis are very different. For suitable parameters this ensures that the ions form a planar zig-zag structure as shown in Fig. 6. The use of a zig-zag configuration allows for the direction and magnitude of a laser field to control the size and sign of effective interactions between different rungs of the zig-zag spin ladder [72] and we will exploit this below. Furthermore, we assume that every third ion in our setup is selectively hidden, requiring individual addressing of ions [73], so that they do not participate in the simulation. The resulting pattern is sketched in Fig. 6. To ensure a uniform distance between ions, we only consider the central part of a crystal containing N ∼ 100 ions and assume that all other ions do not participate in the simulation, e.g. because they have been optically pumped to a different internal level [73], see Fig. 6. Alternatively, a more uniform distribution could also be obtained by carefully adjusting the local trapping potential [71].
For convenience, we first employ a simple rotation on the Hamiltonian in Eq. 1 so that where the subscript "sim" denotes the axes of the simulated Hamiltonian. We aim to realize the Hamiltonian with the purpose of simulating the effective Hamilto-nianĤ Ising δ =Ĥ KH (θ = π/2) +V Ising δ from the previous section. We note that the indices j correspond to the effective indicies in the simulated Hamiltonian, which correspond to the active (not hidden) ions inside the ion trap, as shown in Fig. 6  We assume that all ions outside the central region are hidden. so that the distance between different rungs is roughly uniform for the participating ions. Alternatively, a similar uniform spacing can be achieved by shaping the trapping potential [71].
As we discuss below a major challenge is to remove additional interactions induces by the coupling mechanism corresponding to next-nearest neighbor (or higher) interactions in each of the effective ZZ and XX interactions. The residual interactions will cause decoherence and need to be reduced enough for the beating mechanism to be observed at different system sizes L. We find that a suitable parameter regime to aim for is δ eff ∼ 1.
The signature of the MEM phase ofĤ Ising δ requires the edge spin to be initially prepared in an eigenstate ofŜ y 1 . This can be achieved by optical pumping and subsequent ±π/2 rotations around the x-axis [74,75]. The remaining spins can be prepared into any mixed state, but to be consistent with the previous sections we here assume their initial states to also be eigenstates ofŜ y i . With the initial spin-state s y 1 known, the autocorrelation function Γ y i (t) in Eq. 8 can be evaluated by measuring the spin along y at a later time t.
We now proceed to discuss the realization of the two different interactions inĤ sim , starting with the ZZ interaction. For a detailed derivation and discussion of parameters of the ion-trap simulation, see the Supplemental materials XIII.

A. Simulating ZZ and XX interactions
For the ZZ interaction we consider a two-photon Λ scheme where the lasers coupling two stable ground levels | ↑ | ↓ to an excited state are far-off resonant with the dipole allowed transition and the detuning is given by D which is much larger than the spontaneous decay rate of the system. In line with Refs. [64,72] we employ a pair of laser fields with effective Rabi frequencies Ω 1,ς and Ω 2,ς for ς =↑, ↓ coupling the ground levels ς to the excited state. The effective Raman wavevector of the two fields is given by k a l = k 1 − k 2 = k a l {cos ϕ a , sin ϕ a , 0}, and can be tuned via alignment of the lasers. The laser beatnote ω a l = ω 1 − ω 2 of the fields is chosen close to the ions' collective vibrational motion in the transversal direction y, with mode energies ω y p , and far-off resonance with the vibrational modes in the zig-zag plane with mode energies ω xz p . The transverse vibrational modes act as mediators of an effective spin-spin interaction of the canonically transformed Hamiltonian [64], and by carefully choosing the detunings and alignments of the laser fields a nonisotropic effective ZZ interaction with tunable strength and range can be realized [72].
The overall strength of the interaction is controlled by the magnitude and direction of the laser fields, affording some freedom in choosing the parameters in our effective Hamiltonian. Crucially, the factor V a ij is dependent the alignment of the field and the relative equilibrium positions r 0 ij of the interacting ions, so that the nonhomogeneity of the interaction can be tailored via the fields.
The effective two-photon Rabi frequency is given by Ω a = Ω 1,↓ Ω * 2,↓ + Ω 1,↑ Ω * 2,↑ /2D. For the zig-zag configuration of the ions, it is convenient to choose k a l in the xy-plane, i.e. perpendicular to the direction of the rungs, as shown in Fig. 6 so that Same rung: cos k a l ·˜ r 0 ij = 1 Different rungs: By choosing the angle ϕ a we can now eliminate the interaction between different rungs, even when the optical wavelength is small relative to the mutual ion distances [72]. We numerically calculate the equilibrium configuration of the ion trap with N = 70 to find the ions' positions and their transverse vibrational eigenmodes. We employ these quantities to evaluate the full expression for the effective interaction given in the Supplemental materials XIII. Choosing the detunings of the fields, relative to the vibrational mode energies along y such that ω y p − ω a l /ω y p ∼ 0.05, we find that the typical distance dependence of the interaction strength becomes ∼ 1/|i − j| R with R ∼ 3. Since this falls off of quickly with the distance, the interaction will be dominated by NN interactions [76]. We see from Fig. 6 a that by hiding every third ion, we can map the same rung NN interaction to odd indices in the effective system whereas different rungs corresponds to NN interactions starting on even indices, i.e.
for the effective indices i of the active ions. For the effective XX interaction we need to drive a transition between two internal levels. This can either be done directly or as a two-photon Raman transition.
The effective spin coupling is implemented via the vibrational sidebands of the transition [20]. The angular frequencies of the driving are given by ω b Here ω b is the transiton frequency between the considered internal levels of the ions. The detuning ω b l is roughly matched to the transverse trapping frequency ω y so that ω b l ≈ ω y . We aim to virtually excite the vibrational sidebands, and we employ the sideband detuning γ b = ω b l − ω y . The sideband detuning is chosen to be positive, so that the ω b l lies above the highest vibrational mode along y (out of plane), which is the centre of mass mode ω c.o.m. . We employ the same vibrational branch as for ZZ, but we assume there are no interference effects between the processes implementing ZZ and XX. This can be achieved by ensuring that the frequencies are incommensurate. The virtual phonon exchange between ions induces an effective interactioñ where the vibrational mode eigenvectors M p and mode energies ω p are for the strongly confined direction y [77]. We again employ the calculated vibrational modes and ion positions to explicitly calculate the effective interaction. The overall strength of the XX interaction can be controlled by the effective Rabi frequency Ω b [20]. Choosing the detuning ω y p − ω b l /ω y p ∼ 0.05, we may realize approximate interactions ∼ 1/|i − j| R with R ∼ 3 for the XX interaction, so that residual terms are of similar order as for ZZ. The detuning is chosen to be different than that for the ZZ interaction, but on the same order of magnitude. Since there is no angular factor in Eq. 28 the effective interaction connects both sites within the same rung and sites between rungs.
Putting everything together, we obtain for the zig-zag indicies i, ĵ We now use the mapping in Eq. 27 to transform into the indices i j for the active ionŝ This form of the Hamiltonian agrees with the desired model in Eq. 24 apart from the residuals in the last line. According to the arguments above these residuals can be rather small. Furthermore, for a translationally invariant system, the desired coefficients will be identical between units cells as desired and can be adjusted to the desired values ratio to realize the MEM phase of the Kitaev-Heisenberg model.
For a real ion trap, the parameters can suffer from numerous imperfections. In particular, for a standard ion trap, the density will be higher near the center of the trap and the ions will not be equidistant, as shown in Fig. 6. In principle this can be overcome by carefully designing the trapping potential [78,79], but below we explore the limitations imposed by operating in a standard ion trap with harmonic confinement.
To investigate the imperfections in a real ion trap we consider the situation depicted in Fig. 6, consisting of 70 trapped ions, see the Supplemental materials XIII for further details. To evaluate the role of imperfections the coefficients V a ij ,V b ij and residuals are evaluated for the numerically calculated ion positions and eigenstates. We specifically do this by first fixing |k a l |,|k b l | and then the angle ϕ a so that different rung interactions disappear for ZZ, and proceed to choose the exact detunings for all fields. We can then finally match the Rabi frequencies Ω a and Ω a so that V a 2i',2i'+1 ∼ V b 2i',2i'+1 ∼ V b 2i',2i'+1 /2, mapping the simulated Hamiltonian 30 ontoĤ Ising δ = H KH (θ = π/2) +V Ising δ with δ eff ∼ 1, with additional residual terms. Using the numerically coupling coefficients we can then proceed to evaluate the dynamics in the MEM phase corresponding to the time evolution in Fig. 4.
The results of a numerical simulation, performed for ∼ 40 sampled initial states, of the ion trap are shown in Fig. 7. For the simulated setup the residual terms, the largest of which correspond to ∼ 10% of V a and V b , clearly have a large effect on the results. The beating pattern visible for the idealized Hamiltonian in Eq. 24 is hard to observe due to the rapid decay of oscillations in the autocorrelations. We can however see the sizedependent revivals for the long-time coherent spin along one of the axes, here along x (Note that when comparing to the idealized Hamiltonian we use the rotation 23). In Fig. 4, showing the dynamics of the idealized Hamiltonian, we further had two precessing spin components which were enveloped by a long-time beating.
Here we also saw that the decoherence of one of the precessing spin components (z) was enveloped by the long-time coherent autocorrelation function for the ycompoenent, while the other precessing spin component (x), evolves more independently from the long-time coherent spin. This is behavior is somewhat visible also in the ion-trap simulation when comparing Fig. 7 a and  Fig. 7 b. We see that the oscillation of |Γ z 1 (t)|, which does not correspond to a zero mode, is less dependent on system size than the other spins. We conclude this by noticing that the oscillation of |Γ z 1 (t)| is not suppressed in the region around the first node in |Γ x 1 (t)| for L = 8. However, the amplitude of the oscillating |Γ y 1 (t)| is suppressed in this region. For L = 10, the oscillation in |Γ y 1 (t)| is instead suppressed at slightly later times. This is due to the longer coherence time of |Γ x 1 (t)|, which envelopes the oscillation, as shown in section XII. It is however evident that the precession in |Γ y 1 (t)| also suffers from the same type of decoherence as |Γ z 1 (t)|, since there is no visible revival for |Γ y 1 (t)|. This gives an estimate for the influence of the residual interactions in the simulated Hamiltoninan, since |Γ z 1 (t)| should not decohere in the limit of zero residuals. (Compare to the evolution of the corresponding |Γ x 1 (t)| of the idealized Hamiltonian in Fig. 3 a ) The behavior of the spins away from the edge (not shown) is similar for the simulated Hamiltonian as for the idealized Hamiltonian, where the long time coherence is seen also for s 2 but not for s 3 . There is also no spin precession away from the edge for either spin.
In the ideal case, the nodes in the beating and longtime coherent spin autocorrelation can be used as a direct measure of the finite-size Majorana gap ∆ L = E 1 − E GS , which effectively simulates finite-size scaling in the MEM phase. The residual interactions in the quantum simulation however cause the oscillation of the edge spin to decohere rapidly, generally before the first node. Furthermore, the residual interactions destroy the global degeneracy corresponding to the zero modes, so the gaps throughout the spectrum are no longer homogeneous for the simulated Hamiltonian. By making a more homogeneous distance between the ions and implementing additional fields at different detunings and angles into the quantum simulation, the residuals can be reduced. This would allow for the decoherence time t z c to be increased so that the first node, and subsequent revival of the oscillation in |Γ y 1 (t)| can be observed for successively larger systems.

VI. CONCLUSIONS AND OUTLOOK
We have showcased dynamical features of the Kitaev-Heisenberg Hamiltonian, particularly focusing on the behavior around the so-called Kitaev point at which a multi-degenerate set of Majorana modes appear. By perturbing the model with a nearest-neighbour termŜ yŜy on even sites, we see that the model can be mapped directly onto the Majorana edge-mode (MEM) phase of Kitaev's model for p-wave paired superconductors, with one additional degree of freedom per unit site. By perturbing with a uniform Ising termŜ yŜy between all sites, the additional degrees of freedom are lost but the system retains the Majorana edge modes.
Studying the direct time-development of local spins under this Hamiltonian we show that the MEM phase of the latter system can be identified via the precession of edgesite spins, which oscillate with two dominant frequencies.
This gives rise to a beating pattern corresponding to the indicates that by increasing the number of active ions in the simulation, one can observe a heavily damped beating pattern in |Γ y 1 (t)|, corresponding to the undamped beating pattern observed for |Γ z 1 (t)| in Fig. 3 a (for the idealized Hamiltonian). Note that the simulated Hamiltonian is rotated relative to the idealized Hamiltonians via Eq. 23, so that the oscillation for |Γ y 1 (t)| indeed corresponds to the oscillation for |Γ z 1 (t)| in Fig. 3 a. We see that |Γ x 1 (t)| exhibits revivals of the long-time coherence scaling with L. This is due to the presence of a strong zero mode in the Hamiltonian. The coherence time t z c of the oscillation in |Γ z 1 (t)|, on the other hand, is more independent of system size since it does not correspond to a zero mode. Its decoherence is explained by the residual terms in the quantum simulation and gives an overall estimate of the influence of impurities in the simulation. The oscillation of |Γ y 1 (t)| should exhibit a beating pattern in the limit of an idealized Hamiltonian, but decoheres on the same time-scale as |Γ z 1 (t)|. There is however an indication of the first node in the beating in |Γ y 1 (t)| for L = 8, which occurs before the system has fully decohered. This is seen from the node in |Γ y 1 (t)| at t ≈ 6, which is pushed towards longer times for larger system sizes L.
finite-size energy gap between the semidegenerate Majorana edge-mode states. The precession frequency of the spin is set by the interaction within the outermost unit cell of the system, while the beating is enveloped by the long-time coherent spin dynamics, which depends on the system size L. This is analogous to the long-time coherence of edge-spins in the Ising model, an effect caused by strong zero modes present in the system [35,36]. We show that zero modes affect the long-time coherence of the spins measured along two different axes in our model, while the dynamical evolution of a spin prepared along the third axis is independent of zero modes, and therefore independent of system size L. These characteristics are not present in other parts of the phase diagram, or for spins in the bulk. Crucially, this method of studying the dynamical properties does not rely on the repeated and deterministic preparation of a single initial state. It instead only re-quires deterministic preparation of a single spin, while the remaining spins can be randomly distributed.
We sketch an ion-trap quantum simulation, in which the steady-state zig-zag configuration of harmonically confined ions is exploited to realize the MEM phase. We see that some finite-size scaling properties of the spin dynamics can be observed, even for a setup with relatively large residual interaction terms. If the residuals were to be further reduced, our setup and readout mechanism could realize a quantum simulation of an interesting numerical challenge; the finite-size scaling of collective Majorana edge modes.

VII. MAJORANA MODES AND SPECTRAL PROPERTIES AT THE KITAEV POINT
Here, we discuss the decomposition and transformation of the Hamiltonian 1 referenced in section II of the main text. These transformations are performed to highlight the appearance of Majorana modes in the model. For an intuitive understanding of the spectrum of the Hamiltonian 1 we may rewrite it aŝ We can fermionize this Hamiltonian by applying the Jordan-Wigner transformation and writinĝ wheref † andf are the fermionic creation and annihilation operators respectively. In the spirit of Ref. [6] we rewrite the sum in terms of the unit cell index j, where each cell contains a black(b) site and a white(w) site. This transforms the Hamiltonian intoĤ with two termsĤ Like in Ref. [6] we define Majorana operators aŝ We can now writeĤ Finally, we define two independent non-local fermion operatorŝ This allows us to finally writeĤ which describes two independent p-wave superconductors at the boundary point of the Majorana edge-mode phase. At the Kitaev point (θ = π/2), only one of the Kitaev chains contributes energy in the Hamiltonian, and the system gets one free spin per unit cell, leading to degeneracies 2 N/2 and 2 N/2−1 for a non-periodic and periodic chain, respectively [6]. Slightly tilting θ away from the Kitaev point this degeneracy is removed, but we still retain the multiplet structure.

VIII. RETRIEVING THE GAPPED EDGE MODE PHASE IN THE SPIN-REPRESENTATION
Here, we construct the appropriate form of perturbations to the Kitaev-Heisenberg model to study Majorana edge modes, as discussed in section II of the main text. For definitions and background consult sections II and VII. Starting with the fermionic Hamiltonians 41 and 42 derived in the previous section we note that going away from the critical point into the edge-mode phase requires a decrease in the relative size of the local term(s)d † id i andd † id i , respectively. We therefore perturb each of the two Kitaev chains. After mapping back to the spin-representation we find that for θ = π/2 the appropriate operator isĤ where index j denotes the unit cell index. We map this change back to the fermion operators.
Neglecting constant terms and applying the Jordan-Wigner transform we are left witĥ and similarlyĤ At the point θ = π/2, the Hamiltonian is equivalent toĤ A . Note that the appropriate sign of the perturbation δ depends on which Kitaev point we consider. As shown in Ref. [6], the absent Majorana operators map back to the spin operatorsb each of which has one free index per unit cell, giving a degeneracy of the ground state.

IX. PHASES AROUND THE KITAEV POINT
We here further characterize the phases around the critical point δ = 0, θ = π/2 forĤ KH + V Ising δ . In addition to the spectral properties discussed in the main text, we here calculate the static spin structure factor where the unit cell length is taken as 1. We also consider the von Neumann entanglement entropy of a subsystem with length lS with the reduced density matrix Tr L−l ρ for the full density matrix ρ. The relevant phases are sketched in the phase diagram in Fig. 1, and they are found to largely agree to those found in Ref. [6].

X. ZERO MODES
In this section we perform a derivation for the zero modes, i.e. low-energy modes which approximately commute with our Hamiltonian (see section IV A in the main text). The presence of such modes largely explains the intriguing features of the autocorrelation simulations for the edge spins. Derivations follow the procedures in Refs. [4,35,36], where higher-order corrections to a first-order approximation of the mode are sequentially added by commuting the mode with Hamiltonian.
A. Zero modes forŜ y We start with by considering the commuation between the operatorŜ y 2 , our first-order approximation for the zero mode, and the Hamiltonian. This gives . This entanglement measure is calculated for a range of system sizes L in each of the relevant phases around the Kitaev point. In the TLL and Spiral XY phases the von Neumann entropy scales with the sizeSL(L/2) ∝ L 1/4 , whereas it oscillates with a periodicity of 4 in the other phases. We also note that (not shown here) the interactions in Eq. 47, controlled by δ, introduce large fluctuations in the finite-size scaling for both the TLL phase and the Spiral XY phase, thereby eventually breaking the TLL phase.
Since we want our mode to commute with the Hamiltonian, we offset this by introducing higher-order terms to the mode Ψ (1) x 2Ŝ y 3 K/(K + δ), so that NowV and Ψ (2) B do not commute, and produces a higher order term V , Ψ In the spiral-XY phase the unit cell still encompasses two sites, while the order parameter has periodicity of four lattice sites. The Majorana edge-mode phase has an antiferromagnetic structure in Sx whereas no particluar structure is detected in the other phases. Inset: Spiral structure of the spiral-XY phase with arrows indicating the direction starting at the first site i = 1. Right panel: Spin structure factorsP (q) = P x (q) + P y (q). The periodicity of the spiral-XY phase is indicated by a large peak at q = π/2. Like in the case of for Ψ A , all terms in Ψ anticommute with each other, as well as with G z , so the normalization constant are the same as for both Ψ A and Ψ B .
B. Zero modes forŜ z 1 As before, we note that Ψ C =Ŝ z 1 commutes withĤ 0 but not withV , giving the commutation Ĥ MEM ,Ŝ z 1 = −i K K + δŜ y 1Ŝ x 2 which we offset by introducing Ψ giving for the next order V , Ψ Again, all terms in Ψ anticommute with each other so the normalization constant are the same as for both Ψ A and Ψ B . Also, each term has an odd number of σ z so Ψ C anticommutes with G x and G y .

XI. SPIN COHERENCE
Here, we derive the influence of zero modes on the dynamical development of single spins in our model. A schematic derivation of the autocorrelation is performed in the main text in section IV B, and this section is included as a complement. We start by noting that and proceed to write down the autocorrelation function Γ y 1 (t) for an eigenstate |S y , with corresponding eigenvalue s y 1 , of the edge spin Pauli matrix σ y 1 along direction y , where |n , |m are eigenstates of the Hamiltonian. Since Ψ A almost commutes with the Hamiltonian H we may divide all energy states into two sectors denoted by positive or negative signs, corresponding positive or negative eigenvalues of Ψ A so that We can now re-write the autocorrelation function for the new indices Γ y 1 (t) = s y For long times t and large system size L, terms with n A = m A add up incoherently while terms with n A = m A add up coherently, so the double sum may be approximated with T 1 = S y |n A − n A − |σ D 1 |n A − n A − |S y T 2 = S y |n A − n A − |σ y 1 |n A + n A + |S y T 3 = S y |n A + n A + |σ y 1 |n A − n A − |S y T A = S y |n A + n A + |σ y 1 |n A + n A + |S y For T 2 and T 3 we exploit the spin-flip operator in Eq. 7 which anticommutes with Ψ A so that G z |n A ± = |n A ∓ and we get T 2 = S y |n A − n A + |S y · n A − |σ y 1 |n A + = S y |n A − n A + |S y · n A − |σ y 1 |n A + + n A + |σ y 1 |n A − 2 = S y |n A − n A + |S y · n A + |G z σ y 1 + σ y 1 G z |n A + /2 = S y |n A − n A + |S y · n A + |{σ y 1 , G z }|n A + /2 T 3 = S y |n A + n A − |S y · n A + |{σ y 1 , G z }|n A + /2 For T 1 and similarly for T 4 we get

XII. SPIN BEATING
Here, we derive the origin for the beating mechanism in the temporal evolution of edge spins, as discussed in section IV C of the main text. See section IV C for definitions and background. We use the same technique as in the previous section, but with an additional trick to derive the short-time dynamics in the autocorrelation function. We start withĤ =Ĥ +Ĥ P whereĤ =Ĥ KH (θ = π/2) Here we could reduce the general perturbation to a pure Ising term with strength δ by puttinĝ V Ising δ =Ĥ P ∀A i = δ but now we keep the more general form of the perturbation. We note that two terms in the Hamiltonian commute Ĥ ,Ĥ P = 0 and we proceed to study the time evolution