A proposal to extract and enhance four-Majorana interactions in hybrid nanowires

We simulate the smallest building block of the Sachdev-Ye-Kitaev (SYK) model, a system of four interacting Majorana modes. We propose a 1D Kitaev chain that has been split into three segments, i.e., two topological segments separated by a non-topological segment in the middle, hosting four Majorana Zero Modes at the ends of the topological segments. We add a non-local interaction term to this Hamiltonian which produces both bilinear (two-body) interactions and a quartic (four-body) interaction between the Majorana modes. We further tune the parameters in the Hamiltonian to reach the regime with a finite quartic interaction strength and close to zero bilinear interaction strength, as required by the SYK model. To achieve this, we map the Hamiltonian from Majorana basis to a complex fermion basis, and extract the interaction strengths using a method of characterization of low-lying energy levels and then finding the differences in energies between odd and even parity levels. We show that the interaction strengths can be tuned using two methods - (i) an approximate method of tuning overlapping Majorana wave functions (without non-local interactions) to a zero energy point followed by addition of a non-local interaction, and (ii) a direct parameter space optimization method using a genetic algorithm. We propose that this model could be further extended to more Majorana modes, and show a 6-Majorana model as an example. Since eigenspectral characterization of one-dimensional nanowire devices can be done via tunneling spectroscopy in quantum transport measurements, this study could be performed in experiment.


I.A. Our goal
Our goal is to design the building block of the Sachdev-Ye-Kitaev (SYK) model realizabe in an experimental setup, i.e., a single four-Majorana Zero Mode interacting system. This is inspired by the idea of an experimental device that implements the SYK Hamiltonian [1]. However, our motivation is not to build an SYK model, rather it is to detect the presence of four-body Majorana interactions in experiments, and to construct a method which enables us to enhance and extract the interaction strength in an experimental setup. The key initial steps to constructing such a device are: (1) finding the experimental signatures of a quartic interactions between the MZMs and (2) tuning the device to minimize bilinear terms coupling pairs of MZMs with respect to four-body interaction strength. Here, we theoretically consider a simple setup composed of a Kitaev chain nanowire hosting four MZMs, all interacting with each other via twobody and four-body interactions. First and foremost we develop a method of characterizing the interaction strengths using the eigenspectrum of our Hamiltonian. Next, we explore the conditions under which we can enhance the quartic interactions and suppress the bilinear interactions. Finally, lookinf ahead but not as our main focus, we consider how to scale up the setup to N=6 Majoranas.

I.B.1. What is the SYK model?
SYK model is a 0+1 dimensional model of N γ MZMs with random, all-to-all four body interactions. Kitaev [1] proposed this Majorana-based model as a variant of the original SY model of Sachdev and Ye [2] that described spins with random all-to-all couplings. The Hamiltonian for this model is: where γ i 's are the Majorana operators that obey the canonical anti-commutation relations: J ijkl are independent and identically distributed (i.i.d.) real random numbers drawn from a Gaussian distribution with zero mean and standard deviation given as:

I.B.2. Significance of the SYK model
There are many aspects that make the SYK model particularly interesting [1,3,4]. It is a strongly interacting model with symmetry properties that closely resemble quantum gravity. A remarkable property of this model is that it is maximally chaotic in the large N γ limit. Black holes scramble information at a fast rate, i.e., they are maximally chaotic, characterized by the upper limit of the Lyapunov exponent. It has been shown theoretically that the four-point out of time ordered correlators (OTOC) [5] in the SYK model also saturate the upper bound on the Lyapunov exponent [6]. Thus, from a quantum chaos perspective the SYK mimics the physics of black holes, i.e., it has holographic duality to black hole physics [7]. Besides this, the SYK model being exactly solvable at large N γ limit presents itself as a great tool towards simplifying understanding of the physics of quantum gravity in general. All these combined makes SYK the perfect toy model for engineering experimentally realizable black hole models.

I.B.3. Theoretical proposals towards experimental SYK devices
There have been a couple of proposals in the recent past that address challenges towards a practically realizable Sachdev-Ye-Kitaev (SYK) device which include (i) designing a system with a large number of MZMs with random interactions and (ii) formulating a method to suppress unwanted bilinear interactions that come along with the quartic interaction in a large-N γ system. One set of proposals features multiple nanowires coupled via a central quantum dot [8,9], where they use a time reversal symmetry argument to suppress the bilinear terms. In another proposal, MZMs are coupled around a vortex by engineering a nanoscale hole in the superconducting film on the surface of a three-dimensional topological insulator [10]. At the neutrality point, it is hypothesized that the bilinear terms are suppressed. A graphene based model [11] has also been put forward. Other than condensed matter systems, there are have been attempts to realize SYK in ultracold atoms [12], optical lattice systems [13], nuclear spins [14]. Digital quantum computers were also proposed to simulate the SYK model [15][16][17].

