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

Tasnum Reza  Sergey M. Frolov  David Pekker  Department of Physics and
Astronomy, University of Pittsburgh, 15260, USA

###### Abstract

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.

††preprint: APS/123-QED

##  I Introduction

###  I.A Our goal

Our ultimate goal is to construct an experimental device that implements the
Sachdev-Ye-Kitaev (SYK) Hamiltonian [1], which consists of a collection of
four-body (quartic) interactions between Majorana zero modes (MZM), with each
four-body interaction having a random valued strength. The key initial steps
to constructing such an SYK 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 strengths. Here, we theoretically consider a simple setup
composed of a Kitaev chain nanowire hosting four MZMs, all interacting with
each other via two-body 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, we consider how to scale the setup to include more than
four MZMs.

###  I.B Context

####  I.B.1 What is the SYK model?

SYK model is a 0+1 dimensional model of $N_{\gamma}$ 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:

$H=\sum_{i<j<k<l}^{N_{\gamma}}J_{ijkl}\gamma_{i}\gamma_{j}\gamma_{k}\gamma_{l},$
(1)

where $\gamma_{i}$’s are the Majorana operators that obey the canonical anti-
commutation relations:

$\\{\gamma_{i},\gamma_{j}\\}=\delta_{ij},\qquad\gamma_{i}^{\dagger}=\gamma_{i}.$
(2)

$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:

$\overline{J^{2}_{ijkl}}=\frac{3!J^{2}}{N_{\gamma}^{3}},\qquad\overline{J_{ijkl}}=0.$
(3)

####  I.B.2 Significance of the SYK model

There are many aspects that make the SYK model particularly interesting [3, 4,
1]. 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_{\gamma}$ 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_{\gamma}$ 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_{\gamma}$ 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).

