Theory of topological spin Josephson junctions

We study the spin transport through a 1D quantum Ising$\mid$XY$\mid$Ising spin link that emulates a topological superconducting$\mid$normal$\mid$superconducting structure via Jordan-Wigner (JW) transformation. We calculate, both analytically and numerically, the spectrum of spin Andreev bound states and the resulting $\mathbb{Z}_2$ fractional spin Josephson effect (JE) pertaining to emerging Majorana JW fermions. Deep in the topological regime, we identify an effective time-reversal symmetry that leads to $\mathbb{Z}_4$ fractional spin JE in the $\it{presence}$ of interactions within the junction. Moreover, we uncover a hidden inversion time-reversal symmetry that it is showed to protect the $\mathbb{Z}_4$ periodicity in odd chain sites even in the $\it{absence}$ of interactions. We also analyze the entanglement between pairs of spins by evaluating the concurrence in the presence of spin currents and highlight the effects of the JW Majorana states. We propose to use a microwave cavity setup (cQED) for detecting the aforementioned JEs by dispersive readout methods and show that, surprisingly, the $\mathbb{Z}_2$ periodicity is immune to $\it{any}$ local magnetic perturbations. Our results are relevant for a plethora of spin systems, such as trapped ions, photonic lattices, electron spins in quantum dots, or magnetic impurities on surfaces.


I. INTRODUCTION
Condensed matter system is an endless playground for emergent exotic phenomena and quasi-particles. In particular, the concept of topological phases associated with the band structure of solids has seen tremendous developments over the past decades [1]. Topological insulators and superconductors are probably among the most scrutinized, notably because they can host Majorana fermions, quasi-particles that are their own antiparticle, which occur as excitations in such materials [2][3][4]. Thanks to non-Abelian statistics, Majorana fermions are crucial ingredients for a functional topological quantum computer: a set of distant, non-interacting Majorana fermions allow, through the process of braiding, to implement a category of topologically protected gates, albeit not universal unless supplemented by other conventional gates [5][6][7].
Genuine topological superconductors are rare, for example, Sr 2 RuO 4 is believed to be one [8]. However, quantum engineering through the use of proximity effects in fermionic systems, can lead to such special superconductors, i.e. 1D nanowire and 2D topological insulator with strong spin-orbit interaction (SOI) proximitized with the conventional s-wave superconductor [9][10][11]. On the other hand, quantum magnets can mimic electronic systems without above proximity requirements [12][13][14]. Specifically, a 1D quantum Ising model can emulate a Kitaev p-wave superconductor, via a renowned Jordan-Wigner transformation (JWT) [15][16][17][18]. In particular, the topological phase transition and the occurrence of Majorana fermions as low-energy modes are all mapped into the spin system when the applied transverse magnetic field * trifmircea@gmail.com is varied, where the ferromagnetic (paramagnetic) phase in the spin chain corresponds to the topological (trivial) phase of the fermionic system [19].
However, one should not be mislead: although there are some analogies of low-energy excitations between fermionic system and spin space, some topological properties will be lost after transformation [20][21][22]. In the spin space, Majorana fermions are not localized objects anymore, and they can be mixed simply by a magnetic field along the Ising axis, i.e. the parity of the ground state is fragile. Nevertheless, it is of crucial importance to investigate how many of topological properties associated with Majorana fermions survive in the spin chain, and provide possible experimental witnesses of their manifestations, such as the existence of Majorana modes. To achieve that, in this work we propose and study the spin transport through an Ising|XY|Ising (IXI) spin link in the presence of Ising axes misalignment. Borrowing from the electronic description, such spin chain system emulates a phase biased topological superconducting|normal|superconducting (SNS) junction that hosts Andreev bound states (ABSs), with supercurrents flowing through the middle part [23][24][25][26].
The symmetries of a system play an essential role in the topological phase classification. Nowadays, noninteracting fermionic systems are classified into ten classes by means of three fundamental symmetries: timereversal symmetry (TRS), particle-hole symmetry (PHS) and sublattice symmetry (SLS) [27][28][29][30]. Besides, crystalline symmetries (e.g. inversion symmetry) [31][32][33][34], as well as many-body interactions [35,36], can also lead to different topological classes which, among other things (e.g. impurities [37][38][39]), may result in various periodicities of Josephson effects (JEs). Roughly speaking, periodicities of JEs in fermionic systems are 2π in the trivial phase, 4π in the topological phase, 8π in topological phase with many-body interactions or impurities (see Sec. V for rigorous descriptions). The latter two cases are known as Z 2 and Z 4 fractional JEs pertaining to contributions from Majoranas and parafermions respectively [39][40][41]. While in this paper, we would like to see how symmetries are transformed between spins and fermions, especially how these fractional JEs affect the spin current transport [42][43][44][45], even under the influence from manybody interactions and quasi-particle poisoning.
One of the most counter-intuitive characteristics in the quantum world is the entanglement, whose non-locality provides another instructive insight to understand topological phases [46]. Nowadays, there is still no universal way to quantify the entanglement of a mixed state shared by arbitrary subsystems [47]. However, one could analytically compute the entanglement of a mixed state in a bipartite spin-1 /2 systems by virtue of concurrence [48]. The variation of entanglement across the quantum phase transition point has been investigated in the anisotropic XY spin chain with periodic boundary conditions [49]. Here instead we evaluate the entanglement between spins in the presence of a spin flow, and analyze how the spin supercurrent pertaining to the Majorana fermions affect the entanglement.
The experimental method of choice for detecting spin currents in insulating (quantum) magnets is via the inverse spin-Hall effect: a metal with strong SOI is coupled to the insulating magnet, and the spin current injected in the metal, via the SOI, results in a charge current that can be measured by usual means [50]. While this method is effective for large spin systems, the signal might be too small for quantum spin chains. Thus, we propose to detect the spin current flow in a cavity QED setup, where such a spin flow affects the cavity frequency and Q-factor which can then be detected by measuring the spectral features of the cavity.
The paper is organized as follows. In Sec. II we introduce the spin system and the model Hamiltonian. There we perform the mapping from spins to fermions via the JWT. In Sec. III we analyze the symmetries of the two representations appearing at the lattice level, and figure out their resulting degenerate properties in the spectrum. In Sec. IV we consider the low-energy continuum theory, both close to the quantum phase transition point as well as in the deep topological regime, to solve for the ABSs spectra analytically, and compare to those found numerically. In Sec. V we discuss different scenarios of fractional JEs regarding effective TRS in the continuum limit and inversion TRS at the lattice level, respectively. In Sec. VI we calculate the texture of spin entanglement quantified by concurrence in the presence of a persistent spin supercurrent in the XY sector. In Sec. VII we propose and analyze the coupling of the spin chain to a microwave cavity for readout of the spin current and the periodicities, along with examining the robustness of fractional JEs under perturbations of in-plane magnetic fields. Finally, in Sec. VIII we end up with some conclusions and an outlook on future directions.