II. THE TWO-COMPLEX-FERMION MODEL
In this section we show that by characterizing the energy level spacings between the odd and even parity states we can back out the interactions that are present in the system. This is the first step for tuning the system to the SYK point which also requires control over the interaction strength (see Section IV).
Our starting point is the two-complex-fermion model that hosts four MZMs. We analyze this model with the goal of establishing the relationship between the fermionic spectrum that can be probed in transport experiments and the interaction terms in MZM Hamiltonians.
The model consists of two fermion operators, each of which can be decomposed into a pair of MZMs. In terms of the four MZM operators γ 1 − γ 4 , the interactions between the four MZMs can be of two types -(i) six different two-body, or bilinear, interactions K ij and (ii) a single four-body , or quartic, interaction J 1234 that includes all four MZMs in this model (Fig. 1). The Hamiltonian of this model is: The representation that we have used for γ's in this paper is described in Appendix A. At the same time, the spectrum of the two-complex fermion system can also be described by the Hamiltonian: where n 1 = (f † 1 f 1 ) and n 2 = (f † 2 f 2 ) are the quasiparticle number operators, λ 1 and λ 2 are the quasiparticle energies and u is the interaction strength of the two quasiparticles (and we have dropped the constant term as we are only interested in energy differences). We note that there are multiple different complex fermion Hamiltonians that result in the same spectrum that are connected by means of Bogoliubov transformations. We choose the representation of Eq. (5) because this Hamiltonian is in the diagonal form. We establish the relation between the MZM and the complex fermion representations of the Hamiltonian in Appendix B.

II.A. Energy spectra
In this section, we investigate the effect of the bilinear and quartic interaction terms on the spectral properties of the Hamiltonian. From the complex fermion representation Eq. (5), it is clear that the Hamiltonian is parity preserving with four energy levels, a pair of odd parity and a pair of even parity states. Switching focus to the MZM representation Eq. (4), we observe that tuning the bilinear and the quartic interaction terms results in a shift of the even and odd energy levels. Crucially, the shifts in the energy levels depend on which terms were tuned. In Fig. 2 we show some examples of how tuning terms in the MZM representation of the Hamiltonian affects the energy level spacings between odd (shown with black dashed lines) and even (shown with red solid lines) FIG. 1. Two fermion sites (ovals) with four Majoranas γ1-γ4 (circles). All possible bilinear terms are shown as lines above, and the single quartic term as the line below the circles. parity states. We note that in Fig. 2 we label the states by their complex fermion occupation numbers using the connection between representations from Appendix B.
In the top row we show the plots for the dependence of the energy levels on J 1234 with K ij fixed. In Fig. 2(a), we set K ij = 0, and observe the splitting of states of different parity as a function of J 1234 , while the same parity states stay degenerate throughout. This shows that J 1234 causes repulsion between odd and even parity energy levels. In the consecutive plots in Fig. 2(b), (c) and (d), as we increase K ij from 0.5 to 50, the strength of K ij becomes progressively more dominant over the strength of J 1234 and the energy level spacings between same parity states become significantly greater than that of different parity states. Similarly, in the bottom set of figures we plot the dependence of energy levels on K ij with J 1234 fixed in each plot. In Fig. 2(e) we set J 1234 = 0, and observe that states of the same parity split as function of K ij . This shows that K ij causes repulsion between same parity energy levels, and in the consecutive plots in Fig. 2(f), (g) and (h), as we increase J 1234 from 0.5 to 50 and as the strength of J 1234 becomes dominant over K ij , the energy level spacings between different parity states increase significantly over that of the same parity states. We conclude that K ij causes repulsion between same parity energy levels whereas J 1234 causes repulsion between different parity energy levels. In the case of the SYK model, since there are only quartic interactions, i.e., only J 1234 term is non-zero, the energy levels should look like Fig. 2(a), where states of the same parity remain degenerate, while those of opposite parity split, with the energy gap set by J 1234 . Thus, the spectroscopy of the energy differences between the odd and even parity states could tell us which interactions are present in the system.
Note: The criteria for energy level crossings in an N cf > 2 complex-fermion SYK model might be different than what is shown in Fig. 2(a). For example, for an N cf = 3 complex-fermion SYK model, we get energy crossings/degeneracies in opposite parity states rather than the same parity states as in N cf = 2 case. This is discussed in Section VI.B (and shown in Fig. 9(b) and Fig. 10(b)).

