On the conclusive detection of Majorana zero modes: conductance spectroscopy, disconnected entanglement entropy and the fermion parity noise

Majorana zero modes (MZMs) [1–4] in condensed matter systems have been a crucial talking point, tending to their exotic attributes, such as non-Abelian quantum statistics, which maybe leveraged for futuristic quantum technologies like fault-tolerant topological quantum computation [5]. Semiconductor–superconductor heterostructures composed of Rashba nanowires proximitised with superconductors are a prime solid-state candidate for realizing MZMs [1–4, 6]. Howbeit, conclusive and unambiguous detection of MZMs in these systems has been a matter of intense scrutiny and open deliberation due to controversies involving non-topological origins of MZM-like signatures. Major experimental efforts for the detection MZMs have been based on observing zero-bias conductance peaks (ZBCPs) in the local differential conductance of normal-topological superconductor (N-TS) setups [7–12]. However, disorder-induced near zero-energy Andreev bound states (ABS), often known as quasi-MZMs, imitate local conductance signatures of MZMs to a great extent [13–21], rendering ZBCP based detection of MZMs ambiguous.

Recent protocols have thus focused on exploiting the non-locality of MZMs through non-local transport spectroscopy using the normal metal-topological superconductor-normal metal (N-TS-N) hybrid setups [6, 22–28]. The non-local differential conductance probes the bulk-gap closing and reopening and is presumed to be uninfluenced by the trivial modes [23]. Ideas proposing measurements of the full conductance matrix, which includes both local and non-local measurements, to identify topological regimes have also been presented [29, 30]. Nonetheless, numerical [28, 31] and experimental [25] studies have portrayed several issues with the non-local conductance signatures as well, with the signals becoming faint and sputtered [28] accompanied by a premature gap closing [25, 31] in the presence of disorders. This motivates exploring other possibilities that exploit non-local correlations to classify topological orders. One such attribute which can be leveraged to ascertain non-locality in topological phases [32–34] is the degree of entanglement in the system, which may be quantified using metrics derived from the von Neumann entropy, such as the disconnected entanglement entropy (DEE) [31].

In this work, we provide an in-depth analysis of various entanglement entropy signatures of MZMs in SM-SC heterostructures. We illustrate that the DEE [35] distills out the topological contributions of $ln2$ units to the entanglement entropy [35–37] and is hence apt for identifying topological phase transitions, unlike the bipartite entanglement entropy (BEE) which is padded with volume and area contributions [31, 38]. In the process, we develop a new pedagogical model, 'the spinfull Kitaev chain', to elucidate the entropy signatures of $ln2$ units in Rashba nanowires from a more first-principles perspective. We demonstrate that the entanglement entropy remains robust and quantized over a wide range of controllable parameters and withstands the test of disorder and quasi-MZMs. Further, we inquire into fermion parity noise, a physical observable closely related to entanglement entropy [39]. With recent advances in the experimental measurement of fermion parity of SC-SM heterostructures [40–42], the fermion parity noise maybe potentially used for gauging entanglement entropy. We conclude by demonstrating the reliability of the DEE over the local and non-local conductance spectroscopy [31] for the conclusive detection of MZMs.

In this section we briefly highlight the system Hamiltonian that we consider. We then elucidate the techniques used for calculating the three major quantities presented in this paper—bi-partite and DEE, fermion parity noise, and local and nonlocal conductance.

The three terminal N-TS-N setup [25, 29, 30] that we consider is shown in figure 1(a). The setup comprises a semiconducting nanowire (in blue) with strong Rashba type spin–orbit coupling with an epitaxial layer of s-wave superconductor (in purple), placed in contact with metallic leads 'L' and 'R' (in pink). The Hamiltonian of the nanowire system in the Bogoliubov–de Gennes (BdG) representation is

$$\begin{equation} \begin{aligned} &\mathcal{H} = \frac{1}{2}\int {\mathrm {d}}x \Psi_{\mathrm {BdG}}^{\dagger}\left(x\right) H_{\mathrm {BdG}} \Psi_{\mathrm {BdG}}\left(x\right)\\ \end{aligned} \end{equation} \tag{ 1 }$$

where $\Psi_{\mathrm {BdG}}(x) = (\psi_\uparrow(x), \: \psi_\downarrow(x), \psi_\uparrow^\dagger(x), \: \psi_\downarrow^\dagger(x))^T$ represents the Nambu spinor, ψ_σ is the electron annihilation field operator with spin $\sigma = \{\uparrow, \downarrow\}$ and $H_{\mathrm {BdG}}$ represents the BdG Hamiltonian [18] described by:

Zoom In Zoom Out Reset image size

$$\begin{equation} \begin{aligned} H_{\mathrm {BdG}} = \left(-\frac{\hbar^2}{2 m}\partial_x^2 -\mu + V\left(x\right)\right) \tau_z \sigma_0+i \alpha_{\mathrm {R}} \tau_z \sigma_y \partial_x + V_z\tau_z \sigma_x+\Delta \tau_y \sigma_y, \end{aligned} \end{equation} \tag{ 2 }$$

where, σ_i are the set of Pauli-matrices in the spin basis, with σ₀ being the identity matrix and τ_i are the set of Nambu-matrices in the particle-hole basis. µ is the electrochemical potential, $\alpha_{\mathrm {R}}$ is the Rashba spin–orbit coupling strength, $V_z = \frac{g\mu_{\mathrm {B}} B}{2}$ is the Zeeman field, B is the magnetic field and Δ is the proximity-induced superconducting gap. In this work, we elucidate two cases. First, the pristine nanowire case with $V(x) = 0$, portrayed in figure 1(b). Here, both the superconducting gap and the on-site potential are homogeneous and constant throughout the nanowire. As stressed upon earlier, accounting for disorder is crucial for studying conclusive MZM-identification techniques. In the second case, we study the one with disorder where V(x) is non-zero and inhomogeneous. Figure 1(c) shows the disordered setup where the gap-parameter remains homogeneous but the on-site potential varies along the nanowire due to disorder-induced local inhomogeneities.

We analyze smooth inhomogeneous potentials at the ends of the nanowire [28, 43, 44]. In experiments, such disorders occur as Schottky barrier formations due to Fermi energy mismatches with the contacts [45] or due to charge inhomogeneities in the environment or the leads [9, 46]. For our numerical simulations, we model this smooth disorder as a half-Gaussian at each end of the nanowire with a peak value of $V_{\mathrm {max}}$ and variance σ. For our numerical calculations, we discretize (see appendix A) [47, 48] the Hamiltonian (1) which yields a tight-binding Hamiltonian:

$$\begin{equation} \begin{aligned} \mathcal{H} = \sum_i \mathcal{C}_i^\dagger \left[\left(2t -\mu + V_i\right) \tau_z \sigma_0 + V_z\tau_z \sigma_x+\Delta \tau_y \sigma_y \right]\mathcal{C}_i + \sum_{\langle ij \rangle} \mathcal{C}_i^\dagger \left[-t\tau_z \sigma_0+i t_{SO} \tau_z \sigma_y \right]\mathcal{C}_j, \end{aligned} \end{equation} \tag{ 3 }$$

where, $\mathcal{C}_i = (c_{i\uparrow}, c_{i\downarrow}, c_{i\uparrow}^\dagger, c_{i\downarrow}^\dagger)$ is the Nambu spinor at site 'i' with $c_{i\uparrow(\downarrow)}$ being the annihilation operator at site 'i' with spin $\uparrow(\downarrow)$, $t_{SO} = \alpha_R/2a$, $t = \hbar^2/2ma^2$, $V_i = V(x = ai)$ and a is the lattice spacing. A direct diagonalization of the tight-binding Hamiltonian (3) is used to obtain the eigenspectra.

Entanglement entropy: For calculating the entanglement entropy, we leverage its intricate connection with the correlation matrix (see appendix B) [32], which is defined as :

$$\begin{equation} C_{n m} = \left(\begin{array}{cc}\left\langle c^{\dagger} c\right\rangle & \left\langle c^{\dagger} c^{\dagger}\right\rangle \\ \langle c c\rangle & \left\langle c c^{\dagger}\right\rangle\end{array}\right), \end{equation} \tag{ 4 }$$

where, $\langle c^\dagger c \rangle_{nm} = \langle c_{n}^\dagger c_{m} \rangle$, $\langle c^\dagger c^\dagger \rangle_{nm} = \langle c_{n}^\dagger c_{m}^\dagger \rangle$, $\langle c c^\dagger \rangle_{nm} = \langle c_{n} c_{m}^\dagger \rangle$, $\langle c c \rangle_{nm} = \langle c_{n} c_{m} \rangle$ are sub-matrices with $n,m$ being matrix indices representing combined spin and site degree of freedom i.e. $n(m) = i\sigma(j\sigma^{^{\prime}})$ where $i,j$ are sites and $\sigma,\sigma^{^{\prime}} = \uparrow\downarrow$. The matrix elements are calculated using the expectation values of fermionic operators in the BCS ground-state. The truncated correlation matrix of any sub-region of the system, say A with M sites, can be defined as:

$$\begin{equation} \tilde{C}_{n m} = C_{\left(nm\right)\in A}. \end{equation} \tag{ 5 }$$

With a comprehensive calculation (see appendix B), we can prove that the von Neumann entropy of sub-region A can be expressed in terms of the eigenvalues ξ_k of the truncated correlation matrix as:

$$\begin{equation} S = -\sum_{k = 1}^{2M}\left[\xi_k \text{log} \xi_k+\left(1-\xi_k\right) \text{log} \left(1-\xi_k\right)\right] \end{equation} \tag{ 6 }$$