II. MODEL HAMILTONIAN
The N -site anisotropic XY spin chain in a transverse field is presented schematically in Fig. 1(a) with open boundary conditions, whose Hamiltonian readŝ is a spin vector constructed by Pauli matrices at site i, m i = (cos φ i , sin φ i , 0), n i = (− sin φ i , cos φ i , 0), φ i is the spin anisotropic angle with respect to the z axis, 0 ≤ γ i ≤ 1 marks the degree of anisotropy in the xy-plane, J > 0 is the spin exchange constant, 0 ≤ t i ≤ 1 is the coupling strength, g i = g is the relative magnitude of the global transverse field along the z axis, lengths are measured in units of the lattice spacing a.
By tuning the value of parameters, the chain is split into three regions: Ising|XY|Ising, connected with two spin Josephson junctions (JJs). The number of sites in left, middle and right part is N L , N M , N R , respectively. The middle chain, along with two left and right interfaces (x L = N L , x R = N L + N M ), are described as the isotropic XY model by setting γ i = φ i = 0, ∀i ∈ [x L , x R ]. The left and right parts refer to the quasi-Ising (anisotropic XY) model γ i = γ = 0, together with a global twisting angle φ i = φ imposed on the right part, while φ i = 0 in the left part. Although the coupling strength is set as t i = t globally, the parameters at two interfaces can be specifically adjusted to t x L = t x R = t. When t = t, the connections between different parts are perfect, while if t = 0 they are decoupled from each other.
We perform the JWT, (1) and obtain the fermionic Hamiltonian whereĉ † i (ĉ i ) is the creation (annihilation) operator of the JW electron at site i. It turns out the IXI emulates a topological superconducting|normal|superconducting (SNS) junction. Fig. 1(b) presents a schematic of the mapping from the IXI model to the SNS structure. Since Eq. (2) is quadratic, it can be expressed in Bogoliubov-de T is a 2Ndimensional spinor and |i = (0, . . . , 1, 0, . . . ) T is an Ndimensional basis vector corresponding to the i-th site of the chain, ρ x,y,z are Pauli matrices acting on the Nambu particle-hole space. By use of the Bogoliubov quasi- with a set of singleparticle energy n .
When the twisting angle φ of the right Ising part is nonzero, there is a spin supercurrent flowing through the middle sector, whose coupling Hamiltonian is XY-typê Note that Ĵ i remains a constant ∀i ∈ (x L , x R ) in the middle part due to current conservation. After adiabatic evolution of φ for a full period, there can be a net spin pumped between the left and right Ising parts.

III. LATTICE SYMMETRY ANALYSIS
The symmetries of a system are independent of representations, although they can be interpreted differently in the spin and fermionic pictures. In following subsections we will identify the symmetries occurring in the spin system and find out their fermionic counterparts through the JWT. To be more general, we introduce two types of interacting Hamiltonians: the ZZ-typê acting on the spin space and the NN-type (Coulomb) in terms of fermions, which are connected by the JWT as 4n ini+1 ⇔ 1 +σ z i +σ z i+1 +σ z iσ z i+1 . Although they are equivalent after global renormalization of the magnetic field, will have significant implications as we see below.
A. Spin Z2 Symmetry The spin chain contains Z 2 symmetry as [Ĥ S G ,P S ] = 0 by an operatorP S = iσ z i ,P 2 S = +1, which acts on Pauli operators aŝ By the JWT, the corresponding operator in the fermionic system is identified as a parity operatorP F = i (2n i −1), which transforms fermionic operators aŝ Since Eq. (2) is a sum of terms containing even number of fermionic creation and annihilation operators, the system is required to preserve the parity as [Ĥ F G ,P F ] = 0 at any time, although the number of fermions is not conserved. One can easily verify that Z 2 symmetry holds for the aforementioned two types of interacting Hamiltonians. Considering the pure Ising chain with δ i = g i = 0 and γ i = t i , the spin ground states will simultaneously break above Z 2 symmetry, which in turn gives two degenerate ground states in the Kitaev model characterized by Majorana zero modes.
where K is an anti-unitary complex conjugate operator. SinceT 2 S = (−1) N , according to Kramers theorem, all many-body spectra must be at least doubly degenerate when N is odd. Through the JWT, Eq. (2) fulfills [Ĥ F G ,T F ] = 0 inherited from the spin space, This can be interpreted as the charge conjugation in the fermionic language. Based on non-interacting Eq. (3), we can rewriteT F in a first-quantized form which renders [H F G , T F ] = 0. Note that Eq. (10) is more general than Eq. (11) since it can deal with Eq. (6) and find out [Ĥ F I ,T F ] = 0, yet the ZZ-type interactions in Eq. (5) retain rTRS due to [Ĥ S I ,T S ] = 0. When N is odd, the twofold degeneracies in the manybody spectrum are protected by the second-quantized rTRS operator withT 2 F = −1, which enforces intrinsic zero modes in the single-particle spectrum. Under the fermionic picture, as the coupling strength t increases, the amplitudes of aforementioned zero modes in the middle part will exponentially leak into the superconducting parts, and fully merge with Majorana zero modes in the thermodynamic limit, whose wavefunctions are well localized at the edges of the chain and cause no effect on the sub-gap spectrum.