III. EXTRACTING THE INTERACTION STRENGTH FROM TUNNELING SPECTROSCOPY
In this section we discuss the extraction of bilinear and quartic interaction strengths (K ij 's and J 1234 in Eq. (4)) from tunneling transport measurements on hybrid superconductor-semiconductor nanowire devices that host four Majorana zero modes. Generically, transport measurements involve adding or removing electrons from the device and hence these measurements can be used to determine the energy differences between states of different parities, which, in turn, can be used to reconstruct the energy level spectra like those in Fig. 2.
Therefore, our starting point is the spectrum of energy eigenvalues E |00 , E |01 , E |10 , and E |11 . The eigenstates are labeled by the occupancy of the two quasiparticle states in the complex fermion representation. For example, the state |00 corresponds to both states being empty, while the state |01 corresponds to the first state empty and the second quasiparticle state occupied. Following the discussion in appendix B, we can relate the parameters of the Hamiltonian in Eq. 5 to the energy eigenvalues via: where 0 is the overall offset of the energy eigenvalues.
Since we are only interested in the energy level differences, identifying each eigenstate by their quasiparticle occupancies becomes redundant. Hence, we refer to each energy level by their parities. From hereon, we will de- where E e 's and E o 's are the quasiparticle energy levels corresponding to the even and odd parity states respectively). The odd and even energy levels can be used interchangeably amongst each other as long as the definition is followed throughout in all the equations.
Inverting the relation in Eq. (6), we can use the spacings between the even and odd parity energy levels to find the Hamiltonian (Eq. (5)) parameters , i.e., From the derivations in Appendix B, we also show that u gives us a direct relation to the quartic interaction strength J 1234 , but we find that λ's do not have a oneto-one relation with K ij 's (as seen from Eq. (B5)), i.e, a multitude of possible values for K ij 's can lead to the same values of λ's. However we find that it is possible to set a bound on K ij i.e, make all K ij = 0 by setting all λ i = 0. This allows us to find an SYK point at which only quartic interactions are present (a non-zero value of J 1234 and all K ij = 0). From Eq. (7), (8) and (9), we see that u is the difference in energies between the two even parity states and the two odd parity states, whereas λ's are made up of energy differences among even states and among odd states. A region with non-zero u but λ 1 , λ 2 = 0 will have degenerate E e 's and E o 's with some value of energy gap between them as shown in Fig. 2 (a). Therefore, from these relations we see that characterizing the eigenspectra can help us distinguish between the bilinear and quartic interaction strengths and thus guide us towards an SYK point in an experiment.

IV. KITAEV CHAIN WITH INTERACTIONS
To model an experimentally realizable form of the two complex fermion model, we introduce a 1D Kitaev chain model that hosts four MZMs. Specifically, our model consists of a quantum wire with N sites that is divided into three segments -two topological segments separated by a non-topological segment in the middle, as shown in Fig. 3. In order to induce four-MZM interactions, we supplement the Kitaev chain model with a non-local interaction term that is described in the next section. The total Hamiltonian thus consists of the Kitaev chain part H Kitaev-chain and the non-local interaction part H nl-int . In this section we describe the model. In the next section, we demonstrate that it is possible to tune this model to the point where bilinear interactions become zero while the quartic interaction remains finite.

IV.A. Kitaev chain
The Hamiltonian of a spinless p-wave Kitaev chain of length N is: where c † j , and c j are the complex fermion creation and annihilation operators on site j. The physical parameters governing this system are the site-dependent electrochemical potential µ, hopping amplitude t, and the superconducting pairing field ∆.
As shown in Fig. 3, the first segment of the wire (labeled t1) runs from site 1 to site n nt1 − 1 and is biased to the electrochemical potential µ t1 . The second segment (labeled nt) runs from site n nt1 to site n nt2 , with electrochemical potential µ nt . The final segment (labeled t2) runs from site n nt2 + 1 to site N with electrochemical potential µ t2 . The two topological segments, t1 and t2, have their electrochemical potentials set to ensure that they are in the topological phase: |µ t1 | < 2|t| and |µ t2 | < 2|t|; while the non-topological segment, nt, has its electrochemical potential set to |µ nt | > 2|t| to ensure that it is in the trivial phase. MZMs appear at the interface between topological and non-topological segments as indicated in Fig. 3 (here, the vacuum at the ends of the wire can be regarded as trivial).
FIG. 3. A Kitaev chain nanowire separated into three segments -two topological segments (red) and a non-topological segment (white) in between. The total number of sites is N and the non topological segment ranges between sites nnt 1 to nnt 2 . There are four MZM's (yellow circles) on this wire -two pairs localized at the edges of the two topological segments.