where k runs over the eigenvalues of the truncated correlation matrix (the summation runs till 2M because the spin doubles the degrees of freedom). When the system is partitioned into two subsystems, the von Neumann entropy of one of the sub-systems is called the bi-partite entanglement entropy (BEE). In this work, we calculate the BEE by partitioning the system into two equal halves, as shown in figure 2(a).

Zoom In Zoom Out Reset image size

The DEE [31, 35] is computed using linear combinations of von Neumann entropy of disconnected regions of the system, an idea similar to the topological entanglement entropy in 2D-systems [49]. The disconnected regions used in this work can be seen in figure 2(b), we call the green region A and the orange region B. The DEE S_D is then calculated as:

$$\begin{equation} S_D = S_A + S_B - S_{A \cap B} - S_{A \cup B} \:\:, \end{equation} \tag{ 7 }$$

where, S_A , S_B , $S_{A \cap B}$, $S_{A \cup B}$ are the von Neumann entropies of the corresponding regions. The DEE dissects out area and volume contributions [31, 33, 35] to the entanglement entropy hence distilling out only the topological contributions to it. As we will see further, BEE signatures are plagued with the issue of volume and area terms but DEE signatures are clean and quantized [31, 35].

Fermion parity noise: We also relate the entanglement entropy to a well known conserved physical observable of the system—Fermion parity [39]. The parity operator is given by $\mathcal{P} = (-1)^{\sum_{i\sigma} c_{i\sigma}^\dagger c_{i\sigma}}$. The fermion parity noise is the statistical measurement uncertainty in $\mathcal{P}$ of the BCS ground-state, defined by:

$$\begin{equation} \langle P \rangle_2 = \langle \mathcal{P}^2 \rangle - \langle \mathcal{P} \rangle^2 \end{equation} \tag{ 8 }$$

where, the expectation value is calculated in the BCS ground-state. It can be shown that the fermion parity noise of any sub-region A of a superconducting system with M sites can be related to the eigenvalues of the truncated correlation matrix [39] of that region by:

$$\begin{equation} \langle P \rangle_2 = \sum_{k = 1}^{2M} \xi_k\left(1 - \xi_k\right) \end{equation} \tag{ 9 }$$

where again ξ_k are the eigenvalues of the truncated correlation matrix and the summation runs till 2M to account for spin degree of freedom. We then define disconnected fermion parity noise $\delta P$ on the lines of the DEE, as a linear combination of fermion parity noise of distinct disconnected regions of the nanowire, given by:

$$\begin{equation} \delta P = \langle P_A \rangle_2 + \langle P_B \rangle_2 - \langle P_{A \cap B} \rangle_2 - \langle P_{A \cup B} \rangle_2. \end{equation} \tag{ 10 }$$

Conductance: The local and non-local conductance signatures are calculated using the Keldysh non-equilibrium Green's function (NEGF) formalism [31]. This requires evaluating the retarded (advanced) Green's function $G^{r(a)}$ of the system by incorporating the contacts as self energies $\Sigma_{L/R}$ along with the device Hamiltonian $H_{\mathrm {BdG}}$. In this work, we use the wide-band approximation to define the broadening matrices $\Gamma_{L/R}$ and subsequently the self energies $\Sigma_{L/R}$ with coupling energy $\gamma = 0.01 \Delta$. The retarded (advanced) Green's functions along with the broadening matrices are then used to calculate the energy resolved Andreev transmission ($T_A(E)$), the crossed-Andreev transmission ($T_{CAR}(E)$), and the direct transmission ($T_D(E)$) (For details see appendix C), which are consequently used to calculate the local and non-local conductance in the zero-temperature limit using the Landauer–Büttiker formalism as:

$$\begin{equation} \begin{array}{r} \left.G_{L L}\left(V\right)\right|_{T \rightarrow 0} \equiv \frac{e^2}{h}\left[T_A\left(E = e V\right)+T_A\left(E = -e V\right)+T_{C A R}\left(E = e V\right)+T_D\left(E = e V\right)\right] \end{array} \end{equation} \tag{ 11 }$$

$$\begin{equation} \left.G_{L R}\left(V\right)\right|_{T \rightarrow 0} \equiv \frac{e^2}{h}\left[T_D\left(E = e V\right) - T_{C A R}\left(E = -e V\right)\right], \end{equation} \tag{ 12 }$$

where V is the bias voltage applied at either contact in the three terminal configuration [26].

We present the results obtained from numerical simulations of the system described in figure 1 based on the experimentally relevant system parameters presented in table 1. We start by calculating the entanglement entropy signatures for the pristine nanowire. Figure 3 focuses on the entanglement entropy as a function of the Zeeman field for short ($4 ~\mu$m) and long ($8 ~\mu$m) nanowires.

Zoom In Zoom Out Reset image size

Table 1. Parameters used in the analysis, unless otherwise stated.

Parameter		Value
Effective mass	$m^*$	$0.015 m_{\mathrm{e}}$
Induced order parameter	Δ	$0.25 \ \mathrm{meV}$
Tight-binding hopping parameter	t0	$10 \ \mathrm{meV}$
Rashba spin–orbit coupling	$\alpha_{\mathrm{R}}$	$20 \ \mathrm{meV} \mathrm{nm}$
Chemical potential	µ	$0.5 \ \mathrm{meV}$
Magnetic g-factor	g	40
Critical Zeeman field	$V_z^c$	$0.56 \mathrm{~T}$

We start by evaluating the BEE signatures for the $4 ~\mu$m and $8 ~\mu$m nanowires. The BEE shows no signs of quantization and exhibits a minimal change in its value over the phase transition for both $4 ~\mu$m and $8 ~\mu$m nanowires, as shown in figures 3(a) and (b), respectively, where the red dashed lines mark the critical Zeeman field for the topological phase transition. The critical Zeeman field $V_z^c$ is evaluated using the well-known relation: $V_z^c = \sqrt{(\mu^2 + \Delta^2)}$. The non-zero values of the BEE in the trivial phase and lack of quantization can be attributed to the area and volume contributions that arise from the bulk of the nanowire over the topological contributions from the ends. These contributions render it ineffective for the identification of the topological phase transition.

Figures 3(c) and (d) present the DEE signatures for the $4 ~\mu$m and $8 ~\mu$m nanowires respectively. The red dashed line, as earlier, represents $V_z^c/\Delta$. Unlike the BEE, the DEE remains zero in the trivial regime $V_z\lt V_z^c$. At the critical point, the DEE shows a sudden jump and attains a non-zero quantized value of $ln2$ in the topological regime $V_z\gt V_z^c$. The DEE signature is robust to the system size and shows similar behavior for both short and long nanowires.

The quantization towards the value of $ln2$ becomes more prominent for larger system size, as is clear from figure 3(d). We also observe other prominent features of the SM-SC hybrid system in the DEE signatures. The oscillations in the DEE, as seen in figure 3(c), correspond to the well-known parity crossings of the MZMs near zero energy [50]. The amplitude of the parity crossingsreduces as the system size increases [47]. This aligns with the reduced oscillatory behavior of the longer nanowire in figure 3(d) than the shorter nanowire in figure 3(c). The near-zero value of the DEE in the trivial regime and a quantized non-zero value in the topological regime, irrespective of system size, indicates its capability to faithfully detect the topological phase transition in the nanowire and distinguish the topological phase from the trivial phase.

In figures 3(e) and (f), we highlight the robustness of DEE as a metric to distinguish between the topological and trivial phases and detect a topological phase transition by presenting the topological phase diagram for the system. Figure 3(e) depicts the value of DEE versus an experimentally relevant range of the electrochemical potential and the Zeeman field for a pristine wire of length 4 $~\mu$m. The dashed black curve denotes the theoretically evaluated topological phase boundary for a fixed value of Δ and is defined by the relation $V_{z}^2 - \mu^2 = \Delta^2$. The region to the right of the curve represents the theoretically expected topological phase. On calculating DEE for each value of µ and V_z in this range, we observe that DEE attains a non-zero quantized value of $\text{ln} 2$ in the region to the right of the dashed curve, whereas it stays close to zero to the left of the dashed curve. This shows that the transition in the value of DEE exactly matches the theoretically expected topological phase boundary, and the value of DEE shows a quantization in the topological regime for a wide range of experimentally relevant parameters. A precisely similar feature is observed for a longer nanowire of length $8 ~\mu$m, as shown in figure 3(f), highlighting the robustness of DEE as a metric to the system size.

A DEE quantization value of $ln2$ units is a trademark signature of MZMs [31]. In appendix D, we elucidate upon topological phase transitions in the Rashba nanowires from a much more elementary and first-principles perspective. We highlight the parallels between the Rashba nanowire and a toy model consisting two distinct p-wave superconductors interacting through an s-wave pairing potential by demonstrating that the nanowire Hamiltonian can be decomposed into this new toy model with appropriate basis transformations. We term this toy model as the spinfull Kitaev chain. We then shed light on how the interplay of various interactions in this toy model lead to the trivial and topological phases. Extending the well established entanglement entropy signatures in topical p-wave superconductors like the Kitaev chain [31] we explain how the nanowire system effectively hosts two-MZMs, thus leading to a DEE value of $ln2$ units in the topological phase.

Next, we focus on the effects of disorder in the nanowire system and its entanglement entropy. As mentioned earlier, the disorder that we study in this work, i.e. smooth inhomogeneous potentials at the ends of the nanowire, causes the appearance of near zero-energy sub-gap states called Andreev bound state (ABS) or quasi-MZMs. These states are trivial in nature and the near zero-energy behavior of these trivial states over wide ranges of V_z and µ hinders the conclusive detection of MZMs through present protocols.