C. Inversion Time-Reversal Symmetry
Although the first-quantized rTRS operator T 2 F = +1 cannot reflect any degenerate properties in the singleparticle spectrum, it gives us a hint to find out a hidden inversion TRS (iTRS) which leads to an odd-even effect (see discussions in Sec. V B). We first introduce a lattice inversion operator which will transform matrix elements of the nearest neighboring sites with additional minus sign after applying on the lattice space, e.g.
we denote parameters with tilde are elements inverted from original position. With the help of I, we can define the iTRS operator Apply By comparing Eq. (3) and Eq. (14), it turns out that in order to retain iTRS as [H F G (φ), T I ] = 0 in the IXI chain, not only should we setg i = g i = 0,t i = t i ,γ i = γ i , but also φ is required to be specific values: with l ∈ Z. Note that T 2 I = −1 in both odd-even cases, according to Kramers theorem, all single-particle states at above specific φ should contain twofold degeneracy.
More generally, we rewrite Eq.(12) as second-quantized form acting on fermionic operators asT I iT −1 and its actions on spin space areT Iσ z iT which can be understood as the charge-parity symmetry. When we deal with ZZ-type interactions in Eq. (5), once given δ i =δ i symmetrically, the system Hamiltonian always commutes with iTRS operator at specific φ illustrated in Eq. (15), which ensures twofold degeneracies of many-body states in the interacting case. As for the NN-type interactions in Eq. (6), we obtain and expand to iχ i (n ini+1 , whose last three terms will break iTRS at any φ, even if we set χ i =χ i symmetrically.
Note that in the above proof all parameters are required to hold strict inversion symmetry under N L = N R , thus the odd-even effect only depends on N M . However, by the fact that the ABSs decay exponentially in the two superconducting parts, as long as their lengths are much larger than superconducting coherence length, the degenerate properties are still robust within the energy gap regardless of the parity and the equality of N L and N R , which in turn underlines the dominance of N M .

IV. LOW-ENERGY CONTINUUM THEORY
In following subsections we will develop the low-energy continuum theory with the aid of fermionic descriptions. Given translation symmetry under periodic boundary conditions, the bulk spectrum of the isolated anisotropic XY spin chain reads [16]: where k is the wave number after the Fourier transformation. When γ = 0 the spectrum is always gapped except at |g| = 2t where the system undergoes a quantum phase transition. In the case of |g| < 2t, the fermionic chain will be in a topological phase where Majorana fermions appear at the edges if we cut off the chain, and the corresponding topological invariant is characterized by the topological winding number W = 1 (see Appx. A for details). However, if |g| > 2t such edges modes will disappear, the chain enters the trivial phase and the value of the topological winding number goes to zero. Fig. 1(c) depicts the wavefunction of the JW Majorana bound state (MBS) in the presence of a phase bias between two superconducting parts. Note that the middle sector is gapped in the trivial regime |g| > 2t, which hinders the occurrence of the supercurrent and makes the chain insulated. Since we are interested in the JEs pertaining to the supercurrent, we will only focus on the topological regime in the whole paper.

A. Near the Critical Point
On account of the long wavelength excitations dominating the low-energy properties near the critical point [51], we can replace the fermionic operators in Eq. (2) by a continuous Fermi field operatorĉ i = √ aψ(x) and expand it to second order in the spatial gradients to obtain the single-particle continuous Hamiltonian: The coefficient in front of the second and first derivative indicates the effective mass m * i = 2 /(4Jta 2 ) and velocity respectively [19]. In order to mimic the imperfect connections between different parts, we introduce a fictitious potential λaδ(x − x L, R ) at two interfaces x L, R with barrier strength λ. When λ → ∞, the three parts of the chain are decoupled from each other.
After solving out the two-component wavefunction , we can impose boundary conditions at the two interfaces to obtain the left and right scattering matrices and the entries are defined as , Ω = t(2t + g) to simplify the expression. The waves at two interfaces only contain different factors caused by the middle wave number K ± M = √ Ω ± Ξ/ta, which is described by scattering matrix S C M = exp(iρ z K ρz M L), L = (N M + 1)a is the length of the middle part. Notice that such wavefunction factor will be canceled out due to Andreev reflection after traveling for one loop, which enforces det(1 − S C M S C R S C M S C L ) = 0 and gives the solvability equation: The energy spectrum of the ABSs is determined by the above equation. In leading order series expansion around zero energy, the spectrum is given by which is plotted in Fig. 2(a) against spectra from the exact continuum and lattice. The explicit wavefunctions and technical details are illustrated in Appx. B.

B. Deep Topological Regime
In the deep topological regime g → 0, the energy gap gap = 2Jγ 4 − g 2 /(t 2 − γ 2 ) → 4Jγ occurs around ±k F = ± arccos(−g/2t)/a ≈ ±π/2a with the proviso of γ t. Accordingly, we can expand the lattice fermionic operator around two Fermi points asĉ , whereψ R, L are right and left mover field operators. We substitute the above transformation into Eq. (2), expand it to the leading order in the spatial gradients and neglect the fast oscillating terms. By defining a continuous Fermi field spinor

the deep topological Hamiltonian can be expressed in the BdG form
where Υ = 2ta sin(k F a) is the effective velocity, ∆ i = 2γ i sin(k F a) is the effective pairing potential [23], τ x,y,z are Pauli matrices acting on the mover space. Note that the phase is globally shifted by π/4 in order to keep ∆ i a real number. The above Hamiltonian emulates JJs created at the edge of a quantum spin Hall (QSH) insulator [52][53][54], and thus we can define an effective TRS (eTRS) there must be spectrum degeneracies at those specific phases due to Kramers theorem.
Owing to [H D G , τ z ] = 0, it is more convenient to decompose the Hilbert space in two τ z eigensectors τ = ±1 and solve H D Applying continuity conditions at two interfaces x L, R on the wavefunctions of each eigensector, we can obtain the left and right scattering matrices One can use the above transcendental equation to solve out the continuum spectrum. Under the low-energy leading approximation, the energy becomes which is plotted in Fig. 2(b) against the exact continuum spectrum and the full lattice spectrum. The index τ indicates the slope of the spectrum as a function of φ: τ = ±1 for downward (upward) branches respectively. In the case of the point contact limit L → 0, Eq. (26) is reduced to E → τ ∆ cos φ [25]. Appx. B presents the explicit wavefunctions and calculation details.

C. Numerical Comparison
With the single-particle spectrum n at hand, one can construct the many-body spectrum E n : the ground state is built with all the negative-energy single-particles filled, the following excited states are obtained by adding the corresponding quasi-particles to the ground state, whose total number characterizes the parity of the system. Fig. 2 displays the exact numerical single-particle and the many-body spectra near the critical point and in the deep topological regime, compared with results from two low-energy continuum models respectively. It is clear that both continuum theories show great agreement with solutions from the numerical lattice model in the singleparticle spectrum (a)(b), which can be interpreted as follows. When the spin chain is near the critical point Ω → 0 with Γ Ω, the energy gap 2J |2t − |g|| will always happen around k = 0, where the long wavelength continuum theory takes charge of the system. While if the chain is in deep topological regime with Γ Ω, the energy gap gap ≈ 4Jγ reaches around two Fermi points ±k F , which supports the validity of the deep topological continuum theory. From the perspective of fermionic language, the superconducting coherence length is defined as ξ = Υ/∆ = ta/γ [23], while the continuum theory requires the coherence length to be much larger than the wave length, i.e. ξ 2π/k F , which also leads to the validity condition γ t. In spite of the excellent accordance in the singleparticle spectra (a)(b), there are only fair agreements between numerical and analytical results in the manybody spectra (c)(d), where we have globally shifted the energies to make the ground state energy zero at φ = 0. Since the many-body spectrum of the low-energy continuum theory can only be constructed by few singleparticle energies of ABSs, the attributions from propagating states outside the gap will not be captured in the analytical continuum theory. Yet it is still worthy to point out the spectrum near the critical point matches better than that in the deep topological regime due to weaker φ-dependence of propagating state energies.

V. FRACTIONAL SPIN JOSEPHSON EFFECT
Historically, the original JE was used to described the supercurrent tunneling through a weak link between the conventional s-wave superconductors, following 2π periodicity of the system Hamiltonian [26]. Nevertheless, JJs designed between the topological p-wave superconductors are predicted to exhibit 4π-periodic supercurrent, a hallmark manifestation for the existence of MBSs [56][57][58]. Notably, a variety of JEs can be inferred from JJs based on the edges of QSH insulators: under TRS and parity conservation, there is a dissipative current varied as 2π periodicity in the presence of a dc voltage bias; once TRS is broken, the current becomes non-dissipative and evolve as 4π periodicity protected by the PHS from MBSs [54]. Furthermore, given TRS with Coulomb interactions [40] or impurities [37][38][39], the current can even be dissipationless with 8π periodicity, while the s z -conserving interactions will lead to dissipation with original 2π periodicity [38]. Such 4π (8π) periodicity is called Z 2 (Z 4 ) fractional JE for the sake of e (e/2) electron charge being transferred in 2π period of the system Hamiltonian, instead of Cooper pairs 2e in the conventional superconductors. However, in Ref. [59] it was shown that such 8π periodicity can be achieved without Coulomb interactions, based on a p-wave superconductor lattice ring interrupted by one weakly coupled normal site. Before analyzing the spin JEs in our setup, we want to make a key observation. The spin twisting angle φ has been mapped into the superconducting phase 2φ, i.e. it was doubled, which makes all periodicities of the fermionic JEs twice as large as the spin JEs. Explicitly, periodicities of trivial, Z 2 , Z 4 JEs become π, 2π, 4π in the spin chain, respectively, compared with 2π, 4π, 8π in the fermionic systems. To avoid confusion, in following discussions, we will use trivial, Z 2 , Z 4 terms to illustrate various JEs in two representations.
Having seen the scheme of JEs in fermionic systems, a question naturally arises: except for the alteration at the phase φ by a factor of two, are there any topological similarities and differences between the fermionic JEs and the spin JEs? In the following subsections, we will investigate various spin JEs from two perspectives: continuum model and lattice level. Moreover, in order to reveal the influence of many-body interactions on spin fractional JEs, we can add ZZ-type interactions into Eq. (1), and NNtype interactions into Eq. (2), respectively, both of them acts only within the middle sector. Note that previous BdG diagonalization fails due to the non-quadratics of the fermionic Hamiltonian, we should apply brute-force diagonalization on a 2 N -dimensional matrix in the spin space, which is limited by the number of sites.

A. Continuum Scenarios
In the low-energy continuum limit, both Eq. (21) and Eq. (25) obey the PHS: {H C G , C C } = 0, C C = ρ x K near the critical point and {H D G , C D } = 0, C D = ρ y τ y K in the deep topological regime, which guarantees the crossings of MBSs and switches the parity of the ground state at φ = π/2. Additionally, as we have shown in Sec. IV, crossings at φ = lπ/2 are protected by the eTRS of Eq. (25) in the deep topological regime, which is indeed equivalent to JJs attached to the edge of QSH insulators. Therefore, adiabatically advancing the spin twisting angle φ will pump each ABS into the bulk and lead to dissipative current with trivial periodicity, as displayed in Fig. 2(b)(d). Nonetheless, when the system is tuned close to the critical point where the eTRS is broken, there are anti-crossings at φ = lπ/2 in Fig. 2(a)(c), with the exception of φ = π/2 crossing protected by Majorana PHS. Under this circumstance, every ABSs are detached from the bulk and give rise to dissipationless spin current with Z 2 periodicity.
Here we briefly discuss the effect of the many-body interactions on the spin JEs. In Fig. 3, we show the manybody spectra in the deep topological regime, taking into account interactions of ZZ-type and NN-type, respectively, both still with the eTRS maintained. Compared with Fig. 2(d), prior fourfold degeneracy at φ = π/2 is lifted in Fig. 3(b) via the Coulomb interactions, a dissipationless Z 4 spin current raises as expected in Ref. [40]. While ZZ-type interactions in Fig. 3(a) can only shift crossings away as opposed to breaking them. On account of the energy levels moving up into the bulk as φ is increased, the spin current retain dissipative with trivial periodicity as the aforementioned non-interacting case. This phenomenon basically resembles QSH JJs accompanied with s z -conserving interactions in Ref. [38].
Although there are small gaps at φ = π caused by finitesize effects (e.g. slowly oscillatory umklapp or Friedel terms), they can be fairly suppressed under the continuum limit as γ → 0 [55].

B. Lattice Odd-Even Effect
The eTRS in the continuum limit requires the transport through JJs to be highly transparent, any imperfect connections t = t are able to break such symmetry and open gaps at the lattice level, which leads to following odd-even effects. As we have proven in Sec. III, there is an iTRS appearing at the lattice level when all parameters are set inverted symmetrically, bringing about different crossing properties for odd-even sites. Especially for all single-particle states illustrated in Fig. 4(a)(b), there must be Kramers pairs at φ = lπ for odd N and φ = π/2 + lπ for even N , according to conclusions in Eq. (15). This remarkable result can be used to exactly attach or detach crossings at specific φ of the many-body spectra, shown in Fig. 4(c)(d), by means of changing the parity of sites. As a consequence, adiabatically following the ground states will eventually lead to Z 2 (Z 4 ) spin currents for the even (odd) sites, as displayed in Fig. 4(e)(f) calculated by Eq.(4). Alternatively, the spin current can be analytically computed by Ĵ z n = −2∂E n /∂φ, if we apply a phase-shifted JWTĉ † i = e −iφ i−1 j=1 (−σ z j )σ + i on the right part and move φ into the tunneling Hamiltonian [20], which gets along with conventional results for the fermionic Josephson current [26].
In addition, our conclusion reveals the unusual Z 4 fractional JE in Ref. [59] is actually protected by the iTRS. In fact, their model Hamiltonian is equivalent to ours for N M = 1 after applying the phase-shifted JWT [60]. The reason why in their case the Z 4 periodicity cannot survive under the Coulomb interactions is that NN-type interactions do not commute with iTRS, whereas ZZ-type interactions do, as it happens in spin chains. Namely, the spectra may be shifted under ZZ-type interactions while crossings are still protected. Therefore, Z 4 spin current originated from iTRS does not depend on whether there are ZZ-type interactions or not. in both odd-even cases. Solid (dashed) lines in the singleparticle spectra are the energies of the particles (holes), solid (dashed) lines in the many-body spectra refer to the even (odd) parity, solid (dashed) circles are crossings protected by the iTRS (PHS), the gaps specified by the arrows are lifted by the imperfect couplings t < t. Specifically: (a) and (b) are the single-particle spectra for NM = 10 and NM = 11 respectively, (c) and (d) are their corresponding many-body spectra, whose single-particle occupations are shown in plot labels, (e) and (f) are spin currents of their lowest two and four states evaluated from Eq.(4), showing Z2 and Z4 periodicities respectively.

VI. TEXTURE OF SPIN ENTANGLEMENT
In this section we evaluate various spin correlation functions in the presence of the spin supercurrent carried by JW Majoranas in the XY sector. Specifically, we are interested in the single spin expectation value p α i ≡ p α i (φ) = σ α i , as well as the spin-spin correlation function p αβ ij ≡ p αβ ij (φ) = σ α iσ β j with α, β = x, y, z. Their knowledge allows us to derive the reduced density matrices for an arbitrary single i and pair of spins ij, respectively, whose expressions are explicitly shown as: Since the Hamiltonian conserves the parity of the system, we can readily infer that p x i = p y i = 0, thus the spin texture has only one non-zero component p z i , along the z-direction. Based on the same arguments, p xz ij = p zx ij = p yz ij = p zy ij = 0 also holds for the two-spin correlators. Nevertheless, p xy ij and p yx ij cannot be zero when there are spin supercurrents flowing through the middle XY part, as a result of nonzero value in Eq. (4). Under this circumstance, regular determinant stratagems used in Refs. [16,19,49] fail to solve out correlators of two arbitrary spins. However, such correlators, together with nonzero p xx ij and p yy ij , can be obtained by computing the Pfaffian of their corresponding well-organized 2k-dimensional skew-symmetric matrices [61,62], where k = |i − j|, and technical details are given in Appx. C.
With all spin correlators at hand, we are able to establish the reduced density matrices, and then evaluate the degree of entanglement in the system. There are two simple ways to quantify the entanglement between two subsystems [49]: (i) a single site and the rest of the lattice; (ii) two arbitrary spins in the chain. For the former case, the entanglement can be calculated via the von Neumann entropy S i (φ) = −tr[ρ i (φ) log ρ i (φ)], assuming the whole chain in a pure state. For the two sites case in a mixed state, the amount of entanglement shared between the spins is quantified by the concurrence C. In particular, for two arbitrary spin-1 /2 sites at the positions i and j in the chain, the concurrence is given by [48]: where the λ k ij are the eigenvalues sorted in descending order of the Hermitian matrix R ij = √ ρ ijρij √ ρ ij with ρ ij = (σ y i ⊗σ y j )ρ * ij (σ y i ⊗σ y j ). The concurrence increases from C = 0 for a separable state to C = 1 for a maximally entangled state. In Ref. [49] it was showed that the single-site entropy, as well as the concurrence between two arbitrary spins peak around the quantum phase transition as a function of the transverse magnetic field. Here we pose a different question: assuming the system in the topological regime, how is the entanglement in the XY sector affected by the presence of spin supercurrents pertaining to a finite twist angles φ?
In Fig. 5 we plot the texture of the spin concurrences as a function of φ for odd-even cases in different regimes, following the ground states in Fig. 4. It is apparent to see that there are two different textures of spin entanglement for odd-even cases depicted in (a)(b), not only evolving with two kinds of periodicities, but also taking peaks (nadirs) at different φ. Such phenomena are due to the fact that through increasing φ, the many-body levels have been shifted to higher values, which makes them more susceptible to upper excited states. On account of finite size effects with open boundary conditions, the entanglement oscillates with frequency ∼ 2k F as site varies [63], which can be enhanced by larger susceptibilities close to anti-crossing points.
Furthermore, by comparing (a)(b) to (c)(d), one might wonder why concurrences near the critical point are less than that in the deep topological regime, since the chain should be more entangled around quantum phase transition. The reason is as follows: in the deep topological regime, only nearest-neighbor concurrences are nonzero, which means the entanglement is well confined in nearestneighbor spins; while as the system approaches the critical point, the entanglement will be spread out into nextnearest-neighbor (and so on) spins [49], which makes the initial nearest-neighbor concurrence decrease.

VII. DETECTION AND ROBUSTNESS
In this section we address the detection of the spin supercurrents pertaining to the JW Majoranas in the Ising|XY|Ising spin junction. While the method of choice for measuring spin currents is through the use of the spin Hall effect [50], in which case a spin current is converted to a charge current that can be measured by usual techniques, via the SOI in the adjacent material, here we propose a less invasive method based on microwave detection. Such an approach has been found suitable for measuring both the statics and dynamics of ABSs in electronic systems [64][65][66]. The idea is to couple the field of a nearby resonator to various observables of the system, the interaction terms reads: where a (a † ) is the annihilation (creation) operator for the photon in the resonator (assuming one mode only), whileÔ are the observables of the system, e.g.Ô =σ α i (or the sum of a string of spins), with coupling strength β. This coupling will alter the properties of the resonator, which in turn can be measured in dispersive readout. Following Ref. [67], we can write the equation of motion for the cavity field in the Heisenberg picture as: whereĤ ph = ω 0 a † a is the cavity Hamiltonian, κ quantifies the decay rate of the cavity and b in (t) is the input field sent to probe it. Note that the output field, exiting from the cavity b out (t), and the input one satisfy b out (t) = b in (t)+ √ κa(t), which is used to infer the cavity response. In leading order in the cavity-system coupling and in the frequency space, we find [67]: where a(ω) = dt e −iωt a(t) and being the retarded correlation function associated with the observableÔ over the stationary state of the system . . . 0 . Above, |n and E n are the many-body eigenstates and eigen-energies of the system, respectively, F n is the many-body occupation, while the index selects only the states n = m in the summation. Note that all quantities are expressed in the interaction picture, and Ô I (ω) 0 is the expectation value of the observableÔ in the frequency space in the absence of the cavity. Since the energies E n , as well as the matrix elements m|Ô|n are functions of φ, the entire correlation function will carry such a dependence too. In typical spectroscopic experiments, the input field b in (ω) Ô I (ω) 0 (large number of photons are sent into the cavity), and we can neglect this term in the following. Nevertheless, such contribution can become relevant in out of equilibrium situation, when it affects the photon number and photon statistics in the cavity. We will not discuss such regimes here, but refer to Ref. [66] for some details (along with the schematic of cQED setups). The effect of the correlation function on the cavity is as follows: the real (imaginary) part renormalizes the cavity frequency (decay rate).
In this work, we consider a capacitive-like coupling between the spin chain and the cavity magnetic field (through the Zeeman coupling), following Ref. [66]. Moreover, we assume the magnetic field of an microwave cavity couples to the spins in the XY part over a length l < L, orÔ =Ŝ l · n, withŜ l = i∈lσ i . Here, n is the direction of the cavity magnetic field at the position of the wire, which can be different from the z direction, and the coupling is assumed to take place from site l 0 to site l 0 + l − 1. The susceptibility can be written as Π S (ω) = Π z S (ω) + Π ⊥ S (ω), where the first and second term corresponding to the matrix element m|Ŝ z l |n (longitudinal) and m|Ŝ l · n ⊥ |n (transverse), respectively, with n ⊥ = n − e z . There are no cross terms between the z (does not change parity) and x, y (flip parity) spin components as all the states in the system have a definite parity. The above susceptibilities have a simple interpretation in the fermionic language: the first contribution stems from the cavity probing particle number operator over the length l, while the second one effectively represents electronic tunneling into the spin chain over the same distance, thus accessing the transport properties of the spin chain. However, as we see in the following discussions, the analogy is only partial for the second coupling because of the non-locality of the JW string.

A. Longitudinal Susceptibility
The longitudinal susceptibility can now be numerically evaluated from the lattice model by including all possible states. However, in order to understand the behavior, it worth analyzing the limit of small ω ∆ in which case the cavity probes mostly the low-energy ABSs (truncated up to the twelfth state in calculation), including the MBSs. We transform the spins into fermions in the latticeĉ i , and eventually in terms of quasi-particles describing the Andreev statesd n , with i and n specifying the lattice and eigen-energy index, respectively. By usingĉ i = n [u n (i)d n + v * n (i)d † n ] with coefficients u n (i) and v n (i) found from wavefunctions of numerical diagonalization (see Appx. A for details), we write downŜ z l in the form of quasi-particles: with a s (i) = u s (i) + v s (i), b s (i) = u s (i) − v s (i), where r, s are single-particle indices of their corresponding manybody states in Eq. (33), given in the labels of Fig. 4(c)(d).
There are two types of m|Ŝ z l |n : quasi-particle conserving type S c r,s and non-conserving type S n r,s , which are shown explicitly as With single-particle occupation f s ≡ d † sds , the longitudinal susceptibility is written in the single-particle form: where the first (second) line accounts for the quasiparticle conserving (non-conserving) contributions.
In Fig. 6 (a)(b) we show the real and the imaginary parts of Π z S (ω) as a function of φ for odd and even cases respectively, evolving adiabatically in their initial ground states at φ = 0, whose peaks indicate the resonances between the cavity and the low-energy levels in Fig. 4(c)(d). They present different periodicities and reach peaks at different φ, as a result of odd-even effect. Particularly, one can distinguish Z 4 spin current from Z 2 case, by dint of opposite signs near φ = π in the imaginary parts. Moreover, even taking into account the relaxation effects such that the system always follows the ground state, the real part still exhibits a singularity at φ = π in Fig. 6(c), which is again a signature for Z 4 crossing of the levels. We note that while the magnetic coupling to each individual spin is typically small (a few Hz in cQED setups), by coupling the cavity to many spins in the chainsŜ l , the response function is enhanced by an order ∼ l 2 as compared to the single spin scenario.

B. Transverse Susceptibility and Spin Noise
Borrowing from the fermionic parity-flipping picture due to quasi-particle poisoning, one may conjecture that the transverse susceptibility Π ⊥ S (ω) has a nonzero value. Surprisingly, we find out numerically that the matrix elements of Π ⊥ S (ω) are exponentially reduced to zero as the length of Ising part increases, which makes transitions between different parities impossible in topological spin JJs. Such phenomenon is because local in-plane spin operatorsσ x i ,σ y i become highly non-local objects with additional JW string in the fermionic space, it is inevitable to alter states of external JW Majoranas, which in turn flips the parity back to itself and thus forbids transitions between them. To verify this, we study the influences from two kinds of in-plane perturbations within the middle part (x L , x R ): where η x i and η y i are perturbation strengths along x and y directions respectively, both set randomly site by site. Eq. (39) indeed emulates the conventional local fermionic perturbations from quasi-particle poisoning and breaks Majorana crossings in Fig. 7(b) as expected. On the other hand, in Fig. 7(a) we see that the local spin perturbations only shift twofold degeneracy (from external JW Majoranas) away and cannot destroy Z 2 periodicity (even we extend the random perturbations to the whole spin chain), in stark contrast to topological JJs in superconducting systems.

VIII. CONCLUSIONS AND OUTLOOK
In this work, we analyzed an Ising|XY|Ising spin link that emulates a topological SNS structure, both analytically and numerically. In a nutshell, our results can be summarized as following comparisons: (i) Odd vs. Even. The iTRS gives rise to the odd-even effect at the lattice level and protect Z 4 (Z 2 ) fractional spin JE in odd (even) chain sites, irrespective of ZZ-type interactions. The resulting texture of spin entanglement highlights the effects of the spin currents carried by JW Majoranas, whose periodicities can be detected by cQED setup through dispersive readout methods.
(ii) Lattice vs. Continuum. By use of the low-energy continuum theory, we analytically solve out the spectra of ABSs and their fermionic wavefunctions. Particularly in the deep topological regime, there can host eTRS that mimics QSH JJs, forming Z 4 fractional spin JE in the presence of Coulomb interactions within the junction, while ZZ-type interactions act like s z -conserving interactions which leads to dissipative trivial currents.
(iii) Spin vs. Fermion. At the lattice level we identified various symmetries emerging from the spin chain and figure out their electronic counterparts, demonstrating that ZZ-type interactions and NN-type interactions act differently in the many-body spectra. One remarkable result is that although Z 2 -periodic current can be broken by local fermionic perturbations, spin Z 2 JEs are robust to local spin perturbations.
Our proposal could be implemented in a plethora of spin systems, such as trapped ions [68], photonic lattices [69,70], electron spins in quantum dots [71], and magnetic impurities on surfaces [72,73]. Since Z 2 fractional spin JEs are immune to any local perturbations from arbitrary directions of magnetic field (as long as the chain is still in the topological phase), the ground state can, together with the first excited state, be used to set up a logical qubit: advancing φ adiabatically for π achieves a quantum X gate [20,44]. Alternatively, we can utilize such robustness for quantum memory. Besides, the middle XY chain will be gapped when |g| > 2t, which prohibits the transport of spin supercurrents. Hence, one may use this feature to engineer a quantum spin transistor based on the JEs [42].
There are several generalizations of our study. First, it would be interesting to consider dissipation (due to, for example, the presence of a magnetic substrate), and evaluate its effects on the various fractional JEs, as well as on the topology of the chain in general. Moreover, the cQED setup proposed here could serve as an engineered environment that can not only monitor the spin flow, but also affect and control it. Second, generalization to multi-quantum spin chains junctions, similar to superconducting systems [74], which could result in emulating various higher dimensional topological structures. And third, generalization to more complex insulating quantum spin systems, such as 2D quantum (anti) ferromagnets insulators or even quantum spin liquids [75], subject to dissipationless spin flows.

IX. ACKNOWLEDGMENTS
We thank Peter Zoller, Dong-Ling Deng, Thore Posske, Tie-Cheng Guo for helpful discussions. This work was supported by the NSFC under the Research Fund for International Young Scientists No.11750110412 and the International Centre for Interfacing Magnetism and Superconductivity with Topological Matter project (MT), carried out within the International Research Agendas program of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund.
Since we are only interested in the ABSs, whose eigenvalues lie within the gap, i.e. | | < 2J(2t + g) ⇔ |Ξ| < Ω, which ensures K ± M to be real. By introducing a new set of coefficients C 5 , C 6 , C 7 , C 8 in the middle region, the explicit wavefunctions are shown as The above wavefunctions have been normalized by the square root of wave numbers in order to maintain the quasi-particle currents [24]. Through imposing continuity conditions and current conservation conditions at two interfaces presented in Appx. B 3, we obtain energy transcendental Eq. (23) for the ABSs, and their wavefunction coefficients are determined by normalization condition |u n (x)| 2 + |v n (x)| 2 dx = 1. Now the Hamiltonian is diagonalized into n n (d † ndn − 1 /2) by Bogoliubond n , whose transformation with field operator is given bŷ It is worthy to point out that it is the branches jumping of V on the Riemann surface that takes great effect on the quantum phase transition, i.e. V → − arccos[(Λ + Ξ)/Γ]/2 with additional minus sign across the critical point, which prohibits the zero mode solution of Majoranas.

Wavefunctions in the Deep Topological Regime
By decomposing the Hilbert space into the two τ z eigensectors τ = ±1, we can solve out the continuous Eq. (25) in the deep topological regime as H D G Φ τ (x) = τ Φ τ (x) with their corresponding eigenfunctions Φ + (x) = [u + (x), 0, v + (x), 0] T , Φ − (x) = [0, u − (x), 0, v − (x)] T , whose explicit expressions are shown as where K = √ ∆ 2 − E 2 /Υ, W = arccos(E/∆)/2, E = /2J are introduced for simplicity. The wavefunctions of left and right parts only contain the exponential decaying branches due to infinite boundary conditions, while the middle part is the case of φ = γ = 0 where K = iE/Υ ≡ iK M . Since gap → 4Jγ, ∆ → 2γ in the deep topological regime, |E| < ∆ will be always valid for the ABSs. The explicit middle wavefunctions are shown as u τ M (x) = C 3 exp (+iτ K M x) , v τ M (x) = C 4 exp (−iτ K M x) with two new coefficients. Applying continuity conditions at two interfaces for each eigensector respectively, we can obtain the energy transcendental Eq. (26) for ABSs. By use of the normalization condition, the full normalized wavefunctions for the whole chain are expressed as where A n = 1/ 2(L + 1/K) is the normalization factor, l(x) = x for x ≤ |L/2| and sgn(x)L/2 for x > |L/2|. With the help of Bogoliubond τ n , Eq. (25) is diagonalized into n,τ τ n (d τ † nd τ n − 1 /2) with transformation: Apart from PHS as {H D G , C D } = 0, C D = ρ y τ y K, Eq. (25) contains eTRS as [H D G (φ), T E ] = 0, T E = iτ y K for φ = lπ/2, l ∈ Z. Since T 2 E = −1, the single-particle spectrum must be at least doubly degenerate at above phases.

Boundary Conditions Near the Critical Point
We can add a fictitious barrier potential λaδ(x − x ± ) into Eq. (21) to emulate the imperfect connections between different parts (we denote +, − for R, L respectively to generalize the expressions of two junction sites in the following statements). Around two interfaces, the stationary Schrödinger equation requires: where the phase φ is absorbed in γ temperately, the anticommutator parentheses [Θ(±x ∓ x ± ), ∂ x ] + can be calculated into 2Θ(±x ∓ x ± )∂ x ± δ(x − x ± ). Move the second order derivative term to the left hand side and integrate the whole equation around the junction sites by an infinitesimal parameter, we find Replacing subscript +, − back into R, L and specify the value of γ, φ in different parts (releasing φ from γ), we obtain the current conservation conditions: of the following well-organized 2k-dimensional skewsymmetric matrices [61,62]: with their corresponding blocks The correlators of σ x iσ y j , σ y iσ x j only differ on the last operator from σ x iσ x j , σ y iσ y j , hence we can calculate Q xy ij , Q yx ij by replacing the last column and its corresponding transpose row · · · of Q xx ij , Q yy ij : · · · G j−1,j · · · M i+1,j · · · . . .
When the twisting angle is zero, the spin supercurrent vanishes with σ x iσ y j = σ y iσ x j = 0. Furthermore, blockdiagonal terms in Q xx ij , Q yy ij are also found out to be zero. In this special case, σ x iσ x j and σ y iσ y j are reduced into det(G xx ij ) and det(G yy ij ) respectively, which agree with previous formulae used in Refs. [16,19,49].