IV.B. Four-Majorana coupling from non-local interactions
Majorana Zero Modes in the Kitaev chain model are characterized by zero energy states separated from the bulk states by an energy gap. The wave functions corresponding to the Majorana modes have an oscillatory behavior with decaying amplitude. Depending on parameters like the length of the wire (N ) and other parametersµ, t and ∆, the MZMs can either be localized at the very ends of the wire (under the condition t = ∆) or spread out into the inner sections of the wire. Interactions between MZMs occur when these wave functions overlap with each other. As the overlap between distant MZMs is typically smaller than between nearby ones, this type of overlap tends to induce bilinear interactions between neighboring MZMs. To introduce a quartic interaction between MZMs, we introduce a non-local interaction term in the Hamiltonian with interaction strength U : This form of interaction is meant to model long-range interactions mediated, e.g. by charge or by phonons [18,19]. This interaction could be of Coulomb/densitydensity origin. It may be possible to enhance the interaction strength by taking certain measures like modifying the device design so as to minimize screening effects from nearby metals. Another idea is to modify the geometry of the device such as instead of having one nanowire with three segments (topo, non-topo, topo), which could have limited long-range interactions, we could have two nanowires (topo segments), with MZMs at their ends, placed parallel to each other, shown in Fig. 4.
In addition to quartic interactions, the non-linear term also induces additional bilinear interactions between the MZMs. In the following, we show that it is possible to eliminate the bilinear interactions by means of tuning the Hamiltonian parameters.

IV.C. Connection between the two-complex-fermion model and the Kitaev chain model
The eigenspectrum of the Hamiltonian of the two complex fermion model consists of four energy levels, that are separable in parity sectors as the total Hamiltonian is parity preserving. Tuning the interaction strengths relative to one another affect the energy level spacings specified by their parities as shown in Eqs. (7), (8), (9).
In the Kitaev chain model, the energy levels corresponding to the MZMs are the lowest four energy levels separated from the bulk states by an energy gap. The energy level structure of these four lowest energy levels is dependent on the Hamiltonian parameters. By applying equations (7), (8), (9) to these four lowest energy levels we can construct an effective low-energy model and extract its interaction parameters in terms of the MZM representation of the Hamiltonian Eq. (4).

V. REDUCING BILINEAR INTERACTIONS WHILE ENHANCING THE QUARTIC INTERACTION STRENGTH
Once we are able to identify the interactions between the MZMs, our goal becomes tuning the system to an SYK point at which the bilinear interactions become zero while the quartic interaction remains finite. In order to achieve this goal we tune the Hamiltonian parameters of our system (defined by µ t1 , µ t2 , µ nt , t, ∆ and U ) in order to zero out the bilinear interactions, i.e. |λ 1 |, |λ 2 | 0, while at the same time maximizing the value of quartic interaction |u|.
In this section, we first sweep the system parameters and show the existence of multiple, approximate SYK points. Next, we use an advanced optimization algorithm (the hybrid Genetic Algorithm described below) to locate SYK points at which bilinear interactions become essentially zero.

V.A. Existence of approximate SYK points
To find approximate SYK points, we begin with the 1D Kitaev chain model of Eq. (11) with no long-range interactions. When the wave functions of the MZMs overlap, the associated energy levels become non-zero. However, as the MZM wave functions are oscillatory, it is possible to tune the overlap integral to be zero in certain parametric regime/points where the overlapping wave functions cancel each other out [20]. We search for such optimal points where the lowest two quasiparticle energies are close to zero (i.e. the lowest four many-body states are degenerate) despite the MZM wave function envelopes having significant overlap. An example for this has been shown in Appendix D where we plot MZM wave functions and show how they overlap at an optimal point (shown in Fig. 13). After finding one of these optimal points, we add the non-local interaction term and show that by sweeping across a range of non-local interaction strength U , we can find approximate SYK points.
Therefore, we optimize for the condition |λ 1 |, |λ 2 | = 0 for the Hamiltonian in Eq. (10) at U = 0, i.e., without the presence of the non-local interaction which suppresses all nearest neighbor Majorana interactions. We then show that if we add the non-local interaction term U at these special optimized points, we can effectively have a dominant quartic interaction strength.

V.A.1. Obtaining the optimal points
We show in the Appendices that there are points in the parameter space of µ t1 , µ t2 , µ nt , t and ∆ at which we can minimize hybridization energy of the overlapping MZMs, i.e, tune their overlap integral to zero. We discuss the role of these parameters for optimization in Appendix C (see, Fig. 12 for a summary). At these optimal points the lowest four many-body energy levels are close to degenerate. Hence, we perform a numerical search for points where E |11 − E |00 is minimized. The values for µ t1 , µ t2 , µ nt , t and ∆ are obtained using a global search algorithm which we discuss in details in Appendix E. The optimal parameters obtained for N = 10 complex fermions and non-topological regionn nt1 = 4, n nt2 = 7 are: t = 0.4025, ∆ = 0.2167, µ t1 = 0.4832, µ t2 = 0.4832, µ nt = 8.5364.