First, in figure 4, we show the energy eigenspectrum of the nanowire in the presence of smooth potentials V(x) at the ends. Figures 4(a) and (b) show the energy eigenspectrum for a 4$~\mu$m and an 8$~\mu$m wire, respectively. As visible in in figure 4(a), the system hosts disorder induced near zero-energy modes prior to the topological phase transition (highlighted by the dashed red line). These states become more prominent as the system length increases, as is visible in the longer nanowire case (figure 4(b)). In figures 4(c) and (d), we show the DEE of the 4$~\mu$m and an 8$~\mu$m disordered nanowires corresponding to the energy spectra in figures 4(a) and (b). For both the lengths, the DEE is zero in the trivial phase and rises to a quantized value of $ln2$ only in the topological phase, similar to the pristine nanowire case. Hence, the DEE remains faithful in signaling MZMs and the topological phase transition with clarity in cases where near-zero energy trivial states are present. We also plot the phase diagram of DEE in the disordered case. Figures 4(e) and (f) depict the phase diagram of DEE for a 4$~\mu$m and an 8$~\mu$m wire, respectively, and highlight that the transition in DEE exactly matches the theoretically estimated topological phase boundary. Thus, the DEE conclusively detects topological phase transitions, even in the presence of disorder, for a wide range of parameters V_z and µ.

Zoom In Zoom Out Reset image size

As mentioned before, we now focus on the behavior of fermion parity noise of the nanowire system. Both the fermion parity noise and entanglement entropy can be derived from the eigenvalues of the ground-state correlation matrix. As evident from equations (6) and (9) the fermion parity noise and entanglement entropy share a closely resembling mathematical form and are expected to behave similarly. Therefore, the disconnected fermion parity noise $\delta P$ is a physical observable that can possibly used to gauge DEE in experiments. In figures 5(a) and (b), we plot the disconnected fermion parity noise $\delta P$ of a pristine $4 ~\mu$m nanowire and a pristine $8 ~\mu$m nanowire against the Zeeman field V_z . For both nanowire lengths, we observe that the fermion parity noise is zero in the trivial regime and rises to a values of 0.5 in the topological regime. The signatures closely match the DEE signature, the only change being the quantization value, thus pointing out the connection between DEE and the fermion parity noise.

Zoom In Zoom Out Reset image size

To study the robustness of this physically observable under disorder, we show the corresponding plots for the disordered nanowire case in figures 5(c) and (d). As expected, the signatures stay robust even in the presence of trivial near-zero energy modes and closely match the DEE signatures. We also plot the phase diagrams of the disconnected fermion parity noise for figure 6(a) a $4~\mu$m pristine nanowire, figure 6(b) an $8~\mu$m pristine nanowire, figure 6(c) a $4~\mu$m disordered nanowire, figure 6(d) and an $8~\mu$m disordered nanowire to show that the fermion parity noise signature remains clean and robust over a good range of controllable parameters µ and B hence acting as a faithful detector of topological phase transitions.

Zoom In Zoom Out Reset image size

Recent experimental progress in the measurement of fermion parity of SC-SM heterostructures includes various techniques. Circuit-QED based methods [41, 42] involve coupling the desired system to a superconducting resonator and then a fermion parity readout is enabled through the microwave spectroscopy of transitions between the even and odd parity states. An alternative method is parity to charge conversion [40] where the system is tuned such that the parity of the sub-gap states correspond to the charge occupancy of a tunnel-coupled quantum dot. The parity of the system is thus converted to charge of the quantum dot and a measurement of the quantum dot's charge state then provides an indirect parity measurement. The challenge that lies in realizing fermion parity noise as a robust metric for MZM detection is the precise measurement of fermion parity itself. For once the parity is measured accurately, the fermion parity noise is just the statistical variance of the parity measurements. The aforementioned experimental methods can be conveniently integrated into nanowire device setups and facilitate fast electrical readout of the parity, thus helping mitigate the challenge of accurate fermion parity noise measurement.

Lastly, we find it worthwhile to demonstrate the superiority of entanglement entropy (and hence observables such as the fermion parity noise) over techniques with current experimental attention, such as local and non-local differential conductance. In figure 7, we depict the local conductance (G_LL ), non-local conductance (G_LR ), and the DEE (S_D ) of a three-terminal N-TS-N setup involving an intermediate length $6 ~\mu$m nanowire with minimal to modest disorder. Even under low disorder, the local conductance fails to detect the phase transition with a premature closure, as shown in figure 7(a). The non-local conductance signatures, on the other hand, become obscure and faint under experimentally relevant measurement precision, unable to pinpoint the critical field as depicted in figure 7(b). The DEE, however, remains faithful and accurately signals the phase transition, as shown in figure 7(c).

Zoom In Zoom Out Reset image size

We proposed and analyzed in depth the use of entanglement entropy and its close relation with topological order for a fool-proof signaling of MZMs in Rashba nanowires. We showed that although the BEE is incapable of determining topological phase transitions, owing to volume and area contributions, the DEE reliably pinpoints topological phase transitions following the distillation of the size contributions. The DEE signatures stay robust to system size and controllable parameters such as chemical potential and magnetic field, even in the presence of disorders such as smooth inhomogeneous potentials at the ends of the nanowire relevant to experimental setups. We explained the quantization of the entanglement entropy in a bottom-up fashion by deconstructing Rashba nanowires as interacting Kitaev Chains. Further, we connected the entanglement entropy signatures to a conserved observable of the system, the fermion parity noise, and illustrated that it displays signatures similar to the entanglement entropy. With active advancement in parity measurement of SC-SM heterostructure, the fermion parity noise emerges as a novel metric for MZM detection. Finally, we concluded by elucidating a simple case with minimal disorder and moderate nanowire lengths where the local conductance signature closes prematurely, and the non-local conductance signature is obscure, but the DEE signature is conclusive and faithfully signals the topological transition. The theoretical methods developed here are of significant relevance toward making experimental progress for the conclusive detection of MZMs, and looking beyond conductance spectroscopy measurements.

The authors acknowledge the Visvesvaraya PhD Scheme of the Ministry of Electronics and Information Technology (MEITY), Government of India, as implemented by Digital India Corporation (formerly Media Lab Asia). This work is also supported by the Science and Engineering Research Board (SERB), Government of India, Grant No. MTR/2021/000388 under the MATRICS scheme, the Ministry of Human Resource Development (MHRD), Government of India, Grant No. STARS/APR2019/NS/226/FS under the STARS scheme.

All data that support the findings of this study are included within the article (and any supplementary files).

The Rashba nanowire Hamiltonian is given by

$$\begin{equation} \mathcal{H} = \int {\mathrm {d}} x \Psi^{\dagger}\left(x\right) H \Psi\left(x\right) \end{equation} \tag{ A1 }$$

where $\Psi(x) = (\psi_\uparrow (x), \psi_\downarrow (x), \psi_\uparrow^\dagger (x), \psi_\downarrow^\dagger (x))^T$ is a Nambu spinor and the H is the Hamiltonian density given by

$$\begin{equation} H = \frac{p_x^2}{2 m}\tau_z-\mu\tau_z-\frac{\alpha_{\mathrm {R}}}{\hbar}p_x\tau_z\sigma_y +V_z \tau_z\sigma_x + \Delta\tau_y\sigma_y \end{equation} \tag{ A2 }$$

where σ_i are the Pauli matrices in the spin basis, and τ_i are the Nambu matrices in the particle-hole basis. For our simulation purposes, we discretize this continuous Hamiltonian through the method of finite differences, under which $p_x = i\hbar \partial_x$ is replaced by $(\mathcal{C}_{i+1}^\dagger - \mathcal{C}_{i}^\dagger)\mathcal{C}_{i}/{a}$, where $(c_{\uparrow i}, c_{\downarrow i}, c_{\uparrow i}^\dagger, c_{\downarrow i}^\dagger) = \mathcal{C}_{i} = \Psi(x = ia)$ is the discretized annihilation operator and a is the lattice constant. Under this transformation, the Hamiltonian becomes:

$$\begin{equation} \begin{aligned} \mathcal{H} = \sum_i \mathcal{C}_i^\dagger \left[\left(2t -\mu + V_i\right) \tau_z \sigma_0 + V_z\tau_z \sigma_x+\Delta \tau_y \sigma_y \right]\mathcal{C}_i + \sum_{\langle ij \rangle} \mathcal{C}_i^\dagger \left[-t\tau_z \sigma_0+i t_{SO} \tau_z \sigma_y \right]\mathcal{C}_j, \end{aligned} \end{equation} \tag{ A3 }$$

where, $t = \hbar^2/2m^*a^2$ and $t_{SO} = \alpha_{\mathrm {R}}/2a$

The entanglement entropy is a measure of entanglement in the system that is derived from the von Neumann entropy. Consider that the concerned quantum system is defined by a Hilbert Space $\mathcal{H}$, the central problem in computing any entanglement entropy of the system (whether bipartite or disconnected) can be mapped to calculating the von Neumann entropy of a partition of the system $\mathcal{H}_A \subset \mathcal{H}$. If the reduced density matrix associated with the subsystem Hilbert space $\mathcal{H}_A$ is ρ_A , then the von Neumann entropy of the subsystem A is given by

$$\begin{equation} S = -\operatorname{Tr}\left[\rho_{A} \text{ln} \rho_{A}\right] = -\sum_{k} \alpha_{k} \text{ln} \alpha_{k}, \end{equation} \tag{ B1 }$$