Figure 1: Two fermion sites (ovals) with four Majoranas
$\gamma_{1}$-$\gamma_{4}$ (circles). All possible bilinear terms are shown as
lines above, and the single quartic term as the line below the circles.

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
$\gamma_{1}-\gamma_{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:

$H=J_{1234}\prod_{i=1}^{4}\gamma_{i}+i\sum_{1\leq i<j\leq
4}K_{ij}\gamma_{i}\gamma_{j}.$ (4)

The representation that we have used for $\gamma$’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:

$\displaystyle H$ $\displaystyle=\lambda_{1}(2n_{1}-1)+\lambda_{2}(2n_{2}-1)$
(5) $\displaystyle-u(2n_{1}-1)(2n_{2}-1),$

where $n_{1}=(f_{1}^{\dagger}f_{1})$ and $n_{2}=(f_{2}^{\dagger}f_{2})$ are
the quasiparticle number operators, $\lambda_{1}$ and $\lambda_{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

Figure 2:  (Top row) Energy levels of the Hamiltonian in Eq. (4) as a function
of $J_{1234}$, fixing all $K_{ij}$’s equal to a constant - (a) $K_{ij}=0$. (b)
$K_{ij}=0.5$. (c) $K_{ij}=5$. (d) $K_{ij}=50$. (Bottom row) Energy levels as a
function of $K_{ij}$ (all $K_{ij}$’s are equal and varied simultaneously),
fixing $J_{1234}$ equal to a constant - (e) $J_{1234}=0$. (f) $J_{1234}=0.5$.
(g) $J_{1234}=5$. (f) $J_{1234}=50$.

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) 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_{\text{cf}}>2$ complex-
fermion SYK model might be different than what is shown in Fig. 2(a). For
example, for an $N_{\text{cf}}=3$ complex-fermion SYK model, we get energy
crossings/degeneracies in opposite parity states rather than the same parity
states as in $N_{\text{cf}}=2$ case. This is discussed in Section VI.B (and
shown in Fig. 8(b) and Fig. 9(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_{\ket{00}}$, $E_{\ket{01}}$, $E_{\ket{10}}$, and $E_{\ket{11}}$. The
eigenstates are labeled by the occupancy of the two quasiparticle states in
the complex fermion representation. For example, the state $\ket{00}$
corresponds to both states being empty, while the state $\ket{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:

$\begin{pmatrix}E_{\ket{00}}\\\ E_{\ket{01}}\\\ E_{\ket{10}}\\\
E_{\ket{11}}\end{pmatrix}=\begin{pmatrix}1&-1&-1&-1\\\ 1&-1&1&1\\\ 1&1&-1&1\\\
1&1&1&-1\\\ \end{pmatrix}\begin{pmatrix}\epsilon_{0}\\\ \lambda_{1}\\\
\lambda_{2}\\\ u\end{pmatrix},$ (6)

where $\epsilon_{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 denote
$E_{\ket{00}}$ as $E^{e}_{1}$, $E_{\ket{10}}$ as $E^{o}_{1}$, $E_{\ket{01}}$
as $E^{o}_{2}$, and $E_{\ket{11}}$ as $E^{e}_{2}$ (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.,

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle(-E^{e}_{1}+E^{e}_{2}+E^{o}_{1}-E^{o}_{2})/4$ (7)
$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle(-E^{e}_{1}+E^{e}_{2}-E^{o}_{1}+E^{o}_{2})/4$ (8) $\displaystyle
u$ $\displaystyle=$
$\displaystyle(E^{o}_{1}+E^{o}_{2}-E^{e}_{1}-E^{e}_{2})/4.$ (9)

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
$\lambda$’s do not have a one-to-one relation with $K_{ij}$’s (as seen from
Eq. (30)), i.e, a multitude of possible values for $K_{ij}$’s can lead to the
same values of $\lambda$’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 $\lambda_{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
$\lambda$’s are made up of energy differences among even states and among odd
states. A region with non-zero $u$ but $\lambda_{1},\lambda_{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

Figure 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 $n_{\text{nt}_{1}}$ to $n_{\text{nt}_{2}}$. There are four MZM’s
(yellow circles) on this wire - two pairs localized at the edges of the two
topological segments.

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

$H_{\text{tot}}=H_{\text{Kitaev-chain}}+H_{\text{nl-int}}.$ (10)

thus consists of the Kitaev chain part $H_{\text{Kitaev-chain}}$ and the non-
local interaction part $H_{\text{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:

$\displaystyle H_{\text{Kitaev-chain}}$ $\displaystyle=$
$\displaystyle-\sum_{j}\mu_{j}c_{j}^{\dagger}c_{j}-t\sum_{j}c_{j}^{\dagger}c_{j+1}+\text{h.c.}$
(11) $\displaystyle+$ $\displaystyle\Delta\sum_{j}c_{j}c_{j+1}+\text{h.c.},$

where $c_{j}^{\dagger}$, 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 $\mu$, hopping
amplitude $t$, and the superconducting pairing field $\Delta$.

As shown in Fig. 3, the first segment of the wire (labeled t1) runs from site
$1$ to site $n_{\text{nt}_{1}}-1$ and is biased to the electrochemical
potential $\mu_{\text{t1}}$. The second segment (labeled nt) runs from site
$n_{\text{nt}_{1}}$ to site $n_{\text{nt}_{2}}$, with electrochemical
potential $\mu_{\text{nt}}$. The final segment (labeled t2) runs from site
$n_{\text{nt}_{2}}+1$ to site $N$ with electrochemical potential
$\mu_{\text{t2}}$. The two topological segments, t1 and t2, have their
electrochemical potentials set to ensure that they are in the topological
phase: $|\mu_{\text{t1}}|<2|t|$ and $|\mu_{\text{t2}}|<2|t|$; while the non-
topological segment, nt, has its electrochemical potential set to
$|\mu_{\text{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).

###  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 - $\mu$, $t$ and $\Delta$, the MZMs can either be localized at the
very ends of the wire (under the condition $t=\Delta$) 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$:

$H_{\text{nl-int}}=U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{j}^{\dagger}c_{j}.$ (12)

This form of interaction is meant to model long-range interactions mediated,
e.g. by charge or by phonons [18, 19]. 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
$\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{nt}}$, $t$, $\Delta$ and
$U$) in order to zero out the bilinear interactions, i.e.
$|\lambda_{1}|,|\lambda_{2}|\simeq 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. 12). 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 $|\lambda_{1}|,|\lambda_{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
$\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{nt}}$, $t$ and $\Delta$ 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. 11 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
$\left|E_{\ket{11}}-E_{\ket{00}}\right|$ is minimized. The values for
$\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{nt}}$, $t$ and $\Delta$ 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 region - $n_{\text{nt}_{1}}=4$, $n_{\text{nt}_{2}}=7$ are:
$t=0.4025$, $\Delta=0.2167$, $\mu_{\text{t1}}=0.4832$,
$\mu_{\text{t2}}=0.4832$, $\mu_{\text{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
$\\{\mu_{\text{t1}},\mu_{\text{t2}},t,\Delta\\}$, we then add the non-local
interaction term $H_{\text{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
$|\lambda_{1}|,|\lambda_{2}|\simeq 0$ and $|u|$ is at a maximum value. This is
shown in Fig. 4(a), where we fix the values of
$\\{\mu_{\text{t1}},\mu_{\text{t2}},t,\Delta\\}$ and sweep the value of $U$.
In this figure we see that $\lambda_{1}$ (solid yellow line), $\lambda_{2}$
(dashed blue line) and $u$ (solid purple line) have several local maxima and
minima at various points of $U$. Both $\lambda_{1}$ and $\lambda_{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$,
$\lambda$’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 $\left|\frac{u}{\lambda}\right|$ has been shown in this
figure with a dotted black line-cut at $U=0.212$ at which we get
$\lambda_{1},\lambda_{2}\simeq 8.28e^{-5}$ and $u\simeq 0.036$
($\left|\frac{u}{\lambda}\right|\simeq 435$).

In Fig. 4(b), we plot the energies of the lowest four states, which were used
to construct Fig. 4(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, $\left|E_{\ket{11}}-E_{\ket{00}}\right|\simeq 0$,
$\left|E_{\ket{10}}-E_{\ket{01}}\right|\simeq 0$), with a finite energy gap
between them, i.e.
$\left|E_{\ket{11}}+E_{\ket{00}}\right|-\left|E_{\ket{10}}+E_{\ket{01}}\right|\neq
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}\neq 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. 13). 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 ($\simeq\Delta$). 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$.

(a)

(b)

Figure 4:  The optimal point parameters obtained for a $N=10$ site Kitaev
chain model (non-topological segment:
$n_{\text{nt}_{1}}=4,n_{\text{nt}_{2}}=7$) obtained following the method of
finding approximate SYK point by minimizing MZM wave function overlap to zero
as described in Section V.A are: $t=0.4025$, $\Delta=0.2167$,
$\mu_{\text{t1}}=0.4832$, $\mu_{\text{t2}}=0.4832$, $\mu_{\text{nt}}=8.5364$.
(a), We sweep across $U$ at the optimal point and find a couple of SYK points
at $U=0.0733$ and $U=0.2118$. At the best SYK point, i.e. at $U=0.2118$, we
have $\lambda_{1},\lambda_{2}\simeq 8.28e^{-5}$ and $u\simeq 0.036$. (b), We
plot the lowest four energy levels vs $U$ labelled by their parities (even-
red, odd-black). At the best SYK point obtained at $U=0.2118$ (shown in
inset), 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.

###  V.B Using a search algorithm to find SYK points

Another way to look for SYK points is to implement a direct search for
$|\lambda_{1}|,|\lambda_{2}|\simeq 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 six parameter space
$\\{\mu_{\text{t1}},\mu_{\text{t2}},t,\Delta,U\\}$ with multiple parametric
constraints ($|\mu_{\text{t1}}|<2|t|$, $|\mu_{\text{t2}}|<2|t|$,
$|\mu_{\text{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) $|\lambda_{1}|$, $|\lambda_{2}|$ $\simeq 0$, and (ii) total parity of
lowest four eigenstates equals zero. This search yields optimized results for
all the parameters $\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{nt}}$,
$t$, $\Delta$ 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 region - $n_{\text{nt}_{1}}=4$, $n_{\text{nt}_{2}}=7$, the
values obtained are $t=0.3229$, $\Delta=0.1$, $\mu_{\text{t1}}=0.5871$,
$\mu_{\text{t2}}=0.5871$, $\mu_{\text{nt}}=6.3944$ and $U=0.2254$.

In Fig. 5(a), we plot $\lambda_{1}$, $\lambda_{2}$, and $u$ as a function $U$
to show a comparison with the results in Fig. 4(a) (from the previous method).
We fix the values of $\\{\mu_{\text{t1}},\mu_{\text{t2}},t,\Delta\\}$ 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 $\lambda$’s and $u$ follow
similar trajectory in both the figures with the extrema of $\lambda$’s located
at different points than from that of $u$, and we can indeed find three SYK
points at $U=0.0176,0.0891,\text{ and }0.2254$. The best SYK point can be seen
at $U=0.2254$ (shown by black dotted line cut) where
$\left|\frac{u}{\lambda}\right|\simeq 530$. Likewise, in Fig. 5(b), we show
that at the SYK point $\left|E_{\ket{11}}-E_{\ket{00}}\right|\simeq 0$,
$\left|E_{\ket{10}}-E_{\ket{01}}\right|\simeq 0$ while
$\left|E_{\ket{11}}+E_{\ket{00}}\right|-\left|E_{\ket{10}}+E_{\ket{01}}\right|\neq
0$, similar to that in Fig. 4(b). We note that the SYK point obtained by this
method (Fig. 5) yields a better result (higher value for
$\left|\frac{u}{\lambda}\right|$) than the previous case (as shown in Fig. 4).

(a)

(b)

Figure 5:  The SYK point parameters obtained for a $N=10$ site fermion Kitaev
chain model (non-topological segment:
$n_{\text{nt}_{1}}=4,n_{\text{nt}_{2}}=7$) obtained using a hybrid genetic
algorithm search as described in Section V.B are: $t=0.3229$, $\Delta=0.1$,
$\mu_{\text{t1}}=0.5871$, $\mu_{\text{t2}}=0.5871$, $\mu_{\text{nt}}=6.3944$,
$U=0.2254$. (a), We sweep across $U$ fixing the other parameters so as to meet
the best SYK point at $U=0.2254$. We find three SYK points at $U=0.0176$,
$U=0.0891$ and $U=0.2254$. At the best SYK point, i.e. at $U=0.2254$, we have
$\lambda_{1},\lambda_{2}\simeq 9.42e^{-5}$ and $u\simeq 0.05$. (b), 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.

####  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 Hamiltonian, 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_{\text{c}}$ at the center of the wire different from $t_{\text{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. 6. In Fig.
6(a) we see that $\lambda$’s and $u$ follow quite a different trajectory than
in the previous figures (Fig. 4(a), Fig. 5(a)) but we are still able to find
an SYK point at $U=0.3477$ where $\left|\frac{u}{\lambda}\right|\approx
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. 6(b) we observe that at $U=0.3477$,
$\left|E_{\ket{11}}-E_{\ket{00}}\right|\simeq 0$,
$\left|E_{\ket{10}}-E_{\ket{01}}\right|\simeq 0$ while
$\left|E_{\ket{11}}+E_{\ket{00}}\right|-\left|E_{\ket{10}}+E_{\ket{01}}\right|\neq
0$.

(a)

(b)

Figure 6:  We add another parameter (tunneling amplitude - $t_{\text{c}}$ at
the center of the wire different than at the edges - $t_{\text{e}}$) for
optimization in order to find better SYK points as described in Section V.B.1.
The SYK point parameters for $N=10$ site Kitaev chain model (non-topological
segment: $n_{\text{nt}_{1}}=4,n_{\text{nt}_{2}}=7$) obtained using a hybrid
genetic algorithm search are: $t_{\text{e}}=0.9028$, $t_{\text{c}}=0.1166$,
$\Delta=0.1041$, $\mu_{\text{t1}}=0.3299$, $\mu_{\text{t2}}=0.3299$,
$\mu_{\text{nt}}=7.2691$ and $U=0.3477$. (a), We sweep across $U$ fixing the
other parameters to meet the SYK point at $U=0.3477$. At this point, we have
$\lambda_{1},\lambda_{2}\simeq 1.42e^{-6}$ and $u\simeq 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.

##  VI Extension to a six MZM model

###  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_{\gamma}=4$ to higher
$N_{\gamma}$. As a proof of concept of our model being extendable to higher
$N_{\gamma}$, we try applying the methods discussed in the sections above to a
six-Majorana model.

Figure 7: 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
$n_{{nt1}_{1}}$ to $n_{{nt1}_{2}}$ for non-topological segment 1, and
$n_{{nt2}_{1}}$ to $n_{{nt2}_{2}}$ for non-topological segment 2. There are
six MZM’s (yellow circles) on this wire - three pairs localized at the edges
of the three topological segments.

###  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 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:

$\displaystyle H$ $\displaystyle=$ $\displaystyle
J_{123456}\prod_{i=1}^{6}\gamma_{i}+\sum_{1\leq i<j<k<l\leq
6}J_{ijkl}\gamma_{i}\gamma_{j}\gamma_{k}\gamma_{l}$ (13) $\displaystyle+$
$\displaystyle i\sum_{1\leq i<j\leq 6}K_{ij}\gamma_{i}\gamma_{j},$

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}\times 2^{3}$ Hilbert space. The eigenstates are in
the form $\ket{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:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\epsilon_{0}+\lambda_{1}(2n_{1}-1)+\lambda_{2}(2n_{2}-1)+\lambda_{3}(2n_{3}-1)$
(14) $\displaystyle+$ $\displaystyle
u_{12}(2n_{1}-1)(2n_{2}-1)+u_{13}(2n_{1}-1)(2n_{3}-1)$ $\displaystyle+$
$\displaystyle u_{23}(2n_{2}-1)(2n_{3}-1)$ $\displaystyle+$ $\displaystyle
v(2n_{1}-1)(2n_{2}-1)(2n_{3}-1),$

where $\lambda$’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
$\ket{n_{1}n_{2}n_{3}}$ are related to interaction strengths $\lambda$’s,
$u$’s and $v$ as:

$E=AI,$ (15)

where $E$ is the column matrix of energy levels for states
$\ket{n_{1}n_{2}n_{3}}$, where $n_{1},n_{2},n_{3}=0,1$ in the sequence as
shown below:

$E=\begin{pmatrix}E_{\ket{000}}\\\ E_{\ket{001}}\\\ E_{\ket{010}}\\\
E_{\ket{100}}\\\ E_{\ket{011}}\\\ E_{\ket{101}}\\\ E_{\ket{110}}\\\
E_{\ket{111}}\end{pmatrix},$ (16)

where $I$ is the column matrix of interaction strengths $\lambda$’s
(bilinear), $u$’s (quartic) and $v$ (sextic), i.e.

$I=\begin{pmatrix}\epsilon_{0}\\\ \lambda_{1}\\\ \lambda_{2}\\\ \lambda_{3}\\\
u_{12}\\\ u_{13}\\\ u_{23}\\\ v\end{pmatrix},$ (17)

and A is the matrix transforming $E$ to $I$ :

$A=\begin{pmatrix}1&-1&-1&-1&1&1&1&-1\\\ 1&-1&-1&1&1&-1&-1&1\\\
1&-1&1&-1&-1&1&-1&1\\\ 1&1&-1&-1&-1&-1&1&1\\\ 1&-1&1&1&-1&-1&1&-1\\\
1&1&-1&1&-1&1&-1&-1\\\ 1&1&1&-1&1&-1&-1&-1\\\ 1&1&1&1&1&1&1&1\\\
\end{pmatrix}.$ (18)

Thus we can extract the interaction strengths by inverting this relation to
solve for $I$:

$I=A^{-1}E.$ (19)

###  VI.D Kitaev chain model with interactions

(a)

(b)

Figure 8: The SYK point parameters for three-complex-fermion model for a
$N=10$ site Kitaev chain model (first non-topological segment:
$n_{\text{nt1}_{1}}=3$, $n_{\text{nt1}_{2}}=4$, second non-topological
segment: $n_{\text{nt2}_{1}}=7$, $n_{\text{nt2}_{2}}=8$) obtained using the
approximate method for SYK point search, i.e., method (i) described in Section
VI.E are: $t=0.0144$, $\Delta=0.0104$, $\mu_{\text{t1}}=0.0100$,
$\mu_{\text{t2}}=0.0100$, $\mu_{\text{t3}}=0.0100$, $\mu_{\text{nt1}}=3.0001$,
$\mu_{\text{nt2}}=3.0001$, $U=0.0039$. (a), At the SYK point (shown by black
dotted line-cut) $\lambda_{1}=5.5203e^{-06},$ $\lambda_{2}=1.1672e^{-06}$,
$\lambda_{3}=1.2026e^{-06}$, $u_{12}=u_{13}=u_{23}=5.3289e^{-04}$ and
$v=8.7008e^{-05}$ (b), At the SYK point, there is degeneracy between one even
and one odd parity level. 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.

(a)

(b)

Figure 9: The SYK point parameters for three-complex-fermion model for an
$N=10$ site Kitaev chain model (first non-topological segment:
$n_{\text{nt1}_{1}}=3$, $n_{\text{nt1}_{2}}=4$, second non-topological
segment: $n_{\text{nt2}_{1}}=7$, $n_{\text{nt2}_{2}}=8$) obtained using a
hybrid genetic algorithm search, i.e., method (ii) described in Section VI.E
are: $t=0.1$, $\Delta=0.1315$, $\mu_{\text{t1}}=0.5746$,
$\mu_{\text{t2}}=0.5743$, $\mu_{\text{t3}}=0.5746$, $\mu_{\text{nt1}}=10$,
$\mu_{\text{nt2}}=10$, $U=0.097$. Fig(a). At the SYK point at $U=0.097$ (shown
by black dotted line-cut) $\lambda_{1},=-1.4910e^{-04},$
$\lambda_{2}=-1.7743e^{-04}$, $\lambda_{3}=-2.5723e^{-04}$,
$u_{12}=u_{13}=u_{23}=-0.0013$ and $v=-1.234e^{-04}$. Fig(b) At the SYK point,
there is degeneracy between one even and one odd parity level. Focusing on the
line-cut 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.

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 (topo-nontopo-topo-nontopo-topo).
Following the conditions for topological phase we have
$|\mu_{\text{t1}}|<2|t|$, $|\mu_{\text{t2}}|<2|t|$ and
$|\mu_{\text{t3}}|<2|t|$ for the topological segments t1, t2, t3, and
$|\mu_{\text{nt1}}|>2|t|$, $|\mu_{\text{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. 7.

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:

$H_{\text{nl-int}}=U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{j}^{\dagger}c_{j}.$ (20)

The total Hamiltonian after adding up the interaction term is:

$H_{\text{tot}}=H_{\text{Kitaev-chain}}+H_{\text{nl-int}}.$ (21)

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
$\left|E_{\ket{111}}-E_{\ket{000}}\right|$ in the parameter space of
$\\{\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{t3}}$,
$\mu_{\text{nt1}}$, $\mu_{\text{nt2}}$, $t$, $\Delta\\}$. 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
$|\lambda_{1}|,|\lambda_{2}|,|\lambda_{3}|,|v|\simeq 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 $\mu_{\text{t1}}$, $\mu_{\text{t2}}$,
$\mu_{\text{t3}}$, $\mu_{\text{nt1}}$, $\mu_{\text{nt2}}$, $t$, $\Delta$ and
$U$. We use a hybrid Genetic Algorithm search to optimize for
$|\lambda_{1}|,|\lambda_{2}|,|\lambda_{3}|,v\simeq 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 $\left|E_{\ket{111}}-E_{\ket{000}}\right|\simeq 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. 8 where all interactions except for the quartic
interactions are well suppressed. In Fig. 8(a), following the black dotted
line-cut, at the SYK point we have
$\left|\frac{u}{\overline{\lambda}}\right|\simeq 200$ ($\overline{\lambda}$:
average of $\lambda_{1,2,3}$, $u=u_{12}=u_{13}=u_{13}$) and
$\left|\frac{u}{v}\right|\simeq 6$.

The energy levels used to compute the parameters $\lambda_{1,2,3}$,
$u_{12,23,13}$, and $v$ in Fig. 8(a) are plotted as a function of $U$ in Fig.
8(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. 9 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, $\left|\frac{u}{\overline{\lambda}}\right|\simeq
2$ ($\overline{\lambda}$: average of $\lambda_{1,2,3}$,
$u=u_{12}=u_{13}=u_{13}$) and $\left|\frac{u}{v}\right|\simeq 10.5$. For
completeness, we plot the energy levels that we used to extract
$\lambda_{1,2,3}$, $u_{12,23,13}$, and $v$ (in Fig. 9(a)) in Fig. 9(b). The
level structure is analogous to the one previously found in Fig. 8(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
$\lambda$’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_{\gamma}$ (particularly $N_{\gamma}=4,6$) SYK model in a 1D nanowire
system. It might be possible to extend this model to a $N_{\gamma}>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 Further Reading

For further background on topological states of matter and Majorana zero modes
we recommend these papers. [22, 23, 24, 25, 25, 26, 27, 28]. Variants of the
SYK model relevant to condensed matter physics: [29, 30, 31].

Discussions of the SYK Theory can be found here [3, 4, 1, 2, 32, 33, 34, 35,
36].

Works on the OTOC and Lyapunov Exponent: [5, 6].

Theories that include interacting MZMs: [37, 38, 19, 39, 40, 18, 41, 42].

Proposed experimental models and simulations of SYK model: [8, 9, 10, 11, 12,
14, 15, 13, 16, 43, 17].

Quantum transport studies related to the SYK model: [44, 45, 46, 47].

##  IX Study Design

This study was undertaken in a period of about nine months. A summary of our
study design is shown in Fig. 10.

Figure 10: Study design diagram

##  X Data Availability

All codes are available on Zenodo [48].

##  XI Acknowledgements

The authors thank J. Stenger for discussions of the SYK model.

##  XII Funding

S.F. and D.P. are supported by NSF PIRE-1743717. S.F. is supported by NSF
DMR-1906325, ONR and ARO.

##  Appendix A Constructing gamma matrices

The matrices representing the $\gamma$ operators are constructed by
transforming the Hamiltonian into spin chain basis (similar to Jordan Wigner
transformation) using the following algorithm.

Based on the Clifford algebra,

$\\{\gamma_{\mu},\gamma_{\nu}\\}=\delta_{\mu\nu}I_{2^{{N_{\gamma}}/2}\times
2^{{N_{\gamma}}/2}};\qquad\mu,\nu=1....N_{\gamma}$ (22)

For any ${N_{\gamma}}$, we can build the matrix representation of
$\gamma_{\mu}$ by taking product of Pauli matrices. Pauli matrices
($\sigma_{i}$’s) satisfy the Clifford algebra and hence forms the
representation:

$\\{\sigma_{i},\sigma_{j}\\}=\delta_{ij}I_{2\times 2}.$ (23)

So starting with ${N_{\gamma}}=2$, we can take $\gamma_{1}=\sigma_{2}$, and
$\gamma_{2}=\sigma_{1}$. We can use an iterative approach to obtain the
matrices for a higher dimension ${N_{\gamma}}$. Let
$\gamma_{\mu}(\mu=1....{N_{\gamma}}-2)$ be a $2^{{N_{\gamma}}/2-1}\times
2^{{N_{\gamma}}/2-1}$ matrices. Then $2^{{N_{\gamma}}/2}\times
2^{{N_{\gamma}}/2}$ matrices
$\gamma_{\widetilde{\mu}}(\widetilde{\mu}=1...{N_{\gamma}})$ are given as:

$\displaystyle\gamma_{\widetilde{\mu}}=\gamma_{\mu}\otimes-\sigma_{3},\qquad\text{for}\enskip\widetilde{\mu}=1,...,{N_{\gamma}}-2$
(24) $\displaystyle\gamma_{{N_{\gamma}}-1}=I_{2^{{N_{\gamma}}/2}\times
2^{{N_{\gamma}}/2}}\otimes\sigma_{2},\enskip\gamma_{{N_{\gamma}}}=I_{2^{{N_{\gamma}}/2}\times
2^{{N_{\gamma}}/2}}\otimes\sigma_{1}$ (25)

##  Appendix B Connecting the four MZM and the two complex fermion
representations

From the two-complex-fermion model [section II], a Hamiltonian with
interacting four MZM can be represented as

$H=H_{\text{bl}}+H_{1234}=i\sum_{1\leq i<j\leq
4}K_{ij}\gamma_{i}\gamma_{j}+J_{1234}\prod_{i=1}^{4}\gamma_{i}.$ (26)

Here $\gamma$’s represent the four MZM’s, $K_{ij}$’s are the bilinear
interaction strengths, where $K_{ij}=-K_{ji}$, and $J_{1234}$ is the quartic
interaction strength.

In the matrix form, this can be written as

$H=A^{\dagger}(iM)A+J_{1234}\prod_{i=1}^{4}\gamma_{i};\qquad
A=\begin{pmatrix}\gamma_{1}\\\ \gamma_{2}\\\ \gamma_{3}\\\ \gamma_{4}\\\
\end{pmatrix},$ (27)

where $M$ is a $4\times 4$ skew symmetric matrix with six unique non-zero
entries, i.e, $K_{ij}$’s. Eigenvalues of the skew-symmetric matrices come in
pairs: $\pm\lambda_{1},\pm\lambda_{2}$ with complex raising and lowering
operators $f$ and $f^{\dagger}$ as their eigenstates.

The eigenvectors of $M$ define the basis transformation $U$ from the MZM to
the complex fermion representation, i.e, from $\gamma$’s to $f$,
$f^{\dagger}$. The Hamiltonian corresponding to the bilinear MZM interactions,
$H_{\text{bl}}$ transforms as

$H_{\text{bl}}=B^{\dagger}\tilde{M}B;\qquad B=UA=\begin{pmatrix}f_{1}\\\
f_{1}^{\dagger}\\\ f_{2}\\\ f_{2}^{\dagger}\\\ \end{pmatrix},$ (28)

where $U$ is an unitary basis transformation matrix from the MZM basis to the
complex fermionic basis and $\tilde{M}$ is the transformed bilinear
interaction matrix

$\tilde{M}=U(iM)U^{-1}=\begin{pmatrix}\lambda_{1}&0&0&0\\\
0&-\lambda_{1}&0&0\\\ 0&0&\lambda_{2}&0\\\ 0&0&0&-\lambda_{2}\end{pmatrix}.$
(29)

The two quasiparticle energies are

$\displaystyle\lambda_{1,2}=\frac{1}{2}\sqrt{\frac{s\pm\sqrt{s^{2}-4d^{2}}}{2}},$
(30)

where $s=K_{12}^{2}+K_{13}^{2}+K_{14}^{2}+K_{23}^{2}+K_{24}^{2}+K_{34}^{2}$
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

$\displaystyle H_{\text{bl}}$ $\displaystyle=$ $\displaystyle
B^{\dagger}\tilde{M}B$ (31) $\displaystyle=$
$\displaystyle\lambda_{1}(f_{1}^{\dagger}f_{1}-f_{1}f_{1}^{\dagger})+\lambda_{2}(f_{2}^{\dagger}f_{2}-f_{2}f_{2}^{\dagger})$
$\displaystyle=$ $\displaystyle\lambda_{1}(2n_{1}-1)+\lambda_{2}(2n_{2}-1),$

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

$H_{1234}=u\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}=-u(f_{1}^{\dagger}f_{1}-f_{1}f_{1}^{\dagger})(f_{2}^{\dagger}f_{2}-f_{2}f_{2}^{\dagger})$
(32)

The total Hamiltonian in the complex fermion representation is therefore

$\displaystyle H$ $\displaystyle=$
$\displaystyle\lambda_{1}(2n_{1}-1)+\lambda_{2}(2n_{2}-1)$ (33)
$\displaystyle-$ $\displaystyle u(2n_{1}-1)(2n_{2}-1),$

where $\\{f_{i},f^{\dagger}_{j}\\}=2\delta_{ij}$ and
$f^{\dagger}_{i}f_{i}=2n_{i}$.

Thus, we can write the energy levels of this Hamiltonian corresponding to
states $\ket{00},\ket{10},\ket{01},\ket{11}$ in the form $\ket{n_{1}n_{2}}$ as

$\displaystyle E^{e}_{1}$
$\displaystyle=\epsilon_{0}-\lambda_{1}-\lambda_{2}-u,$ (34) $\displaystyle
E^{o}_{1}$ $\displaystyle=\epsilon_{0}+\lambda_{1}-\lambda_{2}+u,$ (35)
$\displaystyle E^{e}_{2}$
$\displaystyle=\epsilon_{0}+\lambda_{1}+\lambda_{2}-u,$ (36) $\displaystyle
E^{o}_{2}$ $\displaystyle=\epsilon_{0}-\lambda_{1}+\lambda_{2}+u$ (37)

where $E^{e}$’s are the eigenvalues of even parity states, i.e, $\ket{00}$ and
$\ket{11}$ and $E^{o}$’s are eigenvalues of the odd parity states, i.e.,
$\ket{10}$ and $\ket{01}$. $\epsilon_{0}$ is a constant shift in the energies.

If we know the spectrum of the eigenstates, we can obtain the Hamiltonian
parameters in the complex fermion representation, $\lambda_{1,2}$ and $u$,
using the linear transformation

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle(-E^{e}_{1}+E^{e}_{2}+E^{o}_{1}-E^{o}_{2})/4,$ (38)
$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle(-E^{e}_{1}+E^{e}_{2}-E^{o}_{1}+E^{o}_{2})/4,$ (39)
$\displaystyle u$ $\displaystyle=$
$\displaystyle(E^{o}_{1}+E^{o}_{2}-E^{e}_{1}-E^{e}_{2})/4.$ (40)

We note that multiple Hamiltonians in the MZM representation result in the
same eigenspectrum and hence connect to the same complex representation
because multiple sets of bilinear interactions $K_{ij}$ map onto the same
$\lambda_{1,2}$. However, there are two important features of the mapping
betwen the representations: (1) the quartic interaction $u$ is identical in
both representations (same in magnitude except for a sign change); (2) we are
interested in having all bilinear interactions in the MZM representation being
zero; Since $\lambda_{1}=\lambda_{2}=0$ if and only if
$K_{ij}=0\,\forall\,i,j$ we can verify that we have nulled the bilinear
interactions by verifying that $\lambda_{1}=\lambda_{2}=0$.

##  Appendix C Majorana mode wave function overlap optimization

Figure 11: $E_{\ket{11}}-E_{\ket{00}}$ is a function of parameters
$\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{nt}}$, $t$ and $\Delta$.
Colormap showing dependence of $E_{\ket{11}}-E_{\ket{00}}$ on a pair of
parameters in each sub-figure (keeping others fixed). The fixed parameters in
the sub-figures for $N=10$ site Kitaev chain with non-topological segment:
$n_{\text{nt}_{1}}=4;n_{\text{nt}_{2}}=7$ are as follows. (a), $t=1$,
$\Delta=0.5$, $\mu_{\text{nt}}=5$. (b), $\mu_{\text{t2}}=0.1$, $t=1$,
$\Delta=0.5$. (c), $\mu_{\text{t2}}=0.1$, $\mu_{\text{nt}}=5$, $\Delta=0.5$.
(d), $\mu_{\text{t2}}=0.1$, $\mu_{\text{nt}}=5$, $\Delta=0.5$. (e),
$\mu_{\text{t1}}=0.1$, $t=1$, $\Delta=0.5$. (f), $\mu_{\text{t1}}=0.1$,
$\mu_{\text{nt}}=5$, $\Delta=0.5$. (g), $\mu_{\text{t1}}=0.1$,
$\mu_{\text{nt}}=5$, $t=1$. (h), $\mu_{\text{t1}}=0.1$, $\mu_{\text{t2}}=0.1$,
$\Delta=0.5$. (i), $\mu_{\text{t1}}=0.1$, $\mu_{\text{t2}}=0.1$, $t=1$. (j),
$\mu_{\text{t1}}=0.1$, $\mu_{\text{t2}}=0.1$, $\mu_{\text{nt}}=5$.

Following discussion in Section V.A, we look for optimal points in the
parameter space of $\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{nt}}$,
$t$ and $\Delta$, at which the overlap integral of Majorana modes cancel out,
i.e., points at which $E_{\ket{11}}-E_{\ket{00}}\simeq 0$. We perform a global
search within the total parameter space to look for the global minima of
$\left|E_{\ket{11}}-E_{\ket{00}}\right|$. To visualize this optimization
process better, we have plotted $E_{\ket{11}}-E_{\ket{00}}$ as a function of
two parameters, keeping others fixed as shown in Fig. 11. This shows how each
parameter contributes in the optimization process.

##  Appendix D Spatial distribution of the MZM’s in the Kitaev model

(a)

(b)

Figure 12:  Left and right MZM’s ($\tilde{\gamma_{x}}$ (blue) and
$\tilde{\gamma_{y}}$ (yellow)) for states $\psi_{1}=\ket{01}$ (in (a)) and
$\psi_{1}=\ket{10}$ (in (b)). $\psi_{0}=\ket{00}$. The parameters are the
optimal point at which $\left|E_{\ket{11}}-E_{\ket{00}}\right|\simeq 0$:
$t=0.4025$, $\Delta=0.2167$, $\mu_{\text{t1}}=0.4832$,
$\mu_{\text{t2}}=0.4832$, $\mu_{\text{nt}}=8.5364$ for $N=10$ site Kitaev
chain with non-topological segment: $n_{\text{nt}_{1}}=4,n_{\text{nt}_{2}}=7$.

Following the discussion in Section V.A, overlap integral of MZM wave
functions in a 1D Kitaev chain (with no interactions) can be tuned to zero at
a point where lowest four many-body states are almost degenerate. Here we show
that indeed such MZM wave fucntion overlap can be seen at an optimal point
(which we find by following method in Section V.A.1) as shown in Fig. 12. To
see how these Majorana modes extend throughout the 1D chain, we find weights
of the left and right polarized Majorana modes on $\gamma$’s per site.

Defining Majorana creation and annihilation operators corresponding to the
left and right Majorana modes - $\tilde{\gamma_{x}}$ and $\tilde{\gamma_{y}}$
as:

$\tilde{\gamma_{x}}\ket{\psi_{0}}=\ket{\psi_{1}}\qquad\tilde{\gamma_{y}}\ket{\psi_{0}}=i\ket{\psi_{1}},$
(41)

where $\psi_{0}$ and $\psi_{1}$ are the ground states. We can find the weights
of $\tilde{\gamma_{x}}$ and $\tilde{\gamma_{y}}$ on each site as:

$\tilde{\gamma_{x}}=\sum_{j}{\alpha^{x}_{j}}\gamma_{j}\qquad\tilde{\gamma_{y}}=i\sum_{j}{\alpha^{y}_{j}}\gamma_{j},$
(42)

where $\alpha^{x}_{j}$ and $\alpha^{y}_{j}$ are the weights of
$\tilde{\gamma_{x}}$ and $\tilde{\gamma_{y}}$ operators on the $\gamma$’s per
site:

$\alpha^{x}_{j}=\bra{\psi_{1}}\gamma_{j}\ket{\psi_{0}}\qquad\alpha^{y}_{j}=\bra{\psi_{1}}\gamma_{j}\ket{\psi_{0}}.$
(43)

At the optimal point for $U=0$, $\tilde{\gamma_{x}}$ and $\tilde{\gamma_{y}}$
for $\psi_{0}\equiv\psi_{00}$ and $\psi_{1}\equiv\psi_{01}$ is shown in Fig.
12(a), and for $\psi_{0}\equiv\psi_{00}$, $\psi_{1}\equiv\psi_{10}$ is shown
in Fig. 12(b).

##  Appendix E  Finding optimal point using a global search algorithm

Following Section V.A.1, we search for global minima of the function
$E_{\ket{11}}-E_{\ket{00}}$ depending on the parameters $\mu_{\text{t1}}$,
$\mu_{\text{t2}}$, $\mu_{\text{nt}}$, $t$ and $\Delta$, i.e,
$\left|E_{\ket{11}}-E_{\ket{00}}\right|\simeq 0$. We use a MATLAB global
search algorithm [49] which performs multiple parallel searches (using a
nonlinear programming solver - ’fmincon’ [50]) through the parameter space
with different start points to find multiple local minima and then finalizes
at a global minimum. The search ranges that were used are:
$\\{\mu_{\text{t1}},\mu_{\text{t2}},t,\Delta\\}\in[0.1,1]$ and
$\mu_{\text{nt}}\in[3,10]$. Optimization stopping criteria: Function tolerance
$\sim e^{-10}$, Max Iterations: $3000$.

###  E.1 Three complex fermion, i.e, six MZM case

We search for global minima of the function
$\left|E_{\ket{111}}-E_{\ket{000}}\right|$ dependent on the parameters
$\mu_{\text{t1}}$, $\mu_{\text{t2}}$, $\mu_{\text{t3}}$, $\mu_{\text{nt1}}$,
$\mu_{\text{nt2}}$, $t$ and $\Delta$. The search algorithm is same as the two
complex fermion case, i.e., MATLAB global search algorithm [49]. The search
ranges that were used are:
$\\{\mu_{\text{t1}},\mu_{\text{t2}},\mu_{\text{t3}},t,\Delta\\}\in[0.1,1]$,
$\\{\mu_{\text{nt1}},\mu_{\text{nt2}}\\}\in[3,10]$.

##  Appendix F  Genetic Algorithm search

It starts by generating a random population of individuals (vectors of double
type). The next generation of population is selected based on elites
(individuals with best fitness to selectivity criteria), crossover (combining
certain parents to create children) and mutation (random modification to some
parents in the population). We use MATLAB’s Genetic Algorithm Toolbox [51] to
solve our problem. These are the criterion that we use in our code: Creation
function for creating initial population is the MATLAB default option.
Selection function that decides how to select next generation of population -
’selectiontournament’. Crossover function - ’crossoverscattered’. Mutation
Function - ’mutationadaptfeasible’. Hybrid function that refines search once
Genetic Algorithm search terminates - ’fmincon’ [50]. Optimization stopping
criteria: ’FunctionTolerance’ $\sim e^{-6}$, ’MaxStallGenerations’ (controls
the number of steps the Genetic Algorithm looks over to see whether it is
making progress) - $100$.

Parameters to be optimized are $\mu_{\text{t1}}$, $\mu_{\text{t2}}$,
$\mu_{\text{nt}}$, $t$, $\Delta$ and $U$. Range of search:
$\\{\mu_{\text{t1}},\mu_{\text{t2}},t,\Delta\\}\in[0.1,1]$,
$\mu_{\text{nt}}\in[3,10]$, $U\in[0.001,1]$.

###  F.1 Three complex fermion, i.e, six MZM case

Genetic Algorithm criterion used in the code are same as the two-complex-
fermion case. Parameters to be optimized are $\mu_{\text{t1}}$,
$\mu_{\text{t2}}$, $\mu_{\text{t3}}$, $\mu_{\text{nt1}}$, $\mu_{\text{nt2}}$,
$t$, $\Delta$ and $U$. Range of search:
$\\{\mu_{\text{t1}},\mu_{\text{t2}},\mu_{\text{t3}},t,\Delta\\}\in[0.1,1]$,
$\\{\mu_{\text{nt1}},\mu_{\text{nt2}}\\}\in[3,15]$, $U\in[0.01,1]$.

##  Appendix G Exploring other types of interactions

Figure 13:  Other interaction terms we checked are - (a)
$H_{\text{int}}=U\sum_{i}c_{i}^{\dagger}c_{i}c_{i+1}^{\dagger}c_{i+1}$, (b)
$H_{\text{int}}=U\sum_{i}c_{i}^{\dagger}c_{i}c_{i+2}^{\dagger}c_{i+2}$, (c)
$H_{\text{int}}=U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{i+3}^{\dagger}c_{i+3}$, (d)
$H_{\text{int}}=U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{i+4}^{\dagger}c_{i+4}$. All
other parameters are tuned to the optimal point at which
$\left|E_{\ket{11}}-E_{\ket{00}}\right|\simeq 0$: $t=0.4025$, $\Delta=0.2167$,
$\mu_{\text{t1}}=0.4832$, $\mu_{\text{t2}}=0.4832$, $\mu_{\text{nt}}=8.5364$
for $N=10$ site Kitaev chain with non-topological segment:
$n_{\text{nt}_{1}}=4,n_{\text{nt}_{2}}=7$.

Besides the non-linear interaction term $H_{\text{nl-
int}}=U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{j}^{\dagger}c_{j}$, we have explored
other interaction terms like
$H_{\text{int}}=U\sum_{i}c_{i}^{\dagger}c_{i}c_{i+1}^{\dagger}c_{i+1}$,
$U\sum_{i}c_{i}^{\dagger}c_{i}c_{i+2}^{\dagger}c_{i+2}$,
$U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{i+3}^{\dagger}c_{i+3}$ and
$U\sum_{i<j}c_{i}^{\dagger}c_{i}c_{i+4}^{\dagger}c_{i+4}$ (with interaction
strength $U$) shown in Fig. 13. We sweep across the value of $U$, fixing other
parameters at the optimal point (described in Section V.A). As we follow the
trajectory of $\lambda$’s and $u$ in the sub-figures, we do not find points
where $\lambda_{1,2}\simeq 0$ and $u$ is non-zero, thus showing that no SYK
point exists within the plausible range of $U$.

##  Appendix H Lambda’s follow different path beyond the optimal points

(a)

(b)

Figure 14:  Changing parameters from the optimal point shows that
$\lambda^{\prime}s$ no longer overlap and follow separate trajctories (solid
yellow and dashed blue line). Parameters used in this figure are as follows.
(a), $N=10$ site Kitaev chain with non-topological segment:
$n_{\text{nt}_{1}}=4;n_{\text{nt}_{2}}=7$; $t=0.4049$, $\Delta=0.2207$,
$\mu_{\text{t1}}=0.5513$, $\mu_{\text{t2}}=0.347$, $\mu_{\text{nt}}=7.7363$.
(b), $N=10$ fermion Kitaev chain with non-topological segment:
$n_{\text{nt}_{1}}=4;n_{\text{nt}_{2}}=7$; $t=0.2801$, $\Delta=0.1278$,
$\mu_{\text{t1}}=0.3544$, $\mu_{\text{t2}}=0.4324$, $\mu_{\text{nt}}=8.515$.

As shown in Figs. 4, 5, 6, $\lambda_{1}=\lambda_{2}$ as we tune the non-linear
interaction strength $U$. This is a consequence of starting at either an
optimal point or an SYK point. In Fig. 14, we show that if we change the
parameters slightly, then the $\lambda$’s follow separate paths.

## References

  * Kitaev [2015] A. Kitaev, A simple model of quantum holography, KITP Strings Seminar and Entanglement 2015 Program (2015).
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Gapless spin-fluid ground state in a random quantum Heisenberg magnet, Phys. Rev. Lett. 70, 3339 (1993).
  * Polchinski and Rosenhaus [2016] J. Polchinski and V. Rosenhaus, The spectrum in the Sachdev-Ye-Kitaev model, Journal of High Energy Physics 10.1007/JHEP04(2016)001 (2016).
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the Sachdev-Ye-Kitaev model, Phys. Rev. D 94, 106002 (2016).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet, Nature Physics 10.1038/nphys4119 (2017).
  * Maldacena  _et al._ [2006] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 10.1007/JHEP08(2016)106 (2006).
  * Sachdev [2015] S. Sachdev, Bekenstein-Hawking entropy and strange metals, Phys. Rev. X 5, 041025 (2015).
  * Chew  _et al._ [2017] A. Chew, A. Essin, and J. Alicea, Approximating the Sachdev-Ye-Kitaev model with Majorana wires, Phys. Rev. B 96, 121119 (2017).
  * Haenel  _et al._ [2021] R. Haenel, S. Sahoo, T. H. Hsieh, and M. Franz, Traversable wormhole in coupled Sachdev-Ye-Kitaev models with imbalanced interactions, Phys. Rev. B 104, 035141 (2021).
  * Pikulin and Franz [2017] D. I. Pikulin and M. Franz, Black hole on a chip: Proposal for a physical realization of the Sachdev-Ye-Kitaev model in a solid-state system, Phys. Rev. X 7, 031006 (2017).
  * Chen  _et al._ [2018] A. Chen, R. Ilan, F. de Juan, D. I. Pikulin, and M. Franz, Quantum holography in a graphene flake with an irregular boundary, Phys. Rev. Lett. 121, 036403 (2018).
  * Danshita  _et al._ [2017] I. Danshita, M. Hanada, and M. Tezuka, Creating and probing the Sachdev-Ye-Kitaev model with ultracold gases: Towards experimental studies of quantum gravity, Progress of Theoretical and Experimental Physics 2017, 10.1093/ptep/ptx108 (2017), 083I01, https://academic.oup.com/ptep/article-pdf/2017/8/083I01/19650704/ptx108.pdf .
  * Wei and Sedrakyan [2021] C. Wei and T. A. Sedrakyan, Optical lattice platform for the Sachdev-Ye-Kitaev model, Phys. Rev. A 103, 013323 (2021).
  * Luo  _et al._ [2019] Z. Luo, Y.-Z. You, J. Li, C.-M. Jian, D. Lu, C. Xu, B. Zeng, and R. Laflamme, Quantum simulation of the non-Fermi-liquid state of Sachdev-Ye-Kitaev model, npj Quantum Information 5, 10.1038/s41534-019-0166-7 (2019).
  * García-Álvarez  _et al._ [2017] L. García-Álvarez, I. L. Egusquiza, L. Lamata, A. del Campo, J. Sonner, and E. Solano, Digital quantum simulation of minimal $\mathrm{AdS}/\mathrm{CFT}$, Phys. Rev. Lett. 119, 040501 (2017).
  * Behrends and Béri [2020] J. Behrends and B. Béri, Sachdev-Ye-Kitaev circuits for braiding and charging Majorana zero modes, arXiv:2012.09210 [cond-mat, physics:hep-th, physics:quant-ph] (2020).
  * Babbush  _et al._ [2019] R. Babbush, D. W. Berry, and H. Neven, Quantum simulation of the Sachdev-Ye-Kitaev model by asymmetric qubitization, Phys. Rev. A 99, 040301 (2019).
  * Chiu  _et al._ [2015] C.-K. Chiu, D. I. Pikulin, and M. Franz, Strongly interacting Majorana fermions, Physical Review B 91, 165402 (2015).
  * Ghazaryan and Chakraborty [2015] A. Ghazaryan and T. Chakraborty, Long-range Coulomb interaction and Majorana fermions, Phys. Rev. B 92, 115138 (2015).
  * Domínguez  _et al._ [2017] F. Domínguez, J. Cayao, P. San-Jose, R. Aguado, A. L. Yeyati, and E. Prada, Zero-energy pinning from interactions in Majorana nanowires, npj Quantum Materials 2, 13 (2017).
  * [21] J. D’Errico, Eigenshuffle, https://www.mathworks.com/matlabcentral/fileexchange/22885.
  * Majorana [1938] E. Majorana, Teoria simmetrica dell’elettrone e del positrone, Il Nuovo Cimento (1924-1942) 10.1007/BF02961314 (1938).
  * Kitaev [2001] A. Kitaev, Unpaired Majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Signatures of Majorana fermions in hybrid superconductor-semiconductor nanowire devices, Science 336, 1003 (2012).
  * Leijnse and Flensberg [2012] M. Leijnse and K. Flensberg, Introduction to topological superconductivity and Majorana fermions, Semiconductor Science and Technology 27, 124003 (2012).
  * Wang and Zhang [2017] J. Wang and S.-C. Zhang, Topological states of condensed matter, Nature Materials 16, 10.1038/nmat5012 (2017).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Helical liquids and Majorana bound states in quantum wires, Phys. Rev. Lett. 105, 177002 (2010).
  * Lantagne-Hurtubise  _et al._ [2021] É. Lantagne-Hurtubise, V. Pathak, S. Sahoo, and M. Franz, Superconducting instabilities in a spinful Sachdev-Ye-Kitaev model, Physical Review B 104, L020509 (2021).
  * Fu  _et al._ [2017] W. Fu, D. Gaiotto, J. Maldacena, and S. Sachdev, Supersymmetric Sachdev-Ye-Kitaev models, Phys. Rev. D 95, 026009 (2017).
  * Can and Franz [2019] O. Can and M. Franz, Solvable model for quantum criticality between the Sachdev-Ye-Kitaev liquid and a disordered Fermi liquid, Physical Review B 100, 045124 (2019).
  * Rosenhaus [2019] V. Rosenhaus, An introduction to the SYK model, Journal of Physics A: Mathematical and Theoretical 52, 323001 (2019).
  * Li  _et al._ [2017] T. Li, J. Liu, Y. Xin, and Y. Zhou, Supersymmetric SYK model and random matrix theory, Journal of High Energy Physics 10.1007/JHEP06(2017)111 (2017).
  * Trunin [2021] D. A. Trunin, Pedagogical introduction to the Sachdev-Ye-Kitaev model and two-dimensional dilaton gravity, Physics-Uspekhi 64, 219 (2021).
  * García-García and Verbaarschot [2016] A. M. García-García and J. J. Verbaarschot, Spectral and thermodynamic properties of the Sachdev-Ye-Kitaev model, Physical Review D 94, 126010 (2016).
  * Sarosi [2018] G. Sarosi, AdS$_{2}$ holography and the SYK model, in  _Proceedings of XIII Modave Summer School in Mathematical Physics — PoS(Modave2017)_ (Sissa Medialab, Modave, Belgium, 2018) p. 001.
  * Manolescu  _et al._ [2014] A. Manolescu, D. C. Marinescu, and T. D. Stanescu, Coulomb interaction effects on the Majorana states in quantum wires, Journal of Physics: Condensed Matter 26, 172203 (2014).
  * Więckowski and Ptok [2019] A. Więckowski and A. Ptok, Influence of long-range interaction on Majorana zero modes, Phys. Rev. B 100, 144510 (2019).
  * Ng [2015] H. T. Ng, Decoherence of interacting Majorana modes, Scientific Reports 5, 12530 (2015).
  * Rahmani and Franz [2019] A. Rahmani and M. Franz, Interacting Majorana fermions, Reports on Progress in Physics 82, 084501 (2019).
  * Hassler and Schuricht [2012] F. Hassler and D. Schuricht, Strongly interacting Majorana modes in an array of Josephson junctions, New Journal of Physics 14, 125018 (2012).
  * Stoudenmire  _et al._ [2011] E. M. Stoudenmire, J. Alicea, O. A. Starykh, and M. P. Fisher, Interaction effects in topological superconducting wires supporting Majorana fermions, Phys. Rev. B 84, 014503 (2011).
  * Lantagne-Hurtubise  _et al._ [2018] É. Lantagne-Hurtubise, C. Li, and M. Franz, Family of Sachdev-Ye-Kitaev models motivated by experimental considerations, Phys. Rev. B 97, 235124 (2018).
  * Kruchkov  _et al._ [2020] A. Kruchkov, A. A. Patel, P. Kim, and S. Sachdev, Thermoelectric power of Sachdev-Ye-Kitaev islands: Probing Bekenstein-Hawking entropy in quantum matter experiments, Phys. Rev. B 101, 205148 (2020).
  * Gnezdilov  _et al._ [2018] N. V. Gnezdilov, J. A. Hutasoit, and C. W. J. Beenakker, Low-high voltage duality in tunneling spectroscopy of the Sachdev-Ye-Kitaev model, Phys. Rev. B 98, 081413 (2018).
  * Song  _et al._ [2017] X.-Y. Song, C.-M. Jian, and L. Balents, Strongly correlated metal built from Sachdev-Ye-Kitaev models, Phys. Rev. Lett. 119, 216601 (2017).
  * Altland  _et al._ [2019] A. Altland, D. Bagrets, and A. Kamenev, Sachdev-Ye-Kitaev non-Fermi-liquid correlations in nanoscopic quantum transport, Phys. Rev. Lett. 123, 226801 (2019).
  * [48] DOI: 10.5281/zenodo.5555501, https://zenodo.org/record/5555501.
  * The MathWorks [020aa] I. The MathWorks, Global or multiple starting point search, global optimization toolbox version 4.3, https://www.mathworks.com/help/gads/global-or-multiple-starting-point-search.html (MATLAB 2020aa).
  * The MathWorks [020ab] I. The MathWorks, fmincon, optimization toolbox version 8.5, https://www.mathworks.com/help/optim/ug/fmincon.html (MATLAB 2020ab).
  * The MathWorks [020ac] I. The MathWorks, Genetic algorithm, global optimization toolbox version 4.3, https://www.mathworks.com/help/gads/genetic-algorithm.html (MATLAB 2020ac).