V.A.2. Sweeping interactions to find approximate SYK points
Having minimized the local bilinear interactions by tuning the system to an optimal point in the parameter space of {µ t1 , µ t2 , t, ∆}, we then add the non-local interaction term H nl-int to induce quartic interactions between the MZMs. However, this non-local interaction term can also introduce additional non-local bilinear interactions besides the quartic interaction. Thus, we further optimize to cancel out these additional bilinear terms by tuning the non-local interaction strength U to the SYK point, i.e., we sweep across a range of values for U to find a point where |λ 1 |, |λ 2 | 0 and |u| is at a maximum value. This is shown in Fig. 5(a), where we fix the values of {µ t1 , µ t2 , t, ∆} and sweep the value of U . In this figure we see that λ 1 (solid yellow line), λ 2 (dashed blue line) and u (solid purple line) have several local maxima and minima at various points of U . Both λ 1 and λ 2 follow the same trajectory, i.e, their extrema coincide, and u, having a different trajectory, has extrema at different points. At U = 0.0733 and U = 0.2118, λ's have minima and u has maxima, i.e, they are the SYK points (shown by vertical dotted line-cuts). The best SYK point characterized by the greatest value for u λ has been shown in this figure with a dotted black line-cut at U = 0.212 at which we get λ 1 , λ 2 8.28e −5 and u 0.036 ( u λ 435). In Fig. 5(b), we plot the energies of the lowest four states, which were used to construct Fig. 5(a), as a function of U . We see that at the best SYK point (shown with the black dotted line-cut, zoomed-in in the inset) the even parity (red lines) and odd parity (solid black lines) energy levels are degenerate, i.e, , with a finite energy gap between them, i.e. E |11 + E |00 − E |10 + E |01 = 0. Following the discussions in Section II.A, this level structure is similar to the one displayed in Fig. 2 (a) with K ij = 0 and J 1234 = 0. Thus, from the low energy level spectral point of view, this point matches our expectation for an SYK point.
Other than the non-local interaction term Eq. (12), we have also explored other interaction terms but could not find an SYK point as we swept through a range of the interaction strength U at the optimal point. This is discussed in Appendix G and some examples are shown in Fig. 16). We hypothesize that our inability to tune systems with alternative interactions to an SYK point is due to the more local nature of the alternative interactions.
Note: There is a constraint on the feasible range of U that we have access to in the optimization process. This is based on the condition that MZM's appear in the topological regime separated by an energy gap from the bulk states ( ∆). Specifically, varying U beyond a certain bound results in the closing of this energy gap and the penetration of the MZM states into the continuum of bulk quasiparticle states. In order to detect this possibility, as we tune U , we check the total parity of the four lowest energy (many-body) states. The total parity is expected to be zero as there are two even and two odd parity states. However as the MZM states penetrate into the continuum, the parity of the lowest energy states becomes random and we know that we have exceeded the valid range of U .

V.B. Using a search algorithm to find SYK points
Another way to look for SYK points is to implement a direct search for |λ 1 |, |λ 2 | 0 and a maximum |u|, through the entire parameter space. We intend to not only look for better optimized points but also to compare with the results of the previous subsection in which we tuned for zero MZM wave function overlap integral.
Our optimization problem now consists of a five parameter space {µ t1 , µ t2 , t, ∆, U } with multiple parametric constraints (|µ t1 | < 2|t|, |µ t2 | < 2|t|, |µ nt | > 2|t|) and constraints in objective function. Hence, we attempt to solve this problem using an advanced optimization technique -a hybrid Genetic algorithm (a Genetic Algorithm search for roughly locating the SYK points, followed by a derivative-based search for refining the location of the SYK points). Our objective is to maximize |u| under the constraints of (i) |λ 1 |, |λ 2 | 0, and (ii) total parity of lowest four eigenstates equals zero. This search yields op- timized results for all the parameters µ t1 , µ t2 , µ nt , t, ∆ and U corresponding to the best SYK point in the given search range (discussed in details in Appendix F). For N = 10 complex fermions with non-topological regionn nt1 = 4, n nt2 = 7, the values obtained are t = 0.3229, ∆ = 0.1, µ t1 = 0.5871, µ t2 = 0.5871, µ nt = 6.3944 and U = 0.2254. In Fig. 6(a), we plot λ 1 , λ 2 , and u as a function U to show a comparison with the results in Fig. 5(a) (from the previous method). We fix the values of {µ t1 , µ t2 , t, ∆} and sweep the value of U so as to intersect the SYK point found by the hybrid Genetic Algorithm at U = 0.2254. We see that λ's and u follow similar trajectory in both the figures with the extrema of λ's located at different points than from that of u, and we can indeed find three SYK points at U = 0.0176, 0.0891, and 0.2254. The best SYK point can be seen at U = 0.2254 (shown by black dotted line cut) where u λ 530. Likewise, in Fig. 6 Fig. 5(b). We note that the SYK point obtained by this method (Fig. 6) yields a better result (higher value for u λ ) than the previous case (as shown in Fig. 5).