where α_k represents the kth eigenvalue of the reduced density matrix ρ_A . The steps to calculate entanglement entropy involve first, the calculation of the ground state $\left|\Phi_{0}\right\rangle$ of the full system $\mathcal{H}$, the calculation of the complete density matrix $\rho = \left|\Phi_{0}\right\rangle\left\langle\Phi_{0}\right|$ of the full system, taking partial trace of the full density matrix to obtain the reduced density matrix $\rho_{A} = \operatorname{Tr}_{B}(\rho)$ and then finally diagonalizing ρ_A to obtain α_k .

The aforementioned procedure is quite computationally tedious. Therefore, we employ the free fermion method of calculating entanglement entropy. If the Hamiltonian of the system is free fermionic in nature, i.e. it contains no bi-quadratic interactions, then we can use Wick's theorem to express the reduced density matrix as

$$\begin{equation} \rho_{A} = K \text{exp} \left(-\hat{H}_{A}\right), \end{equation} \tag{ B2 }$$

where $K = |\alpha|^{2}$ and $\hat{H}_{A}$ is the effective Hamiltonian for the sub-system A, often known as the entanglement Hamiltonian of the sub-system, which has the same general form as the original Hamiltonian $\hat{H}$, i.e. if the complete system Hamiltonian has N sites and is given in the Nambu space by:

$$\begin{equation} \hat{H} = \sum_{i, j = 1}^{4 N} {\mathcal{C}}_{i}^{\dagger} H_{i j} {\mathcal{C}}_{j} \end{equation} \tag{ B3 }$$

where ${\mathcal{C}^{\dagger}} = \left(c_{1\uparrow}^{\dagger}, c_{1\downarrow}^{\dagger}, \ldots, c_{N\downarrow}^{\dagger}, c_{1\uparrow}, \ldots, c_{N\uparrow}, c_{N\downarrow}\right)$ . then the subsystem A with M sites has an effective Hamiltonian $\hat{H}_{A}$ of a similar form

$$\begin{equation} \hat{H}_{A} = \sum_{i, j = 1}^{4 M} \tilde{\mathcal{C}}_{i}^{\dagger} \tilde{H}_{i j} \tilde{\mathcal{C}}_{j} \end{equation} \tag{ B4 }$$

where ${\tilde{\mathcal{C}}^{\dagger}} = \left(c_{1\uparrow}^{\dagger}, c_{1\downarrow}^{\dagger}, \ldots, c_{M\downarrow}^{\dagger}, c_{1\uparrow}, \ldots, c_{M\uparrow}, c_{M\downarrow}\right)$.

$\hat{H}_{A}$ can now be diagonalized by defining $\tilde{\Psi}^{\dagger} = \tilde{\mathcal{C}^{\dagger}} \tilde{S} = \left(a_{2M}, \ldots, a_{1}, a_{1}^{\dagger}, \ldots, a_{2M}^{\dagger}\right)$ with $\tilde{S}$ such that $\tilde{S}^{\dagger} \tilde{H}\tilde{S} = \tilde{D} = \operatorname{diag}\left(-\varepsilon_{2M}, \ldots,-\varepsilon_{1}, \varepsilon_{1}, \ldots, \varepsilon_{2M}\right)$. Here, ɛ are the eigenvalues and $\tilde{S}$ is the eigenvector matrix of $\tilde{H}$. In this new basis the diagonalized Hamiltonian can be expressed as:

$$\begin{equation} \begin{aligned} \hat{H}_{A} = \sum_{i j = 1}^{2 M} \tilde{\mathcal{C}}_{i}^{\dagger} \tilde{H}_{i j} \tilde{\mathcal{C}}_{j} = \sum_{j = 1}^{2 M} \Psi_{j}^{\dagger} \tilde{D}_{j j} \Psi_{j} = \sum_{j = 1}^{2M} a_{j}^{\dagger} \varepsilon_{j} a_{j}-a_{j} \varepsilon_{j} a_{j}^{\dagger} = 2 \sum_{j = 1}^{2M} \varepsilon_{j} a_{j}^{\dagger} a_{j}-E_{0} \end{aligned} \end{equation} \tag{ B5 }$$

where $E_0 = \sum_{j = 1}^{2M} \varepsilon_{j}$.Once the eigenvalues of $\tilde{H}$ and therefore $\hat{H}_A$ are known, it is straightforward to calculate the eigenvalues of ρ_A from (B2) and get the entanglement entropy from (B1).

However, the issue is that the matrix elements of $\tilde{H}$ are not known. This is where we leverage the idea of correlation matrix. As we will see, the eigenvalues of the truncated correlation matrix are related to the eigenvalues of subsystem Hamiltonian matrix $\tilde{H}$ and hence $\hat{H}_A$. But first, we define the correlation matrix and the truncated correlation matrix and how they are computed.

The correlation matrix for the entire system with N sites is defined as:

$$\begin{equation} C_{n m} = \left(\begin{array}{cc} \left\langle c^{\dagger} c\right\rangle & \left\langle c^{\dagger} c^{\dagger}\right\rangle \\ \langle c c\rangle & \left\langle c c^{\dagger}\right\rangle \end{array}\right)_{n m} = \left\langle\mathcal{C}_{n}^{\dagger} \mathcal{C}_{m}\right\rangle = \operatorname{Tr}\left[\rho \mathcal{C}_{n}^{\dagger} \mathcal{C}_{m}\right]. \end{equation} \tag{ B6 }$$

The expectation values are calculated in the ground-state of the Hamiltonian. The elements of the correlation matrix have the form $\left\langle\mathcal{C}_{n}^{\dagger} \mathcal{C}_{m}\right\rangle = \left\langle c_{i\sigma}^{\dagger} c_{j\sigma^{^{\prime}}}\right\rangle = \left\langle\Phi_{0}\left|c_{i\sigma}^{\dagger} c_{j\sigma^{^{\prime}}}\right| \Phi_{0}\right\rangle$, $i, j$ represents sites and $\sigma, \sigma^{^{\prime}}$ represent spins, $|\Phi_0 \rangle$ is the ground state of the Hamiltonian.

To evaluate these expectation values, we must first calculate the ground-state of the complete system Hamiltonian. Following from our previous technique, we can diagonalize this Hamiltonian matrix with $\Psi^{\dagger} = \mathcal{C}^{\dagger}S = \left(a_{2N}, \ldots, a_{1}, a_{1}^{\dagger}, \ldots, a_{2N}^{\dagger}\right)$, such that $S^{\dagger}H S = D$. The ground state is then constructed with the new BdG operators a as

$$\begin{equation} | \Phi_0 \rangle = \prod_{i}^{2N} a_{i} |0\rangle \end{equation} \tag{ B7 }$$

where $|0\rangle$ is the fermionic vaccum state.

We then use $|\Phi_0\rangle$ to compute the expectation values

$$\begin{equation} \begin{aligned} \left\langle\Phi_{0}\left|{c}^{\dagger} {c}\right| \Phi_{0}\right\rangle_{nm} & = \sum_{k l}\left\langle\Phi_{0}\left|\Psi_{k}^{\dagger} S_{k n}^{T} S_{m l} \Psi_{l}\right| \Phi_{0}\right\rangle \\ & = \sum_{k l \unicode{x2A7D} 2N}\left\langle\Phi_{0}\left|a_{2N+1-k} S_{k n}^{T} S_{m l} a_{2N+1-l}^{+}\right| \Phi_{0}\right\rangle \\ & = \sum_{k = 1}^{2N} S_{k, n}^{T} S_{m, k}. \end{aligned} \end{equation} \tag{ B8 }$$

Similarly we can get the other terms:

$$\begin{equation} \begin{aligned} \langle\Phi_{0}\left|c c\right| &\Phi_{0} \rangle_{nm} \\ & = \sum_{k>2N, l \leq 2N} S_{n, k} S_{m, l}\left\langle\Phi_{0}\left|a_{k-2N} a_{2N+1-l}^{\dagger}\right| \Phi_{0}\right\rangle \\ & = \sum_{k = 1}^{2N} S_{4N+1-k, n}^{T} S_{m,k}, \end{aligned} \end{equation} \tag{ B9 }$$

and

$$\begin{equation} \begin{aligned} \left\langle\Phi_{0}\left|c^{\dagger} c^{\dagger}\right| \Phi_{0}\right\rangle_{nm} = \sum_{k = 1}^{2N} S_{k, n}^{T} S_{m,4N+1-k} \\ \left\langle\Phi_{0}\left|c c^{\dagger}\right| \Phi_{0}\right\rangle_{nm} = \sum_{k = 1}^{2N} S_{k+2N, n}^{T} S_{m, k+2N}. \end{aligned} \end{equation} \tag{ B10 }$$

The above relations give us the entire correlation matrix C. The truncated correlation matrix $\tilde{C}$ for a sub-system A with M sites is a sub-matrix of the net correlation matrix C which is defined as follows:

$$\begin{equation} \tilde{C}_{n m} = \left(\begin{array}{cc} \left\langle c^{\dagger} c\right\rangle & \left\langle c^{\dagger} c^{\dagger}\right\rangle \\ \langle c c\rangle & \left\langle c c^{\dagger}\right\rangle \end{array}\right)_{n, m \in A} = \operatorname{Tr}\left[\rho_{A} \tilde{\mathcal{C}}_{n}^{\dagger} \tilde{\mathcal{C}}_{m}\right]. \end{equation} \tag{ B11 }$$

Now that we have defined the truncated correlation matrix, we return to our original task of deriving a relation between the eigenvalues of the entanglement Hamiltonian matrix $\hat{H}_A$ and the eigenvalues of the truncated correlation matrix. We start by substituting (B2) in (B11) to get:

$$\begin{equation} \tilde{C}_{n m} = \operatorname{Tr}\left[K \text{exp} \left(-\hat{H}_{A}\right) \tilde{\mathcal{C}}_{n}^{\dagger} \tilde{\mathcal{C}}_{m}\right]. \end{equation} \tag{ B12 }$$