V.B.1. Adding another parameter for optimization
In the four-Majorana Hamiltonian (Eq. (4)) we have seven independent terms: six bilinear terms (K ij ) and one quartic term (J 1234 ). In the Kitaev chain Hamil-tonian, it is reasonable then to expect that we find a better optimal point by expanding our parameter space to seven by adding another independent variable. For this we have set the tunneling amplitude parameter t c at the center of the wire different from t e at the edges as two independent parameters. Then, we implemented the same hybrid Genetic Algorithm to find SYK points, and plot the results in Fig. 7. In Fig. 7(a) we see that λ's and u follow quite a different trajectory than in the previous figures ( Fig. 5(a), Fig. 6(a)) but we are still able to find an SYK point at U = 0.3477 where u λ ≈ 60, 000. By this measure, adding a seventh parameter gives a better SYK point than the previous two cases. Similar to the previous cases, by plotting the unprocessed energy levels in Fig. 7

VI.A. Why do we extend the model?
In the sections above we deal with characterizing and enhancing the quartic interaction strength within four-Majorana models. However, in the full SYK model, we need to have a large numer of Majorana zero modes with multiple quartic interaction terms. Hence, to characterize and measure multiple quartic interactions, we need to extend our model from N γ = 4 to higher N γ . As a proof of concept of our model being extendable to higher N γ , , We plot the lowest four energy levels vs U labelled by their parities (even-red, odd-black). At the SYK point obtained (shown in inset), we find that the two even and the two odd parity states are degenerate (level-crossings) while the energy gap between between different parity states is at a maximum. 1.42e −6 and u 0.085. (b), We plot lowest four energy levels vs U labelled by their parities (even-red, odd-black). At the SYK point obtained (shown in inset), we find that the two even and the two odd parity states are degenerate (level-crossings) while the energy gap between between different parity states is at a maximum.
we try applying the methods discussed in the sections above to a six-Majorana model.

VI.B. The three-complex-fermion model
This model consists of three sites, each with two Majoranas, i.e. six Majoranas in total. There can be three kinds of interaction terms between six Majoranas, they FIG. 8. A Kitaev chain nanowire separated into five segments -three topological segments (red) and two non-topological segment (white). The total number of sites is N and the non topological segments range between sites nnt1 1 to nnt1 2 for non-topological segment 1, and nnt2 1 to nnt2 2 for nontopological segment 2. There are six MZM's (yellow circles) on this wire -three pairs localized at the edges of the three topological segments.
are -(i) 15 bilinear interaction terms, (ii) 15 quartic interaction terms and (iii) 1 sextic interaction term -which is the new term that appears with six MZMs. In the MZM representation, the Hamiltonian of this model is: where J 123456 is the sextic interaction strength, J ijkl 's are the quartic interaction strengths and K ij 's are the bilinear interaction strengths.

VI.C. Extracting the interaction strengths
When we diagonalize the Hamiltonian in Eq. 13, we get eight eigenvalues corresponding to the 2 3 × 2 3 Hilbert space. The eigenstates are in the form |n 1 n 2 n 3 where n 1 , n 2 , n 3 = 0 or 1 for empty or filled fermion quasiparticle states respectively. As an extension to the four Majorana model, the Hamiltonian can be written in the quasiparticle basis as: where λ's are the bilinear interaction strengths, u's are the quartic interaction strengths and v is the sextic interaction strength in the quasiparticle basis. Following the arguments as in Appendix B and Section III, we can conclude that in a similar manner for the six MZMs case the energy levels of states |n 1 n 2 n 3 are related to interaction strengths λ's, u's and v as: where E is the column matrix of energy levels for states |n 1 n 2 n 3 , where n 1 , n 2 , n 3 = 0, 1 in the sequence as shown below: where I is the column matrix of interaction strengths λ's (bilinear), u's (quartic) and v (sextic), i.e.
and A is the matrix transforming E to I : Thus we can extract the interaction strengths by inverting this relation to solve for I:

VI.D. Kitaev chain model with interactions
Our goal is to model a Kitaev chain quantum wire (similar to the four MZMs case) that generates six MZMs. It is an N site chain that we divide into five segments -three topological segments separated by two non-topological segments in between two topological segments (toponontopo-topo-nontopo-topo). Following the conditions for topological phase we have |µ t1 | < 2|t|, |µ t2 | < 2|t| and |µ t3 | < 2|t| for the topological segments t1, t2, t3, and |µ nt1 | > 2|t|, |µ nt2 | > 2|t| for the non-topological segments nt1, nt2. As a result, we get six MZMs at the ends of the three topological segments. This has been shown schematically in Fig. 8.
To introduce the quartic interactions in the MZMs, we add the same non-local interaction term as we did in the four MZM case which is:  The total Hamiltonian after adding up the interaction term is: This also introduces additional non-local bilinear and sextic interactions in between the MZMs which we can suppress as shown in the following sections.

VI.E. Search for SYK points for dominant quartic interaction strength and suppressing bilinear and sextic interaction strengths
Following the same approach as in Section V, we enhance the quartic interaction strength using two different methods. (i) By optimizing the MZM wave function overlap integral to zero, i.e., searching for global minima for E |111 − E |000 in the parameter space of {µ t1 , µ t2 , µ t3 , µ nt1 , µ nt2 , t, ∆}. Then by sweeping the value of U , we search for minima of bilinear and sextic interaction strengths and maxima of quartic interaction strengths, i.e., maximize |u 12 |, |u 13 |, |u 23 | and set |λ 1 |, |λ 2 |, |λ 3 |, |v| 0, within a feasible search range such that the energy gap between the low energy states and bulk states is maintained. (ii) By performing a direct search for SYK points within the total parameter space of µ t1 , µ t2 , µ t3 , µ nt1 , µ nt2 , t, ∆ and U . We use a hybrid Genetic Algorithm search to optimize for |λ 1 |, |λ 2 |, |λ 3 |, v 0 and maximize |u 1 |, |u 2 |, |u 3 | with the constraint that the energy gap between the MZM states and bulk states doesn't close.
Using method (i) we suppress all local bilinear interactions by finding optimal points where MZM wave function overlap integral is zero. We find optimal points where E |111 − E |000 0 using a global search algorithm similar to that in four MZM case (described in details in Appendix E.1), then upon sweeping through U , we can find SYK points as shown in Fig. 9 where all interactions except for the quartic interactions are well suppressed. In Fig. 9(a) The energy levels used to compute the parameters λ 1,2,3 , u 12,23,13 , and v in Fig. 9(a) are plotted as a function of U in Fig. 9(b). In this figure, there are eight energy levels with some degenerate levels. The black dotted line-cut corresponds to the SYK point. Focusing on the line-cut in the inset, the bottom most line (red) consists of three degenerate even parity energy levels and the second line from bottom (black) consists of three degenerate odd parity energy levels, the two intersecting lines on the top are a single odd and a single even parity energy level each. Thus, at the SYK point we observe energy crossings between even and odd parity levels, in contrast to the four MZM case, in which states of the same parity were crossing.
Using method (ii) we find more SYK points using a hybrid Genetic Algorithm (described in detail in Appendix F.1). In Fig. 10 we plot the interaction strengths as a function of U to capture the SYK point found by the Genetic Algorithm at U = 0.097. Here, u λ 2 (λ: average of λ 1,2,3 , u = u 12 = u 13 = u 13 ) and u v 10.5. For completeness, we plot the energy levels that we used to extract λ 1,2,3 , u 12,23,13 , and v (in Fig. 10(a)) in Fig. 10(b). The level structure is analogous to the one previously found in Fig. 9(b). In the inset, the bottom two red and black lines are single levels whereas the top line consists of three degenerate even (red) and odd parity levels (black) crossing each other.
We note that the eigenstates are shuffled at the energy level crossings and we need to reorder them in order to maintain consistency of the definitions of λ's, u's and v at different search intervals/points. This was done using an eigenvalues reordering algorithm [21].

VII. FUTURE RELEVANCE
From this study we show that it is indeed possible to design a low N γ (particularly N γ = 4, 6) SYK model in a 1D nanowire system. It might be possible to extend this model to a N γ > 6 system in future work. We show that it is possible to analyze the experimentally accessible spectra of eigenstates with quantum transport measurements, and use this information to assess the strength of bilinear and quartic interaction terms.
In hybrid superconductor-semiconductor nanowire devices, some of the parameters that we use in constructing the multi-segment Kitaev chain model are tunable. For instance, the tunneling amplitudes and chemical potentials can be tuned with gates. Other parameters, such as the induced gap, may be harder to tune in situ, though it is possible in principle. The biggest anticipated challenge is to tune the interaction strength U which may be severely constrained by the device geometry. This crucial term may also turn out to be too small, thus closing the door to future work. Though it can possibly be enhanced by careful design of nanowire devices.