Substituting (B5) in (B12) gives:

$$\begin{equation} \tilde{C} = K \operatorname{Tr}\left[\text{exp} \left(-2 \sum_{j} \varepsilon_{j} a_{j}^{\dagger} a_{j}+E_{0}\right) \sum_{k l} \tilde{\Psi}_{k}^{\dagger} \tilde{S}_{k n}^{T} \tilde{S}_{m l} \tilde{\Psi}_{l}\right]. \end{equation} \tag{ B13 }$$

This can be simplified as:

$$\begin{equation} \begin{aligned} \tilde{C} & = K e^{E_{0}} \sum_{k = 1}^{4 M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \operatorname{Tr}\left[\text{exp} \left(-2 \sum_{j} \varepsilon_{j} n_{j}\right) \tilde{\Psi}_{k}^{\dagger} \tilde{\Psi}_{k}\right] \\ & = K e^{E_{0}} \sum_{k \unicode{x2A7D} 2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \operatorname{Tr}\left[\text{exp} \left(-2 \sum_{j} \varepsilon_{j} n_{j}\right)\left(1-n_{k}\right)\right] \\ &\quad + K e^{E_{0}} \sum_{k>2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \operatorname{Tr}\left[\text{exp} \left(-2 \sum_{j} \varepsilon_{j} n_{j}\right) n_{k}\right], \end{aligned} \end{equation} \tag{ B14 }$$

where $n_{k} = a^{\dagger}_{k}a_{k}$. The value of K can be evaluated by imposing the condition $Tr[\rho_A] = 1$, which is a fundamental property of density matrices. We can now uncover the following:

$$\begin{equation} \begin{aligned} \operatorname{Tr}\left[\rho_{A}\right] & = K e^{E_{0}} \operatorname{Tr}\left[\text{exp} \left(-2 \sum_{j} \varepsilon_{j} n_{j}\right)\right] = 1 \\ & = K e^{E_{0}}\left(\prod_{l = 1}^{2M} \sum_{n_{l} = 0,1}\right)\left\langle\Phi_{0}\left|\prod_{l = 1}^{2M} a_{l}^{n_{l}} e^{-2 \varepsilon_{l} n_{l}} a_{l}^{\dagger n_{l}}\right| \Phi_{0}\right\rangle \\ & = K e^{E_{0}} \prod_{l = 1}^{2M} \sum_{n_{l} = 0,1} e^{-2 \varepsilon_{l} n_{l}} = K e^{E_{0}} \prod_{l = 1}^{2M}\left(1+e^{-2 \varepsilon_{l}}\right) \\ &\Longrightarrow K = e^{-E_{0}} \prod_{l = 1}^{2M} \frac{1}{1+e^{-2 \varepsilon_{l}}}. \end{aligned} \end{equation} \tag{ B15 }$$

Using the above two relations we get:

$$\begin{equation} \begin{aligned} \tilde{C}_{n m} & = \sum_{k \leq 2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k}\left(1-\frac{e^{-2 \varepsilon_{M+1-k}}}{1+e^{-2 \varepsilon_{M+1-k}}}\right) +\sum_{k>2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \frac{e^{-2 \varepsilon_{k-M}}}{1+e^{-2 \varepsilon_{k-M}}} \\ & = \sum_{k \leq 2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \frac{1}{1+e^{-2 \varepsilon_{M+1-k}}}+\sum_{k>2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \frac{1}{1+e^{2 \varepsilon_{k-M}}}. \end{aligned} \end{equation} \tag{ B16 }$$

We can expand $\tilde{H}$ in terms of its eigenvector matrix $\tilde{S}$ and eigenvalues ɛ_k as

$$\begin{equation} \begin{aligned} H_{n m}& = \sum_{k} \tilde{S}_{n k} \lambda_{k} \tilde{S}_{k m}^{T} \\ & = \sum_{k \leq 2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k}\left(-\varepsilon_{2M+1-k}\right) +\sum_{k>2M} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \varepsilon_{k-2M}, \end{aligned} \end{equation} \tag{ B17 }$$

where $\lambda_{k} = D_{kk}$. Thus in terms of λ_k ,

$$\begin{equation} \begin{aligned} H_{n m} & = \sum_{k} \lambda_{k} \tilde{S}_{k n}^{T} \tilde{S}_{m k} \\ \tilde{C}_{n m} & = \sum_{k} \frac{1}{1+e^{2 \lambda_{k}}} \tilde{S}_{k n}^{T} \tilde{S}_{m k} = \sum_{k} \xi_{k} \tilde{S}_{k n}^{T} \tilde{S}_{m k}. \end{aligned} \end{equation} \tag{ B18 }$$

It is clear that $\tilde{H}$ and $\tilde{C}$ have the same set of eigenvectors and their eigenvalues λ_k and ɛ_k are related as:

$$\begin{equation} \begin{array}{l} \xi_{k} = \frac{1}{1+e^{2 \lambda_{k}}}, \mspace{50mu} \lambda_{k} = \frac{1}{2} \text{log} \left(\frac{1-\xi_{k}}{\xi_{k}}\right). \end{array} \end{equation} \tag{ B19 }$$

These relations finally lead us to the equation for entanglement entropy in terms of the eigenvalues of the truncated correlation matrix ɛ_k :

$$\begin{equation} \begin{aligned} S_{\mathrm {ent}} & = -\operatorname{Tr}\left[\rho_{A} \text{log} \rho_{A}\right]\\ & = \sum_{k = 1}^{2M}\left[\frac{2 \varepsilon_{k}}{1+e^{2 \varepsilon_{k}}}+\text{log} \left(1+e^{-2 \varepsilon_{k}}\right)\right] \\ & = -\sum_{k = 1}^{2M}\left[\xi_{k} \text{log} \xi_{k}+\left(1-\xi_{k}\right) \text{log} \left(1-\xi_{k}\right)\right]. \end{aligned} \end{equation} \tag{ B20 }$$

C.1. Current calculations

First, we start with the assumption that the leads are large enough so that they can be considered semi-infinite, and thus be effectively modeled by infinite dimensional matrices. This assumption allows using the NEGF approach [51] by partitioning the channel and the leads and working with the channel Green's function and incorporating the leads via self energies figure 8. The retarded channel Green's function can now be written as:

$$\begin{equation} G^{r}\left(E\right) = \left[\left(E+i \eta\right) I-H_{\mathrm {BdG}}-\Sigma_{L}^{r}-\Sigma_{R}^{r}\right]^{-1}, \end{equation} \tag{ C1 }$$

where E is the free variable energy, and I is the identity matrix of the dimension of the Hamiltonian, η is a small positive damping parameter, and $\Sigma^{r}_L$ and $\Sigma^{r}_R$ represent the retarded self energies for the semi-infinite contacts.

Zoom In Zoom Out Reset image size

The self energies can be represented in two possible choice of basis: (i) in the eigenbasis of the real-space Hamiltonian, (ii) the eigen basis with wide-band approximation, where the self energy is encapsulated via a parameter $\Gamma_{L/R}$. The contact Hamiltonian can be written in its own eigenbasis as:

$$\begin{equation} \hat{H}_{C} = \sum_{k \sigma \alpha \in L, R} \epsilon_{k \alpha} c_{k \sigma \alpha}^{\dagger} c_{k \sigma \alpha} \end{equation} \tag{ C2 }$$

where $c^{(\dagger)}_{k \sigma \alpha}$ represents the creation/annihilation operator of the contact momentum eigenstate of spin $\sigma = \uparrow\downarrow$ with eigenvalue $\epsilon_{k \alpha}$.

The coupling between the contacts and the nanowire is given by the tunneling Hamiltonian:

$$\begin{equation} \hat{H}_{T} = \sum_{k, m, \sigma, \alpha \in L, R}\left[\tau_{k, \alpha, m} c_{k \sigma \alpha}^{\dagger} d_{m, \sigma}+\tau_{k \alpha m}^{*} d_{m, \sigma}^{\dagger} c_{k \sigma \alpha}\right] \end{equation} \tag{ C3 }$$

where $d^{(\dagger)}_{m \sigma}$ represents the annihilation/creation operator of the m_th site of the Kitaev chain and $\tau_{k \alpha, m}$ is an element of the tunnel coupling matrix between the $\alpha = (L/R)$ contact's momentum state k and the m^th site of the Kitaev chain.

In the basis of the real-space Hamiltonian, the self-energy matrices can be evaluated by a straightforward calculation of the matrix elements $\Sigma^{r}_{\alpha}\left(i,i\right) = \sigma_{\alpha \in L/R}$ :

$$\begin{equation} \begin{aligned} \sigma_{\alpha \in L / R}\left(E\right)& = \sum_{k} \frac{\left|\tau_{\alpha, k}\right|^{2}}{\left(E-\epsilon_{k \alpha}+i \eta\right)} \\ & = P \int d \epsilon_{k \alpha} \frac{\left|\tau_{\alpha, k}\right|^{2}}{\left(E-\epsilon_{k \alpha}\right)}-i \pi \int d \epsilon_{k \alpha}\left|\tau_{\alpha, k}\right|^{2} \delta\left(E-\epsilon_{k \alpha}\right), \end{aligned} \end{equation} \tag{ C4 }$$

where P stands for the Cauchy principal value. Recasting this integral into a real and an imaginary part ($\sigma_{\alpha} = \epsilon_{\alpha}\left(E\right)-i\gamma_{\alpha}/2$) provides a better physical picture , with the imaginary part $\gamma_{\alpha} = i\left(\sigma^r_{\alpha}-\sigma^a_{\alpha}\right)$ representing the level broadening, and the real part _α calculated from the principal value integral.

In this work, we stick to the wide-band approximation basis. Under this approximation, $\tau_{k\alpha,m}$ is independent of energy and thus the principal value integral representing the real-part of the self energy vanishes and we are just left with the broadening matrix $\Gamma_{\alpha}\left(E\right)$ such that the self energies can be represented as $\Sigma^r_{\alpha} = -i\Gamma_{\alpha}/2$, which can be treated as an input parameter. Once we are done calculating the self energies, we need to evaluate the 'lesser' Green's function [51, 52]:

$$\begin{equation} G^{<}\left(E\right) = G^{r}\left(E\right)\left(\Sigma_{L}^{<}\left(E\right)+\Sigma_{R}^{<}\left(E\right)\right) G^{a}\left(E\right), \end{equation} \tag{ C5 }$$

where G^a is the advanced Green's function (Hermitian conjugate of the retarded Green's function calculated from (C1)). Calculating the lesser Green's function is equivalent to calculating the electron (hole) correlation function. Here, $\Sigma^{\lt}_{L,R}(E)$ represent the in-scattering functions from leads L and R respectively which can be evaluated using:

$$\begin{equation} \begin{aligned} \Sigma_{\alpha}^{<}\left(E\right) = -\left[\Sigma_{\alpha}^{r}\left(E\right)-\Sigma_{\alpha}^{a}\left(E\right)\right] f_{\alpha}\left(E\right) = i \Gamma_{\alpha}\left(E\right) f_{\alpha}\left(E\right), \end{aligned} \end{equation} \tag{ C6 }$$

where $f_{\alpha} = f(E-\mu_{\alpha})$ represents the Fermi–Dirac distribution in either leads $\alpha = L(R)$ with a chemical potential µ_α . In this formulation, all different components of the currents can then be deduced from the current operator formula. The elements of the current operator matrix can be deduced using the correlation function and the hopping elements of the Hamiltonian. A rigorous but straightforward calculation based on fundamental consideration gives the current operator matrix:

$$\begin{equation} \begin{aligned} I_{L}^{\mathrm {op}}\left(E\right) = \frac{e}{h} \left[G^{r}\left(E\right) \Sigma_{L}^{<}-\Sigma_{L}^{<} G^{a}\left(E\right) +G^{<}\left(E\right) \Sigma_{L}^{a}\left(E\right)-\Sigma_{L}^{r}\left(E\right) G^{<}\left(E\right)\right]. \end{aligned} \end{equation} \tag{ C7 }$$

Taking the trace of the above equation gives the net charge current, which is given by the current formula:

$$\begin{equation} \begin{aligned} I_{L}\left(E\right) = \frac{ie}{h} \textrm{Trace}\left[\Gamma_{L}\left(E\right) f_{L}\left(E\right)\left(G^{r}\left(E\right)-G^{a}\left(E\right)\right)+ \Gamma_{L}\left(E\right) G^{<}\left(E\right)\right]. \end{aligned} \end{equation} \tag{ C8 }$$

The above formula can be further simplified by the notation introduced in with the spectral function $A(E) = i(G^r-G^a)$ and the electron correlation matrix $G^n(E) = -iG^{\lt}(E)$ as

$$\begin{equation} I_{L}\left(E\right) = \frac{e}{h}\textrm{Trace }\left[\Gamma_{L}\left(E\right) f_{L}\left(E\right) A\left(E\right)-\Gamma_{L}\left(E\right) G^{n}\left(E\right)\right]. \end{equation} \tag{ C9 }$$

Using simple manipulations as given below we get:

$$\begin{equation} \begin{aligned} A\left(E\right) & = i\left(G^{r}-G^{a}\right) = iG^{r} \left[\left(G^{a}\right)^{-1}-\left(G^{r}\right)^{-1}\right] G^{a} \\ & = iG^{r} \left[\Sigma^{r}-\Sigma^{a}\right] G^{a} = G^{r} \Gamma G^{a}, \end{aligned} \end{equation} \tag{ C10 }$$

where $\Gamma(E) = \Gamma_L(E)+\Gamma_R(E)$. Using the above in conjunction with (C5) and (C6), along with the properties of the trace operation, yields the Landauer transmission formula for the current as

$$\begin{equation} I_{L} = \int {\mathrm {d}} E \textrm{Trace}\left[\Gamma_{L} G^{r} \Gamma_{R} G^{a}\right]\left(f_{L}\left(E\right)-f_{R}\left(E\right)\right), \end{equation} \tag{ C11 }$$

where $f_{L(R)}(E)$ represents the Fermi–Dirac distribution in either lead.

However, in the case involving superconductors, either functioning as a contact or as the central system or both, the net current through the contact α is the difference between electron and hole currents in Nambu space. We can now use the full fledged current operator in (C7) to evaluate the current from first principles as

$$\begin{equation} \begin{aligned} I_{L}^{\mathrm {op}}\left(E\right) = \frac{e}{h} \tau_{z}\left[G^{r}\left(E\right) \Sigma_{L}^{<}\left(E\right) -\Sigma_{L}^{<}\left(E\right) G^{a}\left(E\right) +G^{<}\left(E\right) \Sigma_{L}^{a}\left(E\right)-\Sigma_{L}^{r}\left(E\right) G^{<}\left(E\right)\right], \end{aligned} \end{equation} \tag{ C12 }$$

where $\tau_z = \sigma_z \otimes I_{N \times N}$ , where is the Pauli-z matrix and is the identity matrix of dimension N × N. The next step is to take a trace of the above equation, which in general need not take the form of (C8), typically when the contacts are superconducting due to the non-diagonal structure of $\Sigma_{L(R)}$, and may lead to erroneous results, specifically when evaluating the Josephson currents. However, in our case, as the contacts are normal, and hence diagonal, (C12) indeed takes the form of (C8) with the net current becoming

$$\begin{equation} I_{\alpha} = \int {\mathrm {d}}E \frac{\left(I_{\alpha}^{\left(e\right)}\left(E\right) -I_{\alpha}^{\left(h\right)}\left(E\right)\right)}{2}, \end{equation} \tag{ C13 }$$

where each current $I^{e(h)}_{\alpha}$ can be evaluated separately using the form in (C8) or (C9). The factor of a half comes since the original Hamiltonian been doubled while writing the BdG Hamiltonian in order to write it consistently in the Nambu space. It is then instructive in our context to note that various matrices defined within the approach will have a matrix structure due to the electron–hole Nambu space. In particular, the contact broadening matrices etc can be written with a general diagonal structure (due to the diagonal structure of normal contacts in Nambu space) as $\Gamma_{\alpha} = \Gamma^{ee}_{\alpha}+\Gamma^{hh}_{\alpha}$, where the superscripts ee(hh) represent the electron (hole) diagonal part of the self energy or broadening matrix (remember that the electron (hole) diagonal part of any matrix only separates sectors in the Nambu space and not the spin space, and therefore both the parts contain both spin types). Following this and using the above observations on the current operator formula in (C9), we obtain the electron (hole) current across the TS as a sum of three components [47, 53]:

$$\begin{equation} \begin{aligned} I_{L}^{e\left(h\right)}\left(E\right) & = \frac{e}{h}\left( T_D\left[f_{L}^{e e\left(h h\right)}\left(E\right)-f_{R}^{e e\left(h h\right)}\left(E\right)\right]\right) \longrightarrow \left(i\right)\\ &\quad +\frac{e}{h}\left(T_A\left[f_{L}^{e e\left(h h\right)}\left(E\right)-f_{L}^{h h\left(e e\right)}\left(E\right)\right]\right) \longrightarrow \left(ii\right)\\ &\quad +\frac{e}{h}\left(T_{CAR}\left[f_{L}^{e e\left(h h\right)}\left(E\right)-f_{R}^{h h\left(e e\right)}\left(E\right)\right]\right) \longrightarrow \left(iii\right) \end{aligned} \end{equation} \tag{ C14 }$$

and

$$\begin{equation} \begin{aligned} T_D &= \textrm{Trace }\left(\Gamma_{L}^{e e\left(h h\right)} G^{r} \Gamma_{R}^{e e\left(h h\right)} G^{a}\right) \\ T_A &= \textrm{Trace }\left(\Gamma_{L}^{e e\left(h h\right)} G^{r} \Gamma_{L}^{h h\left(e e\right)} G^{a}\right) \\ T_{CAR} &= \textrm{Trace }\left(\Gamma_{L}^{e e\left(h h\right)} G^{r} \Gamma_{R}^{h h\left(e e\right)} G^{a}\right) \end{aligned} \end{equation} \tag{ C15 }$$

where the term (i) represents the direct transmission process of either the electron or the hole, (ii) represents the direct Andreev transmission and (iii) represents the crossed Andreev transmission. At this point it is worth noting that $f^{ee}_{\alpha} = f(E-\mu_{\alpha})$, $f^{hh}_{\alpha} = f(E+\mu_{\alpha})$. The matrix structure of the broadening matrix Γ itself is diagonal such that $\Gamma^{ee}_L = \Gamma_L(1,1) = \gamma$, $\Gamma^{hh}_L = \Gamma_L(2,2) = \gamma$, $\Gamma^{ee}_R = \Gamma_R(2N-1,2N-1) = \gamma$, $\Gamma^{hh}_R = \Gamma_R(2N,2N) = \gamma$ and zeros otherwise. The Andreev transmission can then be written as $T_A(E) = \gamma^2 | G^{r,eh}_{11}(E) |^2$ and the crossed Andreev transmission as $T_{CA}(E) = \gamma^2 | G^{r,eh}_{1N}(E) |^2$ the direct transmission can be written as $T_D(E) = \gamma^2 | G^{r,ee}_{1N}(E) |^2$. Here $G^{r,ee(hh)}_{ij}$ represents the i, j element of the electron (hole) diagonal block of the retarded Green's function in Nambu space and $G^{r,eh(he)}_{ij}$ represents that of the off-diagonal block of the retarded Green's function in Nambu space.

C.2. Local and non-local conductances

We start with the current formula for a generic bias situation, (V_L ,V_R ) for the left and right contacts with $\mu_{L} = eV_{L}$ and $\mu_{R} = eV_{R}$

The current through one of the contacts has contributions from both electron and hole flow. The net current is given by:

$$\begin{equation} I^{L} = \frac{I_{L}^{\left(e\right)}-I_{L}^{\left(h\right)}}{2}. \end{equation} \tag{ C16 }$$

The individual electron and hole components can then be derived using:

$$\begin{equation} \begin{aligned} I_{L}^{\left(e\right)} &= -\frac{e}{h} \left\{\int {\mathrm {d}} E T_{A}^{\left(e\right)}\left(E\right)\left[f\left(E-e V_{L}\right)-f\left(E+e V_{L}\right)\right]\right.\\ &\quad +\int {\mathrm {d}} E T_{C A R}^{\left(e\right)}\left(E\right)\left[f\left(E-e V_{L}\right)-f\left(E+e V_{R}\right)\right] \\ &\quad \left.+\int {\mathrm {d}} E T_{D}^{\left(e\right)}\left(E\right)\left[f\left(E-e V_{L}\right)-f\left(E-e V_{R}\right)\right]\right\}, \end{aligned} \end{equation} \tag{ C17 }$$

$$\begin{equation} \begin{aligned} I_{L}^{\left(h\right)}& = -\frac{e}{h} \left\{\int {\mathrm {d}} E T_{A}^{\left(h\right)}\left(E\right)\left[f\left(E+e V_{L}\right)-f\left(E-e V_{L}\right)\right]\right.\\ &\quad +\int {\mathrm {d}} E T_{C A R}^{\left(h\right)}\left(E\right)\left[f\left(E+e V_{L}\right)-f\left(E-e V_{R}\right)\right] \\ &\quad \left.+\int {\mathrm {d}} E T_{D}^{\left(h\right)}\left(E\right)\left[f\left(E+e V_{L}\right)-f\left(E+e V_{R}\right)\right]\right\}. \end{aligned} \end{equation} \tag{ C18 }$$

Here, $I^{(e)}_L$ and $I^{(h)}_L$ can be shown to be equal and opposite using symmetry conditions and especially:

$$\begin{equation} \begin{aligned} T_{CAR}^{\left(e\right)}\left(E\right) \equiv T^{\left(h\right)}_{C A R}\left(E\right), \mspace{50mu} T_{D}^{\left(e\right)}\left(E\right) \equiv T_{D}^{\left(h\right)}\left(-E\right). \end{aligned} \end{equation} \tag{ C19 }$$

It is thus sufficient to consider either one of $I^{(e)}_L$ and $I^{(h)}_L$ while calculating the net current through a contact. Let us now consider the conductance matrix:

$$\begin{equation} \mathrm{G} = \left(\begin{array}{cc} G_{L L} & G_{L R} \\ G_{R L} & G_{R R} \end{array}\right) = \left(\begin{array}{cc} \left.\frac{\partial I_{L}}{\partial V_{L}}\right|_{V_{R} = 0} & \left.\frac{\partial I_{L}}{\partial V_{R}}\right|_{V_{L} = 0} \\ \left.\frac{\partial I_{R}}{\partial V_{L}}\right|_{V_{R} = 0} & \left.\frac{\partial I_{R}}{\partial V_{R}}\right|_{V_{L} = 0}. \end{array}\right) \end{equation} \tag{ C20 }$$

Now let us consider only the electron current for now since $I_L = I^{(e)}_L$. Using this fact we can evaluate the Local conductance as:

$$\begin{equation} \begin{aligned} G_{L L} &\equiv \left. \frac{\partial I_{L}}{\partial V_{L}} \right|_{V_{R}\equiv 0 }\\ & \equiv \frac {e } {h } \left\{\frac {\partial } {\partial V_L } \left[\int {\mathrm {d}} E T_{A}\left(E\right)\left[f\left(E-e V_L\right)-f\left(E+eV_L\right)\right]\right.\right.\\ \mspace{90mu}&\quad +\int {\mathrm {d}} E T_{CAR}\left(E\right)\left[f\left(E-e V_L\right)-f\left(E\right)\right] \\ \mspace{90mu}&\quad +\left.\left.\int {\mathrm {d}} E T_{D}\left(E\right)\left[f\left(E-eV_L\right)-f\left(E\right)\right]\right]\right\} \\ & \equiv \frac{e}{h}\left[\int {\mathrm {d}} E\right. T_{A}\left(E\right)\left[\frac{\partial f\left(E-eV_L\right)}{\partial V_L}-\frac{\partial f\left(E+e V_L\right)}{\partial V_L}\right] \\ \mspace{-50mu}&\quad + \int {\mathrm {d}} E T_{CAR}\left(E\right)\left[\frac{\partial f\left(E-e V_L\right)}{\partial V_L}\right] +\left.\int {\mathrm {d}} E T_{D}\left(E\right)\left[\frac{\partial f\left(E-e V_L\right)}{\partial V_L}\right]\right]. \end{aligned} \end{equation} \tag{ C21 }$$

At zero temperature, the $f\left(x\right) = \Theta\left(-x\right)$ where Θ is the Heaviside step function. This implies $f\left(E-eV \right) = \Theta\left(eV-E\right)$ and similarly we can write the other terms. Since the derivative of the heaviside function is the Dirac Delta function, we get:

$$\begin{equation} \begin{aligned} \left.G_{LL}\left(V\right)\right|_{T \rightarrow 0} &\equiv \frac{e^2}{h}\left[T_{A}\left(E = e V\right)+T_{A}\left(E = -e V\right)\right. + \left.T_{C A R}\left(E = eV\right)+T_{D}\left(E = e V\right)\right]. \end{aligned} \end{equation} \tag{ C22 }$$

Now we can similarly get the non local conductance.

$$\begin{equation} \begin{aligned} G_{L R} = & \left.\frac{\partial I_L}{\partial V_R}\right|_{V_L \equiv 0} \\ = & \frac{e}{h}\left\{\frac {\partial } {\partial V _ {R } } \left[\int {\mathrm {d}} E T_{C A R}\left(E\right)\left[f\left(E\right)-f\left(E+e V_R\right)\right]\right.\right. \\ & +\int {\mathrm {d}} E T_D\left(E\right)\left[f\left(E\right)-f\left(E-e V_R\right)\right] \\ = & \frac{e}{h}\left[\int {\mathrm {d}} E T_{C A R}\left(E\right)\left[\frac{-\partial f\left(E+e V_R\right)}{\partial V_R}\right]\right. \\ & \left.+\int {\mathrm {d}} E T_D\left(E\right)\left[\frac{-\partial f\left(E-e V_R\right)}{\partial V_R}\right]\right]. \end{aligned} \end{equation} \tag{ C23 }$$

Once again, using the property of Fermi–Dirac functions as $T\rightarrow 0$ as $f(x) = \Theta(-x)$ we get:

$$\begin{equation} \begin{aligned} \left.G_{LR}\left(V\right)\right|_{T \rightarrow 0} \equiv \frac{e^2}{h}\left[T_{D}\left(E = e V\right)-T_{CAR}\left(E = -e V\right)\right]. \end{aligned} \end{equation} \tag{ C24 }$$

In order to explain the entanglement entropy signatures for the Rashba nanowire from a much more elementary and first-principles perspective, we highlight the parallels between the Rashba nanowire and a toy model consisting two distinct p-wave superconductors interacting through an s-wave pairing potential. Starting with the normal part of the Rashba nanowire Hamiltonian, without proximity-induced superconductivity, in the continuous momentum space:

$$\begin{equation} \mathcal{H}_0 = \int {\mathrm {d}}k \Psi^{\dagger}\left(k\right) \left[\frac{-\hbar^2 k^2}{2 m}-\mu+i\alpha_{\mathrm {R}} k \sigma_y +V_z \sigma_x\right] \Psi\left(k\right), \end{equation} \tag{ D1 }$$

where, $\alpha_{\mathrm {R}}$ is the Rashba spin–orbit coupling strength, $V_z = \frac{g\mu_{\mathrm {B}}B}{2}$ is the Zeeman field, k is the momentum and $\Psi(k) = (\psi_\uparrow(k), \: \psi_\downarrow(k))^T$ is the fermionic spinor at k. The eigenvalues $\varepsilon_\pm$ of this Hamiltonian are given by $\varepsilon_{k, \pm} = \frac{\hbar^2 k^2}{2 m}-\mu \pm \sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2}$. The eigen-spinors in terms of the original spinor are given by:

$$\begin{equation} \left(\begin{array}{c} \psi_+\left(k\right)\\ \psi_-\left(k\right) \end{array}\right) = \left(\begin{array}{cc} \frac{\left(-i \alpha_{\mathrm {R}} k+V_z\right)}{\sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2}} & 1 \\[8pt] \frac{\left(i \alpha_{\mathrm {R}} k-V_z\right)}{\sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2}} & 1 \end{array}\right) \left(\begin{array}{c} \psi_\uparrow\left(k\right)\\[5pt] \psi_\downarrow\left(k\right) \end{array}\right). \end{equation} \tag{ D2 }$$

The spinor basis, defined by $\Phi^{^{^{\prime}}} = (\psi_+(k), \: \psi_-(k))^T$ is often called the helical basis. Now we express the full Hamiltonian of the nanowire, including the proximity-induced superconducting pairing terms, in the helical basis. The Hamiltonian in the helical basis [54] is given by:

$$\begin{equation} \begin{aligned} \mathcal{H} & = \frac{1}{2} \int {\mathrm {d}}k \left[\varepsilon_+\left(k\right) \psi^\dagger_+\left(k\right)\psi_+\left(k\right) + \varepsilon_-\left(k\right) \psi^\dagger_-\left(k\right)\psi_-\left(k\right)\right] \\ &\quad + \left[\frac{\Delta_{++}\left(k\right)}{2} \psi^\dagger_+\left(k\right)\psi^\dagger_+\left(-k\right) + \frac{\Delta_{-}\left(k\right)}{2} \psi^\dagger_-\left(k\right)\psi^\dagger_-\left(-k\right) \right.\\ &\quad \left. + \Delta_{+-}\psi^\dagger_+\left(k\right)\psi^\dagger_-\left(-k\right) + h.c.\right] \end{aligned} \end{equation} \tag{ D3 }$$

where,

$$\begin{equation} \begin{aligned} \varepsilon_\pm\left(k\right)& = \frac{\hbar^2 k^2}{2 m}-\mu \pm \sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2} \\ \Delta_{-,++}\left(k\right)& = \frac{\pm i \alpha_{\mathrm {R}} k \Delta}{\sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2}}\\ \Delta_{+-}\left(k\right)& = \frac{V_z \Delta}{\sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2}}. \end{aligned} \end{equation} \tag{ D4 }$$

The Hamiltonian in (3) can be interpreted as follows: there are two bands + and − with their normal parts being described by $\varepsilon_+(k)$ and $\varepsilon_+(k)$ respectively. Within each band, there exists an intra-band superconducting pairing given by $\Delta_{-, ++}(k)$. These pairings are an odd function of k and hence manifest as intra-band p-wave pairings. There also exists an inter-band superconducting pairing between the two bands given by $\Delta_{+-}(k)$. This pairing is an even function of k and manifests as an inter-band s-wave pairing. Thus, the helical basis allows us to decompose the nanowire into two p-wave superconductors + and − interacting via an s-wave coupling. We notice two things here. First, the effective onsite potential $\mu_{\mathrm {eff}} = \mu \mp \sqrt{V_z^2+\alpha_{\mathrm {R}}^2 k^2}$ depends on a controllable parameter V_z and therefore the magnetic field B. It decreases with B for the + chain while increases with B for the − chain. Second, the p-wave pairings $\Delta_{-,++}$ are of opposite parity for both bands.

This motivates introducing a new pedagogical model, which we call the 'spinfull Kitaev chain' based on a topical p-wave superconductor—The Kitaev chain. We consider two distinct Kitaev chains, labeled with + and −, interacting via an s-wave pairing potential $\Delta_s$. The chemical potentials, hopping terms, and p-wave pairings are represented by µ_σ , t_σ , $\Delta_\sigma$, where $\sigma = + \text{or } -$. The Hamiltonian, describing this pedagogical model, in the position space (with open boundary conditions) is given by:

$$\begin{equation} \begin{aligned} \mathcal{H} & = \sum_{i = 1}^N \sum_{\sigma = \pm}\mu_\sigma c_{i, \sigma}^\dagger\: c_{i, \sigma} \\ &\quad + \sum_{i = 1}^{N-1} \sum_{\sigma = \pm} t_\sigma c_{i+1, \sigma}^\dagger\: c_{i, \sigma} + h.c.\\ &\quad + \sum_{i = 1}^{N-1} \sum_{\sigma = \pm} \Delta_\sigma c_{i+1, \sigma}^\dagger\: c_{i, \sigma}^\dagger + h.c.\\ &\quad + \sum_{i = 1}^N \Delta_s c_{i, +}^\dagger\: c_{i, -}^\dagger + h.c. , \end{aligned} \end{equation} \tag{ D5 }$$

where $c_{i,\sigma}^\dagger$ creates a fermion at site 'i' in the $\sigma = \pm$ chain. For further analysis of the model, we express our Hamiltonian in the Majorana basis defined by:

$$\begin{equation} \begin{aligned} &c_{i, \sigma} = \gamma_{i, \sigma, A}+i \gamma_{i, \sigma, B} \\ &c_{i, \sigma}^\dagger = \gamma_{i, \sigma, A}-i \gamma_{i, \sigma, B} \\ &\gamma_{i, \sigma, \tau}^\dagger = \gamma_{i, \sigma, \tau} \\ &\left\{\gamma_{i, \sigma_i, \tau_i} , \gamma_{j, \sigma_j, \tau_j} \right\} = 2\: \delta_{i,j} \: \delta_{\sigma_i,\sigma_j}\: \delta_{\tau_i, \tau_j} \end{aligned} \end{equation} \tag{ D6 }$$

where i represents the site index, $\sigma = + \text{or } -$ and $\tau = A \text{or } B$, where A and B represent two distinct types of Majoranas that constitute a single fermion. In the Majorana basis, our model comprises four interacting Majorana chains. The Hamiltonian in this basis is given by:

$$\begin{equation} \begin{aligned} H & = \sum_{i = 1}^N \sum_{\sigma = \pm}\mu_{\sigma} +\mu_{\sigma} \: \gamma_{i, \sigma, A} \: \gamma_{i, \sigma, B}\\ &\quad +\sum_{i = 1}^{N-1}\sum_{\sigma = \pm}\left(t_{\sigma}+\Delta_{\sigma}\right) \gamma_{i+1, \sigma, A} \: \gamma_{i, \sigma, B}\\ &\quad +\sum_{i = 1}^{N-1}\sum_{\sigma = \pm}\left(t_{\sigma}-\Delta_{\sigma}\right) \gamma_{i, \sigma, A} \: \gamma_{i+1, \sigma, B} \\ &\quad +\sum_{i = 1}^N\Delta_s \: \gamma_{i, -, A} \: \gamma_{i, +, B}+\Delta_s \: \gamma_{i, -, B} \: \gamma_{i, +, A}. \end{aligned} \end{equation} \tag{ D7 }$$

In figure 9, we show the decomposition of the spinfull Kitaev chain into four Majorana chains and show the various pairings of Majoranas caused by different terms of the Hamiltonian as solid bonds. In all the figures, the blue spheres with + represent $\gamma_{i, +, A}$, the pink spheres with + represent $\gamma_{i, +, B}$, the blue spheres with − represent $\gamma_{i, -, A}$, the pink spheres with − represent $\gamma_{i, -, B}$. The green bonds in figure 9(a) represent the terms $\mu_+$ and $\mu_-$, the orange bonds in figure 9(b) represent the terms $(t_+ - \Delta_+)$ and $(t_- - \Delta_-)$, the violet bonds in figure 9(c) represent the terms $(t_+ + \Delta_+)$ and $(t_- + \Delta_-)$, the black bonds in figure 9(d) represent the term $\Delta_s$.

Zoom In Zoom Out Reset image size

Next, we study the spinfull Kitaev chain by setting $\mu_- = 0$, $\Delta_+ = -\Delta_- = 1$, $t_+ = t_- = 1$, $\Delta_s = 0.2$ and we keep $\mu_+$ as the varying parameter. The motivation of varying $\mu_+$ comes from the nanowire model where the + band's effective potential is varied through B. Similarly, we set $\Delta_+ = -\Delta_-$ to maintain the opposite parity of the p-wave couplings, just as in the nanowire case. The reason for keeping $\mu_- = 0$ shall be made clear soon. At $\mu_+ = 0$ both chains are in the topological regime and have two unpaired Majoranas at their ends. However, due to the presence of an s-wave inter-chain coupling $\Delta_s$, the Majoranas at each end pair up to form a normal fermion, rendering the spinfull Kitaev chain into a trivial phase. This can be seen in figure 10(a), where the yellow unpaired Majoranas hybridize through the solid black bonds caused by $\Delta_s$.

Zoom In Zoom Out Reset image size

Figure 10(b) represents an intermediate stage where $\mu_+$ is finite and can be seen as the green bonds in the top two chains. Figure 10(c) depicts the $\mu_+ = \infty$ case. Here, $\mu_+$ acts as the only dominating bond in the + chain and decouples the chain + chain from the − chain. The + chain goes into a trivial regime while the − chain stays in the topological regime, allowing there to be effectively two Majoranas at the ends of the wire. This is what precisely happens in the nanowire system, where the effective onsite potential being controlled by the magnetic field B of one band tunes it towards the topological regimes while tunes the other band towards the trivial regime. Before the transition, there are two p-wave superconductors hosting Majoranas. However, their interaction effectively leads to none at the edges. At the topological transition, one band becomes trivial, leading to two unpaired Majoranas at the ends.

Figure 11 shows the entanglement entropy plots for the spinfull Kitaev chain under the settings described above. Figure 10(a) shows the BEE plot for a 100-site chain while figure 11(b) shows the DEE plot for the same chain. Comparing figures 3(a) and 11(a) we see very similar features for both the spinfull Kitaev chain and the nanowire with no quantization and negligible change in the entanglement entropy value over the phase transition. Similarly, comparing figures 3(b) and 11(b) we see that the DEE is zero in the trivial regime and quantized to $ln2$ in the topological regime in both the spinfull Kitaev chain and the nanowire. A DEE value of $ln2$ in the topological regime is logically expected in the spinfull Kitaev chain since chains are now independent of any inter-chain interaction, the + Kitaev chain, being in its trivial regime, contributes 0 units to the DEE while the − Kitaev chain, being topological, contributes $ln2$ units [31] to the DEE.

Zoom In Zoom Out Reset image size