VIII. EXPERIMENTAL PROTOCOL FOR A 4-MZM INTERACTION DEVICE
We propose a semiconductor-superconductor nanowire device with tunnel probes connected to the ends of the wire acting as source and drain channels across which a voltage bias is applied for performing tunneling spectroscopy. This device can be fabricated to have multiple local gates in contact/close to the nanowire that can be individually controlled to tune the system parameters such as µ and t corresponding to different segments of the wire. Other parameters like ∆ can be pre-selected by choosing appropriate material for the superconductorsemiconductor nanowire. We can vary U by playing with the geometry of the device, e.g., by having two different nanowires, that host MZMs, placed parallel to each other, by minimizing screening from nearby metals and by tuning the electron density.
Differential conductance data will tell us the energy level spacings and thus can be used to construct the lowest four quasiparticle eigenspectra similar to Fig. 2. Then, by using Eq. (7), (8) and (9) we can extract the values of λ 1 , λ 2 and u for a particular set of conditions. Varying the system parameters will give several sets of spectroscopy data for different conditions, which can be further analyzed to tune the system to SYK point, i.e., zero λ's and a maximum u. We can vary a couple of parameters at a time to obtain a map of energy level spacings as shown in Fig. 12 , which will give us an estimate of the parametric space for the approximate SYK points. Then we can further narrow down our search by tuning individual parameters and obtaining the interaction strengths as a function of these parameters, as shown in Fig. 5(a), 6(a), 7(a) where we plot λ's and u as a function of U . Alternatively, we can fix U and and obtain the interaction strengths as a function of other parameters as shown in Fig. 16. From this, we can also conclude that it is not necessary to tune all the parameters simultaneously to arrive at the SYK point once we have a good parametric region to work with.

X. STUDY DESIGN
This study was undertaken in a period of about nine months. A summary of our study design is shown in Fig. 11.

XI. DATA AVAILABILITY
All codes are available on Zenodo [48].

XII. ACKNOWLEDGEMENTS
The authors thank J. Stenger for discussions of the SYK model. The matrices representing the γ operators are constructed by transforming the Hamiltonian into spin chain basis (similar to Jordan Wigner transformation) using the following algorithm.
In the matrix form, this can be written as where M is a 4 × 4 skew symmetric matrix with six unique non-zero entries, i.e, K ij 's. Eigenvalues of the skew-symmetric matrices come in pairs: ±λ 1 , ±λ 2 with complex raising and lowering operators f and f † as their eigenstates.
The eigenvectors of M define the basis transformation U from the MZM to the complex fermion representation, where U is an unitary basis transformation matrix from the MZM basis to the complex fermionic basis andM is the transformed bilinear interaction matrix The two quasiparticle energies are where s = K 2 12 + K 2 13 + K 2 14 + K 2 23 + K 2 24 + K 2 34 and d = K 14 K 23 − K 13 K 24 + K 12 K 34 . Thus, we can write the bilinear interaction Hamiltonian in this basis as where n i 's are the quasiparticle number operators. n i is 0 or 1 corresponding to the filled or empty quasiparticle states. Similarly, we can also write the quartic interaction in this basis as The total Hamiltonian in the complex fermion representation is therefore where {f i , f † j } = 2δ ij and f † i f i = 2n i . Thus, we can write the energy levels of this Hamiltonian corresponding to states |00 , |10 , |01 , |11 in the form |n 1 n 2 as where E e 's are the eigenvalues of even parity states, i.e, |00 and |11 and E o 's are eigenvalues of the odd parity states, i.e., |10 and |01 . 0 is a constant shift in the energies.
parameter (µ t 1,µ t 1,t,∆) and keeping the other parameters fixed at the SYK point mentioned in Fig. 6 (except for Fig. (c) where µ t1 = µ t2 are varied together), shown in Fig. 16. As seen from Fig. 16 (a), (b) and (c), we find SYK points as function of t, ∆ and µ t1 = µ t2 , however we couldn't tune to an SYK point as a function of µ t1 , µ t2 and µ nt when swept through independently, shown in Fig. 16 (d), (e) and (f). Hence we show that it's possible to find SYK points if we sweep through the gate-tunable Kitaev chain parameters and keep U fixed, since U can be difficult to tune in experiments.
Other interaction terms we checked are -(a)  FIG. 16. We plot interaction strengths λ1,λ2,u as a function of parameters of Hamiltonian in Eq. 11. The interaction strengths are plotted as a function of one parameter keeping the others fixed to the SYK point paramteric values mentioned in Fig. 6: We plot λ1,λ2,u as a function of Fig.(a). t, Fig.(b). ∆, Fig.(c). µt1 and µt2 where µt1 = µt2 are varied together, Fig.(d). µnt, Fig.(e). µt1, Fig.(d). µt2.