Black hole on a chip: proposal for a physical realization of the SYK model in a solid-state system

TL;DR

This paper proposes a solid-state experimental setup using Majorana zero modes in a topological insulator-superconductor interface to realize the SYK model, connecting quantum chaos, holography, and black hole physics.

Contribution

It introduces a feasible physical realization of the SYK model in a solid-state system using Majorana modes in a topological insulator-superconductor heterostructure.

Findings

Numerical simulations confirm SYK-like thermodynamic behavior.

The setup exhibits characteristic two-point and four-point correlators.

Potential experimental methods for observing SYK physics are discussed.

Abstract

System of Majorana zero modes with random infinite range interactions -- the Sachdev-Ye-Kitaev (SYK) model -- is thought to exhibit an intriguing relation to the horizons of extremal black holes in two-dimensional anti-de Sitter (AdS$_2$) space. This connection provides a rare example of holographic duality between a solvable quantum-mechanical model and dilaton gravity. Here we propose a physical realization of the SYK model in a solid state system. The proposed setup employs the Fu-Kane superconductor realized at the interface between a three dimensional topological insulator (TI) and an ordinary superconductor. The requisite $N$ Majorana zero modes are bound to a nanoscale hole fabricated in the superconductor that is threaded by $N$ quanta of magnetic flux. We show that when the system is tuned to the surface neutrality point (i.e. chemical potential coincident with the Dirac point of the TI surface state) and the hole has sufficiently irregular shape, the Majorana zero modes are described by the SYK Hamiltonian. We perform extensive numerical simulations to demonstrate that the system indeed exhibits physical properties expected of the SYK model, including thermodynamic quantities and two-point as well as four-point correlators, and discuss ways in which these can be observed experimentally.

Tables

Table 1. Table 1: Gaussian ensembles for even N 𝑁 N .

$N\left(mod8\right)$	0	2	4	6
level stat.	GOE	GUE	GSE	GUE
$\beta$	1	2	4	2
$Z$	$\frac{8}{27}$	$\frac{4\pi }{81\sqrt{3}}$	$\frac{4\pi }{729\sqrt{3}}$	$\frac{4\pi }{81\sqrt{3}}$

Equations

$${\mathcal{H}}_{\mathrm{SYK}}=\sum _{i<j<k<l}{J}_{ijkl}{\chi }_i{\chi }_j{\chi }_k{\chi }_l,$$

$$\{{\chi }_i,{\chi }_j\}={\delta }_{ij},{\chi }_j^{†}={\chi }_j.$$

$${\mathcal{H}}_{\mathrm{FK}}=\int d^{2}r{\hat{\Psi }}_{\mathbf{r}}^{†}{H}_{\mathrm{FK}}(\mathbf{r}){\hat{\Psi }}_{\mathbf{r}},$$

$$H_{\mathrm{FK}}={\tau }^z\left[v_F\mathbf{\sigma }\cdot \left(\mathbf{p}-{\tau }^z\frac{e}{c}\mathbf{A}\right)-\mu \right]+{\tau }^x{\Delta }_1-{\tau }^y{\Delta }_2.$$

$$\Delta({\bm{r}})=\left\{\begin{array}[]{ll}0,&r<R(\varphi)\\ \Delta_{0}e^{iN\varphi},&r>R(\varphi),\end{array}\right.$$

$$\mathcal{H}={\mathcal{H}}_{\mathrm{FK}}^{(N)}+\delta {\mathcal{H}}_{\mathrm{FK}}+{\mathcal{H}}_{\mathrm{int}}.$$

$${\mathcal{H}}_{\mathrm{int}}=\frac{1}{2}\int \int d^{2}r{d}^2{r}^{\prime }\hat{\rho }(\mathbf{r})V(\mathbf{r}-{\mathbf{r}}^{\prime })\hat{\rho }({\mathbf{r}}^{\prime }),$$

$${\mathcal{H}}_{\mathrm{FK}}^{(N)}={\sum _n}^{\prime }{\epsilon }_n{\hat{\psi }}_n^{†}{\hat{\psi }}_n+E_g$$

$${\hat{\psi }}_n=\int d^{2}r{\Phi }_n^{†}(\mathbf{r}){\hat{\Psi }}_{\mathbf{r}}$$

$${\sigma }^y{\tau }^y{\Phi }_j^{\ast }(\mathbf{r})={\Phi }_j(\mathbf{r})$$

$${\hat{\Psi }}_{\mathbf{r}}=\sum _n{\Phi }_n(\mathbf{r}){\hat{\psi }}_n,$$

$${\mathcal{H}}_{\mathrm{eff}}=\frac{i}{2!}\sum _{i,j}{\tilde{K}}_{ij}{\chi }_i{\chi }_j+\frac{1}{4!}\sum _{i,j,k,l}{\tilde{J}}_{ijkl}{\chi }_i{\chi }_j{\chi }_k{\chi }_l,$$

$$i{\tilde{K}}_{ij}$$

$${\tilde{J}}_{ijkl}$$

$${\Phi }_j=\left(\begin{array}{c}{\eta }_j\\ i{\sigma }^y{\eta }_j^{\ast }\end{array}\right),$$

$${\rho }_{ij}=\frac{i}{2}({\eta }_i^{†}{\eta }_j-c.c.)=-\mathrm{Im}({\eta }_i^{†}{\eta }_j).$$

$${\mathcal{H}}_{\mathrm{eff}}=i\sum _{i<j}{K}_{ij}{\chi }_i{\chi }_j+\sum _{i<j<k<l}{J}_{ijkl}{\chi }_i{\chi }_j{\chi }_k{\chi }_l,$$

$$K_{ij}=\frac{1}{2}({\tilde{K}}_{ij}-{\tilde{K}}_{ji}),J_{ijkl}=\frac{1}{3}({\tilde{J}}_{ijkl}-{\tilde{J}}_{ikjl}+{\tilde{J}}_{lijk})$$

$$K^2=N\overline{K_{ij}^2},J^2=\frac{N^3}{3!}\overline{J_{ijkl}^2},$$

$${\bar{\eta }}_j({\mathbf{r}}_n)=\left(\begin{array}{c}{\varphi }_j^1(n)+i{\varphi }_j^2(n)\\ {\varphi }_j^3(n)+i{\varphi }_j^4(n)\end{array}\right),$$

$$\overline{{\varphi }_i^{\alpha }(n)}=0,\overline{{\varphi }_i^{\alpha }(n){\varphi }_j^{\beta }(n)}=\frac{1}{8M_s}{\delta }_{ij}{\delta }^{\alpha \beta }.$$

$$J_{ijkl}=-\frac{V_0}{{\zeta }^2}\sum _{n=1}^{M_s}{ϵ}_{\alpha \beta \mu \nu }{\varphi }_i^{\alpha }(n){\varphi }_j^{\beta }(n){\varphi }_k^{\mu }(n){\varphi }_l^{\nu }(n),$$

$$\overline{J_I}=0,\overline{J_I{J}_J}={\left(\frac{3V_0}{8{\zeta }^2}\right)}^2\frac{1}{M_s^3}{\delta }_{IJ},$$

$$G(\tau )=⟨{\mathcal{T}}_{\tau }\chi (\tau )\chi (0)⟩$$

$$G({\omega }_n)=[-i{\omega }_n-\Sigma ({\omega }_n){]}^{-1}.$$

$$\Sigma (\tau )=K^{2}G(\tau )+J^2{G}^3(\tau ).$$

$$G_f({\omega }_n)=\frac{2i}{{\omega }_n+\mathrm{sgn}({\omega }_n)\sqrt{{\omega }_n^2+4K^2}}.$$

$$A_f(\omega )=\frac{1}{\pi K}\mathrm{Re}\sqrt{1-{\left(\frac{\omega }{2K}\right)}^2},$$

$$G_c({\omega }_n)=i{\pi }^{1/4}\frac{\mathrm{sgn}({\omega }_n)}{\sqrt{J|{\omega }_n|}}$$

$$A_c(\omega )=\frac{1}{\sqrt{2}{\pi }^{3/4}}\frac{1}{\sqrt{J|\omega |}}.$$

$$G_f(\tau )=\frac{1}{\pi K}\frac{\mathrm{sgn}(\tau )}{|\tau |},G_c(\tau )=\frac{1}{{\pi }^{1/4}\sqrt{2J}}\frac{\mathrm{sgn}(\tau )}{|\tau {|}^{1/2}}.$$

Full text

I Introduction

Models of particles with infinite-range interactions have a long history in nuclear physics dating back to the pioneering works of Wigner and Dyson Wigner (1951); Dyson (1962) and in condensed matter physics in studies describing spin glass and spin liquid states of matter Sherrington and Kirkpatrick (1975); Sachdev and Ye (1993); Camjayi and Rozenberg (2003) . More recently, Kitaev Kitaev (2015); Maldacena and Stanford (2016) formulated and studied a Majorana fermion version of the model with all-to-all random interactions first proposed by Sachdev and Ye Sachdev and Ye (1993). The resulting SYK model, defined by the Hamiltonian (1) below, is solvable in the limit of large number $N$ of fermions and exhibits a host of intriguing properties. The SYK model is believed to be holographic dual of extremal black hole horizons in AdS2 and has been argued to possess remarkable connections to information theory, many-body thermalization and quantum chaos Sachdev (2015); Maldacena et al. (2016); Hosur et al. (2016); Polchinski and Rosenhaus (2016); García-García and Verbaarschot (2016); You et al. (2016). Various extensions of the SYK model have been put forth containing supersymmetry Fu et al. (2017), interesting quantum phase transition Banerjee and Altman (2017); Bi et al. (2017), higher dimensional extensions Gu et al. (2016); Berkooz et al. (2017), as well as a version that does not require randomness Witten (2016). Given its fascinating properties it would be of obvious interest to have an experimental realization of the SYK model or its variants. Thus far a realization of the Sachdev-Ye model (with complex fermions) has been proposed using ultracold gases Danshita et al. (2016) and a protocol for digital quantum simulation of both the complex and Majorana fermion versions of the model has been discussed García-Álvarez et al. (2016). A natural realization of the SYK model in a solid state system is thus far lacking.

Recent years have witnesed numerous proposals for experimental realizations of unpaired Majorana zero modes in solid state systems Alicea (2012); Beenakker (2013); Leijnse and Flensberg (2012); Stanescu and Tewari (2013); Elliott and Franz (2015), with compelling experimental evidence for their existence gradually mounting in several distinct platforms Mourik et al. (2012); Das et al. (2012); Deng et al. (2012); Rokhinson et al. (2012); Finck et al. (2013); Hart et al. (2014); Nadj-Perge et al. (2014); Xu et al. (2015); Sun et al. (2016). The purpose of this paper is to propose a physical realization of the SYK model in one of these platforms. The SYK Hamiltonian we wish to implement is given by

[TABLE]

where $J_{ijkl}$ are random independent coupling constants and $\chi_{j}$ represent the Majorana zero mode operators that obey the canonical anticommutation relations

[TABLE]

The proposed device, depicted in Fig. 1, employs an interface between a 3D TI and an ordinary superconductor such as Nb or Pb. Fu and Kane Fu and Kane (2008) showed theoretically that magnetic vortices in such an interface host unpaired Majorana zero modes and signatures consistent with this prediction have been reported in Bi2Te3/NbSe2 heterostructures Xu et al. (2015); Sun et al. (2016). Under ordinary circumstances these vortices tend to form an Abrikosov lattice and the low-energy effective theory is dominated by two-fermion terms $iK_{ij}\chi_{i}\chi_{j}$ with the hoping amplitudes $K_{ij}$ decaying exponentially with the distance between vortex sites $|{\bm{r}}_{i}-{\bm{r}}_{j}|$. Four-fermion interaction terms of the type required to implement the SYK Hamiltonian (1) are generically also present but are subdominant and also decay exponentially with distance. Realizing the SYK model in this setup therefore entails two key challenges: (i) one must find a way to suppress the two-fermion tunneling terms and (ii) render the four-fermion interactions effectively infinite-ranged. In addition the four-fermion coupling constants $J_{ijkl}$ must be sufficiently random. In the following we show how these challenges can be overcome by judicious engineering of various aspects of the device depicted in Fig. 1.

The first challenge can be met by tuning the surface state of the TI into its global neutrality point such that the chemical potential $\mu$ lies at the Dirac point. At the neutrality point the interface superconductor is known to acquire an extra chiral symmetry which prohibits any two-fermion terms Teo and Kane (2010). In other words, the symmetry requires $K_{ij}=0$ and the low-energy Hamiltonian is then dominated by the four-fermion terms Chiu et al. (2015). The second requirement of effectively infinite-ranged interactions can be satisfied by localizing all Majorana zero modes in the same region of space. In our proposed device this is achieved by fabricating a hole in the SC layer as illustrated in Fig. 1. If the sample is cooled in a weak applied magnetic field an integer number $N$ of magnetic flux quanta can be trapped in the hole. The SC phase $\theta$ will then wind by $2\pi N$ around the hole, forming effectively an $N$-fold vortex with $N$ Majorana zero modes bound to the hole. If furthermore the hole is designed to have an irregular shape the Majorana wavefunctions will have random spatial structure and their overlaps will give rise to the required randomness in the coupling constants $J_{ijkl}$. This randomness is related to random classical trajectories inside such a hole, or ‘billiards’ as it is commonly called in the quantum chaos literature Gutzwiller (1990); Haake (2013). We note a related proposal to realize the SYK model using semiconductor quantum wires coupled to a disordered quantum dot advanced in the recent work Chew et al. (2017).

In the rest of the paper we provide the necessary background on our proposed system and support its relation to the SYK model by physical arguments and by detailed model calculations. We first review the Fu-Kane model Fu and Kane (2008) for the TI/SC interface and numerically calculate the Majorana wavefunctions localized in a hole threaded by $N$ magnetic flux quanta in the presence of disorder. Assuming that the constituent electrons interact via screened Coulomb potential we then explicitly calculate the four-fermion coupling constants $J_{ijkl}$ between the Majorana zero modes. We finally use these as input data for the many-body Majorana Hamiltonian which we diagonalize numerically for $N$ up to 32 and study its thermodynamic properties, level statistics, as well as two- and four-point correlators. We show that these behave precisely as expected of the SYK model with random independent couplings. We also discuss the effect of small residual two-fermion terms that will inevitably be present in a realistic device and propose ways to experimentally detect signatures of the SYK physics using tunneling spectroscopy.

II SYK model from Interacting Majorana zero modes at the TI/SC interface

II.1 The Fu-Kane superconductor

The surface of a canonical 3D TI, such as Bi2Se3, hosts a single massless Dirac fermion protected by time reversal symmetry. When placed in the proximity of an ordinary superconductor the surface state is described by the Fu-Kane Hamiltonian Fu and Kane (2008)

[TABLE]

where $\hat{\Psi}_{\bm{r}}=(c_{\uparrow{\bm{r}}},c_{\downarrow{\bm{r}}},c^{\dagger}_{\downarrow{\bm{r}}},-c^{\dagger}_{\uparrow{\bm{r}}})^{T}$ is the Nambu spinor and

[TABLE]

Here $v_{F}$ is the velocity of the surface state, ${\bm{p}}=-i\hbar\nabla$ denotes the momentum operator, $\Delta=\Delta_{1}+i\Delta_{2}$ is the SC order parameter and ${\bm{\sigma}}$, ${\bm{\tau}}$ are Pauli matrices in spin and Nambu spaces, respectively. To describe the geometry depicted in Fig. 1 we take

[TABLE]

where $\varphi$ is the polar angle and $R(\varphi)$ denotes the hole radius. The vector potential is taken to yield total flux through the hole $\oint_{C}d{\bf l}\cdot{\bm{A}}=N\Phi_{0}$ with $\Phi_{0}=hc/2e$ the SC flux quantum and the contour $C$ taken to encircle the hole at a radius well beyond the effective magnetic penetration depth of the SC film $\lambda_{\rm eff}=2\lambda_{L}^{2}/d$. (Here $\lambda_{L}$ is the London penetration depth of the bulk SC and $d$ the film thickness.)

Hamiltonian (4) respects the particle-hole symmetry generated by $\Xi=\sigma^{y}\tau^{y}{\cal K}$ ($\Xi^{2}=+1$), where ${\cal K}$ denotes complex conjugation. For a purely real gap function $\Delta$ and zero magnetic field ${\bm{B}}=\nabla\times{\bm{A}}$, it also obeys the physical time reversal symmetry $\Theta=i\sigma^{y}{\cal K}$ ($\Theta^{2}=-1$). In the presence of vortices $\Delta$ becomes complex and the time reversal symmetry is broken. The Fu-Kane model with vortices therefore falls into symmetry class D in the Altland-Zirnbauer classification Altland and Zirnbauer (1997) which, in accordance with Ref. Teo and Kane (2010), implies a Z2 classification for the zero modes associated with vortices. Physically, this means that a system with total vorticity $N_{V}$ will have $N=(N_{V}\!\mod 2)$ exact zero modes, in accord with the expectation that any even number of Majorana zero modes will generically hybridize and form complex fermions at non-zero energy.

When $\mu=0$, Hamiltonian (4) also respects a fictitious time reversal symmetry generated by $\Sigma=\sigma^{x}\tau^{x}{\cal K}$ ($\Sigma^{2}=+1$).

It is important to note that unlike the physical time reversal this symmetry remains valid even in the presence of the applied magnetic field and vortices. At the neutrality point, the two symmetries $\Xi$ and $\Sigma$ define a BDI class with chiral symmetry $\Pi=\Xi\Sigma=\sigma^{z}\tau^{z}$. This, in accordance with Ref. Teo and Kane (2010), implies an integer classification of zero modes associated with point defects. A system with total vorticity $N_{V}$ will thus exhibit $N=N_{V}$ exact zero modes, irrespective of the precise geometric arrangement of the individual vortices and other details. This remains true in the presence of any disorder that does not break the $\Sigma$ symmetry. Specifically randomness in $v_{F}$ and $\Delta$ will not split the zero modes but random contributions to $\mu$ will.

Another way to establish the existence of exact zero modes in the Hamiltonian (4) with $\mu=0$ is to recognize it as a version of the Jackiw-Rossi Hamiltonian Jackiw and Rossi (1981) well known in particle physics. An index theorem for this Hamiltonian, conjectured by Jackiw and Rossi and later proven by Weinberg Weinberg (1981), equates the number $N$ of its exact zero modes in region ${\cal M}$ to the total vorticity $N_{V}={1\over 2\pi}\oint_{\partial{\cal M}}d{\bf l}\cdot\nabla\theta$ contained in that region. A region threaded by $N_{V}$ magnetic flux quanta will thus contain $N$ exact zero modes.

In the geometry of Fig. 1 the Majorana modes discussed above can equivalently be viewed as living at the boundary between a magnetically gapped TI surface on the inside and a SC region on the outside of the hole. The existence of such modes is well known and has been discussed in several papers Fu and Kane (2009); Akhmerov et al. (2009).

The existence and properties of the zero modes in the Fu-Kane Hamiltonian have been extensively tested by analytic and numerical approaches for a single vortex Fu and Kane (2008), pair of vortices Cheng et al. (2009, 2010), periodic Abrikosov lattices Liu and Franz (2015); Murray and Vafek (2015) as well as the “giant vortex” geometry Pikulin et al. (2015) similar to our proposed setup. This body of work firmly establishes the existence of exact Majorana zero modes for $\mu=0$ in accordance with the Jackiw-Rossi-Weinberg index theorem. Away from neutrality it is found that the zero modes are split due to two-fermion tunneling terms $K_{ij}\propto\mu$ where the constant of proportionality is related to the wavefunction overlap between $\chi_{i}$ and $\chi_{j}$. In addition, it has been found that for a singly quantized vortex at neutrality the zero mode is separated from the rest of the spectrum by a gap $\sim\Delta_{0}$ where $\Delta_{0}$ is the SC gap magnitude far from the vortex. We shall see that for a judiciously chosen hole size this convenient hierarchy of energy scales remains in place with $N$ zero modes separated by a large gap from the rest of the spectrum.

II.2 Low-energy effective theory

Having established a convenient platform that hosts $N$ Majorana zero modes with wavefunctions localized in the same region of space we now proceed to derive the effective low-energy theory in terms of the Majorana zero mode operators $\chi_{j}$. To this end we write the full second-quantized Hamiltonian of the system as

[TABLE]

Here ${\cal H}_{\rm FK}^{(N)}$ stands for the part of the Fu-Kane Hamiltonian (4) that obeys the fictitious time-reversal symmetry $\Sigma$ and exhibits therefore $N$ exact zero modes. $\delta{\cal H}_{\rm FK}$ contains all the remaining fermion bilinears that break $\Sigma$ such as the chemical potential term. ${\cal H}_{\rm int}$ defines the four-fermion interactions that have been ignored thus far but will play pivotal role in the physics of the SYK model we wish to study. We assume that electrons are subject to screened Coulomb interactions described by

[TABLE]

where $V({\bm{r}})$ is the interaction potential and $\hat{\rho}({\bm{r}})=c^{\dagger}_{\sigma{\bm{r}}}c_{\sigma{\bm{r}}}$ is the electron charge density operator.

Now imagine we have solved the single-electron problem defined by Hamiltonian ${H}_{\rm FK}^{(N)}$ for the device geometry sketched in Fig. 1 with $N$ flux quanta threaded through the hole. We thus have the complete set of single-particle eigenfunctions $\Phi_{n}({\bm{r}})$ and eigenenergies $\varepsilon_{n}$ of $H_{\rm FK}^{(N)}$. The corresponding second quantized Hamiltonian can then be written in a diagonal form

[TABLE]

where

[TABLE]

is the eigenmode operator belonging to the eigenvalue $\varepsilon_{n}$. The sum over $n$ is restricted to the positive energy eigenvalues and $E_{g}$ is a constant representing the ground state energy. At the neutrality point, according to our preceding discussion, $N$ of the $\hat{\psi}_{n}$ eigenmodes coincide with the exact zero modes mandated by the Jackiw-Rossi-Weinberg index theorem. We denote these $\chi_{j}$ with $j=1\dots N$. Because $\varepsilon_{j}=0$ these modes do not contribute to the Hamiltonian (8). The zero mode eigenfunctions $\Phi_{j}({\bm{r}})$ can be chosen as eigenstates of the p-h symmetry generator $\Xi$. They then satisfy the reality condition

[TABLE]

which implies that $\chi_{j}^{\dagger}=\chi_{j}$; the zero modes are Majorana operators.

As noted before the $N$ zero modes are separated by a gap from the rest of the spectrum. As long as $\delta{\cal H}_{\rm FK}$ and ${\cal H}_{\rm int}$ remain small compared to this gap we may construct the effective low-energy theory of the system by simply projecting onto the part of the Hilbert space generated by $N$ Majorana zero modes. In practical terms this is accomplished by inverting Eq. (9) to obtain

[TABLE]

then substituting $\hat{\Psi}_{\bm{r}}$ into $\delta{\cal H}_{\rm FK}$ and ${\cal H}_{\rm int}$ and retaining only those terms that contain zero mode operators $\chi_{j}$ but no finite-energy eigenmodes. We thus obtain

[TABLE]

where

[TABLE]

and $\rho_{ij}({\bm{r}})={i\over 2}\Phi_{i}^{\dagger}({\bm{r}})\tau^{z}\Phi_{j}({\bm{r}})$ is the charge density associated with the pair of zero modes $\chi_{i}$ and $\chi_{j}$. We observe that at the neutrality point when $\tilde{K}_{ij}=0$ the low-energy effective Hamiltonian (12) coincides with the SYK model. Eqs. (13) and (14) allow us to calculate the relevant two- and four-fermion coupling constants from the knowledge of the Majorana wavefunctions in the non-interacting system. We shall carry out this program in Section IV below for a specific physically relevant model system. Here we finish by discussing some general properties of Hamiltonian (12) that follow from symmetry considerations.

The reality condition (10) for the Majorana wavefunction implies the following spinor structure of $\Phi_{j}({\bm{r}})$ in the Nambu space

[TABLE]

where $\eta_{j}(r)$ is a two-component complex spinor. We thus have

[TABLE]

The charge density is thus purely real and antisymmetric under $i\leftrightarrow j$. In the simplest case the $\Sigma$-breaking part of the Fu-Kane Hamiltonian will simply be $\delta H_{\rm FK}({\bm{r}})=-\mu\tau^{z}$. In this situation Eq. (13) implies that $\tilde{K}_{ij}=4\mu\int d^{2}r\rho_{ij}({\bm{r}})$. Thus $\tilde{K}_{ij}$ is purely real and antisymmetric as required for ${\cal H}_{\rm eff}$ to be hermitian.

Because of the anticommutation property (2) of the Majorana operators it is clear that only the fully antisymmetric part of $\tilde{J}_{ijkl}$ contributes to the Hamiltonian (12). As defined in Eq. (14) $\tilde{J}_{ijkl}$ is already antisymmetric under $i\leftrightarrow j$ and $k\leftrightarrow l$ due to the antisymmetry $\rho_{ij}=-\rho_{ji}$. With this in mind we can rewrite Hamiltonian (12) in a more convenient form

[TABLE]

with

[TABLE]

now fully antisymmetric. In the following we will be interested in situations where coupling constants are random and will characterize the coupling strengths by two parameters $K$ and $J$ defined by

[TABLE]

where the bar represents an ensemble average over randomness.

II.3 Structure and statistics of the coupling constants $J_{ijkl}$

In order to approximate the SYK Hamiltonian the coupling constants $J_{ijkl}$ given in the previous subsection must behave as independent random variables. To asess this condition we now discuss their structure and statistics. We make two reasonable assumptions: (i) that the interaction potential in Eq. (14) is short ranged and well approximated by $V({\bm{r}})\simeq V_{0}\delta({\bm{r}})$, and (ii) that there exists a lengthscale $\zeta$ beyond which Majorana wavefunctions $\Phi_{j}({\bm{r}})$ can be treated as random independent variables.

We coarse-grain the Majorana wavefunctions on the grid with with sites ${\bm{r}}_{n}$ and spacing $\sim\zeta$. This amounts to replacing $\eta_{j}({\bm{r}})\to\bar{\eta}_{j}({\bm{r}}_{n})/\zeta$ and $\int d^{2}r\to\zeta^{2}\sum_{n}$ in Eqs. (15) and (14). The discretized spinor wavefunctions then have the following structure on each site

[TABLE]

where $\phi_{j}^{\alpha}(n)$ are real independent random variables with

[TABLE]

Here $M_{s}=\pi R^{2}/\zeta^{2}$ is the total number of grid sites in the hole and the second equality follows from the normalization of $\Phi_{j}({\bm{r}})$.

Combining Eqs. (14), (16), (18) and (20) it is possible to express the antisymmetrized coupling constants as

[TABLE]

where $\epsilon_{\alpha\beta\mu\nu}$ is the totally antisymmetric tensor and summation over repeated indices is implied. For a general value of $M_{s}$ the manybody Hamiltonian defined by coupling constants Eq. (22) represents a variant of the original SYK model similar to models studied in Refs. Bi et al. (2017); Fu et al. (2017). As such it might be amenable to the large-$N$ analysis using approaches described in those works. Here we focus on the limit $M_{s}\gg N$ which attains when the hole radius $R$ is large and the wavefunctions can be considered random on short scales $\zeta$. In this limit each $J_{ijkl}$ defined in Eq. (22) is given by a sum of a large number of random terms given by products of four random amplitudes $\phi_{j}^{\alpha}(n)$. The central limit theorem then assures us that $J$’s will be random variables with a distribution approaching the Gaussian distribution irrespective of the detailed statistical properies of $\phi_{j}^{\alpha}(n)$. It is furthermore easy to show that

[TABLE]

where the uppercase label represents a group of four indices $I=\{ijkl\}$, etc. The coupling constants given by Eq. (22) are asymptotically independent with the higher order correlators vanishing as higher powers of $M_{s}$, e.g. $\overline{J_{ijkl}J_{klmn}J_{mnij}}\sim M_{s}^{-5}$.

The above analysis suggests that under reasonable assumptions coupling constants defining the manybody Hamiltonian (17) can be considered independent random variables. When additionally $K_{ij}$ can be taken as negligible we expect the Hamiltonian to approximate the SYK model. Building on the experience gained from Refs. Bi et al. (2017); Fu et al. (2017) we furthermore expect our Hamiltonian to describe an interesting non-Fermi liquid phase even away from the limit when $J$’s are independent variables. For instance certain specific correlations present in $J$’s are known to lead to a very interesting supersymmetric version of the SYK model Fu et al. (2017) and a whole family of SYK-like models discussed in Ref. Bi et al. (2017) .

Recent work Chew et al. (2017) performed a mathematical analysis of deviations in $J$’s from ideal random independent variables in a model qualitatively similar to ours. Here we adopt a different approach and proceed by evaluating the effect of such deviations on the observable physical properties of the manybody model defined by Hamiltonian (17). We find that coupling constants that follow from the giant vortex geometry indeed give rise to a phenomenology that is consistent with the SYK model.

III The large-$N$ solution and the conformal limit

When the number of Majorana fermions $N$ is large the SYK model becomes analytically solvable in the low-energy limit. Specifically, the Euclidean space time-ordered propagator defined as

[TABLE]

can be expressed in the Matsubara frequency domain through the self energy $\Sigma(\omega_{n})$ as

[TABLE]

Here $G(\omega_{n})=\int_{0}^{\beta}d\tau e^{i\omega_{n}\tau}G(\tau)$ and $\beta=1/k_{B}T$ is the inverse temperature. At non-zero temperatures the propagator and the self energy are defined for discrete Matsubara frequencies $\omega_{n}=\pi T(2n+1)$ with $n$ integer and taking $k_{B}=1$ here and henceforth. Using the replica trick to average over disorder configurations, or alternately summing the leading diagrams in the ${1\over N}$ expansion, one obtains (see for example Ref. Maldacena and Stanford (2016)) the following expression for the self energy appropriate for the Hamiltonian (17)

[TABLE]

For arbitrary given parameters $K$, $J$ and $\beta$ the selfconsistent Eqs. (25) and (26) can be solved by numerical iteration. Analytical solutions are available in various limits and will be reviewed below. In subsequent sections we shall compare these with numerical results based on the model described above.

III.1 Free-fermion limit

When $J=0$ the theory becomes noninteracting and an analytic solution to Eqs. (25) and (26) can be given for all temperatures. Specifically the self energy in Eq. (26) can be written in the frequency domain as $\Sigma(\omega_{n})=K^{2}G(\omega_{n})$ and substituted into Eq. (25). Solving for $G(\omega_{n})$ then gives

[TABLE]

This implies high-frequency limit $G_{f}(\omega_{n})\simeq i/\omega_{n}$ and low-frequency limit $G_{f}(\omega_{n})\simeq i/\mathop{\rm sgn}\nolimits(\omega_{n})K$.

It is useful to extract the single-particle spectral function from Eq. (27) defined as $A(\omega)={1\over\pi}{\rm Im}G(\omega_{n}\to-i\omega+\delta)$, by analytically continuing from Matsubara to real frequencies to obtain the retarded propagator. We thus find

[TABLE]

the usual semicircle law. For this 0-dimensional system $A(\omega)$ coincides with the local density of states $D(\omega)$ averaged over all Majorana sites which is experimentally measurable in a tunneling experiment. Specifically, the tunneling conductance $g(\omega)=(dI/dV)_{\omega=eV}$ is proportional to the local density of states $D(\omega)$.

III.2 Conformal limit

When $K=0$ and $T\ll J$ the system is strongly interacting but nevertheless asymptotic solution of Eqs. (25) and (26) can be found by appealing to their approximate reparametrization invariance Kitaev (2015); Maldacena and Stanford (2016) that becomes exact in the low-frequency limit when one can neglect the $-i\omega_{n}$ term in Eq.  (25). The conformal limit solution reads

[TABLE]

and the corresponding spectral function

[TABLE]

These expressions are valid for $|\omega|\ll J$ and must cross over to the $1/\omega$ behavior at large frequencies.

It is important to note that the low-frequency behaviors of $A_{f}$ and $A_{c}$ are quite different with the former saturating at $1/\pi K$ and the latter divergent. Thus it should be possible to distinguish the free-fermion and the interaction dominated behaviors, illustrated in Fig. 2a, by performing a tunneling experiment. We will discuss the measurement in more detail in Section VI.

III.3 Crossover region

When both $K$ and $J$ are nonzero, as will be the case in a typical experimental setup, analytical solutions are not available but one can still understand the behavior of the system from approximate analytical considerations and numerical solutions. Let us focus on the $T=0$ limit and study the effect of $K$ and $J$ on the self energy $\Sigma(\omega_{n})$ in Eq. (26). To this end it is useful to consider the propagators $G_{f}$ and $G_{c}$ in the imaginary time domain. For long times $\tau$ one obtains

[TABLE]

Consider the $K=0$ limit and then slowly turn $K$ on. Initially, $G_{c}(\tau)$ is a valid solution. However, for any non-zero $K$ it is clear that the first term on the right-hand side of Eq. (26) will dominate at sufficiently long times $\tau>\tau_{*}$. At such long times one then expects a crossover to the behavior resembling the free-fermion propagator $G_{f}(\tau)$. The corresponding crossover time $\tau_{*}$ can be estimated by equating the two terms on the right-hand side of Eq. (26), $K^{2}G_{f}(\tau_{*})=J^{2}G_{c}^{3}(\tau_{*})$, which gives

[TABLE]

and the corresponding crossover frequency

[TABLE]

We thus expect the spectral function to behave as indicated in Eq. (30) for $\omega_{*}<\omega\ll J$ with the divergence at small $\omega$ cut off below $\omega_{*}$ and saturate to $\sim 1/\pi K$.

To confirm the above behavior we have solved Eqs. (25) and (26) numerically. We found it most convenient to work with Matsubara Green’s functions at very low but non-zero temperatures. To this end we rewrite Eq.  (26) in Matsubara frequency domain where the last term becomes a convolution and substitute the self energy into Eq.  (25). We obtain a single equation

[TABLE]

for $G_{n}\equiv G(\omega_{n})$ that must be solved selfconsistently. Results obtained by iterating Eq. (34) are displayed in Fig. 2b. For very small $K=0.01J$ we observe that numerically calculated $G(\omega_{n})$ coincides with the conformal limit for a range of frequencies consistent with our discussion above. For $K=0.1J$ this range becomes smaller and completely disappears for $K=0.5J$.

We conclude that for any nonzero $K$ the ultimate low-energy behavior is controlled by the free-fermion fixed point, as expected on general grounds. Nevertheless, when $K$ is sufficiently small in comparison to $J$, there can be a significant range of energies in which the physics is dominated by the SYK fixed point. At low temperatures the corresponding range of frequencies is given by

[TABLE]

In this range we expect the spectral function to obey the conformal scaling form given by Eq. (30). Tunneling experiment in this regime should therefore reveal the SYK behavior of the underlying strongly interacting system.

IV Numerical results: the underlying noninteracting system

In this section we provide support for the ideas presented above by performing extensive numerical simulation and modeling of the system described in Sec. II. We start by formulating a lattice model for the surface of a TI in contact with a superconductor. We then find the wavefunctions of the Majorana zero modes by numerically diagonalizing the corresponding Bogoliubov-de Gennes (BdG) Hamiltonian for the geometry depicted in Fig. 1 with $N$ flux quanta threading the hole. In the following Section, using Eqs. (13) and (14), we calculate the coupling constants $K_{ij}$ and $J_{ijkl}$, which we then use to construct and diagonalize the many-body interacting Hamiltonian (17) for $N$ up to 32. The resulting many-body spectra and eigenvectors are used to calculate various physical quantities (entropy, specific heat, two- and four-point propagators) which are then compared to the results previously obtained for the SYK model with random independent couplings.

IV.1 Lattice model for the TI surface

A surface of a 3D TI is characterized by an odd number of massless Dirac fermions protected by time reversal symmetry $\Theta$. The well known Nielsen-Ninomyia theorem Nielsen and Ninomiya (1981a, b) assures us that, as a matter of principle, it is impossible to construct a purely 2D, $\Theta$-invariant lattice model with an odd number of massless Dirac fermions. This fact causes a severe problem for numerical approaches to 3D TIs because one is forced to perform an expensive simulation of the 3D bulk to describe the anomalous 2D surface. A workaround has been proposed Marchand and Franz (2012) which circumvents the Nielsen-Ninomyia theorem by simulating a pair of TI surfaces with a total even number of Dirac fermions. This approach enables efficient numerical simulations in a quasi-2D geometry while fully respecting $\Theta$.

Here, because the physical time reversal symmetry is ultimately broken by the presence of vortices and is therefore not instrumental, we opt for an even simpler model which breaks $\Theta$ from the outset but nevertheless captures all the essential physics of the TI/SC interface. We start from the following momentum-space normal-state Hamiltonian defined on a simple 2D square lattice

[TABLE]

with $M_{\bm{k}}=m[(2-\cos{k_{x}}-\cos{k_{y}})-{1\over 4}(2-\cos{2k_{x}}-\cos{2k_{y}})]$. Here ${\bm{\sigma}}$ are Pauli matrices in spin space and $\lambda$, $m$ are model parameters. The term proportional to $\lambda$ respects $\Theta$ and gives 4 massless Dirac fermions in accordance with the Nielsen-Ninomyia theorem. The $M_{\bm{k}}$ term breaks $\Theta$ and has the effect of gapping out all the Dirac fermions except the one located at $\Gamma=(0,0)$. The resulting energy spectrum

[TABLE]

is depicted in Fig. 3. In the vicinity of the $\Gamma$ point we observe a linerly dispersing spectrum characteristic of a TI surface state. It is to be noted that for small $|{\bm{k}}|$ we have $M_{\bm{k}}\simeq{1\over 8}mk^{4}$ so the amount of $\Theta$-breaking can be considered small in the physically important part of the momentum space near the $\Gamma$ point.

Proximity induced superconducting order is implemented by constructing the BdG Hamiltonian,

[TABLE]

Writing $H_{\rm BdG}$ in terms of $\sigma$ and $\tau$ matrices it can be easily checked that it respects the particle-hole symmetry $\Xi$ defined in Sec. II.A. The $\mu$ and $M_{\bm{k}}$ terms both break the fictitious time reversal $\Sigma$ that protects the Majorana zero modes in our setup. As before $\mu$ must be tuned to zero to achieve protection. On the other hand it is crucial to remember that $M_{\bm{k}}$ has been introduced only to circumvent the Nielsen-Ninomyia theorem and allow us to efficiently simulate a single two-dimensional Dirac fermion on the lattice. Breaking of $\Sigma$ by $M_{\bm{k}}$ is therefore not a concern in the experimental setup: in a real TI tuned to the neutrality point $\Sigma$ is unbroken. Expanding $H_{\rm BdG}({\bm{k}})$ in the vicinity of $\Gamma$ to leading order in small ${\bm{k}}$ we recover the Fu-Kane Hamiltonian $H_{\rm FK}$ defined in Eq. (4). We thus conclude that at low energies our lattice model indeed describes the TI/SC interface and should exhibit the desired phenomenology, including Majorana zero modes bound to vortices. We shall see that this is indeed the case. The only repercussion that follows from the weakly broken $\Sigma$ (present in the higher order terms in the above expansion) is a very small splitting of the zero mode energies that has no significant effect on our results.

IV.2 Solution in the giant vortex geometry

To study the non-uniform system with magnetic field and vortices we must write the Hamiltonian in the position space. The normal-state piece Eq. (36) is most conveniently written in second-quantized form as

[TABLE]

where we defined on each lattice site ${\bm{r}}$ a two-component spinor $\psi_{\bm{r}}=(c_{{\bm{r}}\uparrow},c_{{\bm{r}}\downarrow})^{T}$ and $\alpha=x,y$. The magnetic field is included through the standard Peierls substitution which replaces tunneling amplitudes on all bonds according to $\psi_{\bm{r}}^{\dagger}\psi_{{\bm{r}}+\alpha}\to\psi_{\bm{r}}^{\dagger}\psi_{{\bm{r}}+\alpha}\exp(-i{e\over\hbar c}\int_{\bm{r}}^{{\bm{r}}+\alpha}d{\bm{l}}\cdot{\bm{A}})$. The full second quantized BdG Hamiltonian then reads

[TABLE]

where $\Delta_{\bm{r}}$ is the pair potential on site ${\bm{r}}$ which takes the form indicated in Eq. (5). In accord with our discussion in the previous subsection ${\cal H}_{\rm BdG}$ given in Eq. (40) represents a version of the Fu-Kane Hamiltonian (3) regularized on a square lattice. This lattice model is suitable for numerical calculations and we expect it to reproduce all the low-energy features of the Fu-Kane Hamiltonian. In particular we will see shortly that it yields $N$ Majorana zero modes mandated by the Jackiw-Rossi-Weinberg index theorem that are of central importance for the SYK model.

It is most convenient to solve the problem defined by Hamiltonian (40) on a lattice with $L\times L$ sites and periodic boundary conditions which ensure that no spurious edge states exist at low energies in addition to the expected $N$ Majorana zero modes bound to the hole. To implement periodic boundary conditions it is useful to perform a singular gauge transformation

[TABLE]

which has the effect of removing the phase winding from $\Delta_{\bm{r}}$ and changing the Peierls phase factors to $\exp(i\int_{\bm{r}}^{{\bm{r}}+\alpha}d{\bm{l}}\cdot{\bm{\Omega}})$ with

[TABLE]

We note that $N$ must be even because only for integer number of fundamental flux quanta $hc/e=2\Phi_{0}$ in the system one can impose periodic boundary conditions. For $N$ even the transformation (41) is single valued and the issue of branch cuts that renders the analogous problem with singly-quantized vortices Franz and Tešanović (2000); Vafek et al. (2001) more complicated does not arise here. After the transformation the total effective flux seen by the electrons $\int dS(\nabla\times{\bm{\Omega}})_{z}$ vanishes and numerical diagonalization of the transformed Hamiltonian (40) with periodic boundary conditions becomes straightforward.

As a practical matter it is easiest to define a regular shaped hole and introduce disorder through a replacement

[TABLE]

Here $(\delta\mu_{\bm{r}},\delta\Delta_{\bm{r}},\delta\lambda_{\bm{r}})$ are independent random variables uniformly distributed in the interval $(-w_{\mu}/2,w_{\mu}/2)$ for $\delta\mu_{\bm{r}}$ and similarly for $\delta\Delta_{\bm{r}}$ and $\delta\lambda_{\bm{r}}$. We chose a stadium-shaped hole sketched in Fig. 4a which is known to support classically chaotic trajectories Gutzwiller (1990); Haake (2013). In our quantum simulation we find that much smaller disorder strength is required to achieve sufficiently random Majorana wavefunctions for stadium-shaped hole than e.g. a with circular hole. We furthermore chose magnetic field $B$ to be uniform inside the radius $R_{B}$ that contains the hole and zero otherwise. We find that our results are insensitive to the detailed distribution of $B$ as long as the total flux remains $N\Phi_{0}$ and is centered around the hole (we tested various radii $R_{B}$ as well as a Gaussian profile).

Typical results of the numerical simulations described above are displayed in Fig. 4. In panel (b) we observe the behavior of the energy eigenvalues $E_{n}$ of ${\cal H}_{\rm BdG}$. For zero magnetic flux there are several states inside the SC gap (Andreev states bound to the hole) but no zero modes. For $N=24$ these are converted into 24 zero modes required by the Jackiw-Rossi-Weinberg index theorem. For $\mu=w_{\mu}=0$ used in the simulation their energies are very close to zero ($\sim 10^{-4}\lambda$), where the small residual splitting is attributable to the fact that $\Sigma$ symmetry is weakly broken in our lattice simulation by the $M_{\bm{k}}$ term. For non-zero $\mu$ or $w_{\mu}$ the energy splitting increases in proportion to these $\Sigma$-breaking perturbations. In the following we include these terms in $\delta{\cal H}_{\rm FK}$ and incorporate them in our many-body calculation via $K_{ij}$ terms given by Eq. (13).

Panels Fig. 4 c-f show examples of zero mode wavefunction amplitudes $|\Phi_{j}({\bm{r}})|^{2}$. The wavefunctions are seen to exhibit random spatial structure which depends sensitively on the specific disorder potential realization. Importantly all the zero mode wavefunctions are localized in the same region of space defined by the hole and its immediate vicinity. One therefore expects Eq. (14) to produce strong random couplings $J_{ijkl}$ connecting all zero modes $\chi_{j}$ once the interactions are included.

V Numerical results: the many-body SYK problem

Having obtained the zero mode wavefunctions it is straightforward to calculate couplings $K_{ij}$ and $J_{ijkl}$ from Eqs. (13) and (14) and construct the many-body SYK Hamiltonian (17). In the following we shall assume that the system has been tuned to its global neutrality point $\mu=0$ and include in $\delta{\cal H}_{\rm FK}$ only the random component of the on-site potential $\delta\mu_{\bm{r}}$. For the interaction term we consider the screened Coulomb potential defined as

[TABLE]

where $\epsilon$ is the dielectric constant and $\lambda_{\rm TF}$ denotes the Thomas-Fermi screening length. We furthermore assume that $\lambda_{\rm TF}$ is short so that in the lattice model the interaction is essentially on-site. The expression for $\tilde{J}_{ijkl}$ then simplifies to

[TABLE]

with $V_{0}=\int d^{2}rV(r)=2\pi e^{2}\lambda_{\rm TF}/\epsilon$. Coupling constants $K_{ij}$ and $J_{ijkl}$ are easy to evaluate using Eq. (18) and the Majorana wavefunctions $\Phi_{j}({\bm{r}})$ obtained in the previous Section. To facilitate comparisons with the existing literature we shall quantify the average strength of these terms using parameters $K$ and $J$ defined in Eq. (19). Specifically we will adjust $w_{\mu}$ and $V_{0}$ to obtain the desired values of $K$ and $J$. In the following Section we will connect these values to the parameters expected in realistic physical systems.

V.1 Thermodynamic properties and many-body level statistics

Once the coupling constants $K_{ij}$ and $J_{ijkl}$ are determined as described above one can construct a matrix representation of the many-body Majorana Hamiltonian (17) and find its energy eigenvalues $E_{n}$ by exact numerical diagonalization. From the knowledge of the energy levels it is straightforward to calculate any thermodynamic property. In Fig. 5 we display the thermal entropy $S(T)$ and the heat capacity $C_{V}(T)$. These are calculated from

[TABLE]

where $\langle E^{\alpha}\rangle={1\over Z}\sum_{n}E_{n}^{\alpha}e^{-E_{n}/T}$, $F=-T\ln{Z}$ is the free energy and $Z=\sum_{n}e^{-E_{n}/T}$ the partition function.

The entropy per particle is seen to saturate at high temperature to $S_{\infty}/N={1\over 2}\ln{2}\simeq 0.3465$ as expected for a system of $N$ Majorana fermions. The behavior of $S(T)$ calculated for the giant vortex system is qualitatively similar to that obtained from the SYK model with random independent couplings. The small deviations that exist are clearly becoming smaller as $N$ grows, suggesting that they vanish in the thermodynamic limit. Nonzero two-body coupling $K$ is seen to modify the entropy slightly at low temperature. For large $N$ and $K=0$ the entropy per particle is expected to attain a nonzero value $\sim 0.24$ as $T\to 0$ due to the extensive ground-state degeneracy of the SYK model. Our largest system is not large enough to show this behavior (in agreement with previous numerical results) although Fig. 5a correctly captures the expected suppression of the low-$T$ entropy in the presence of two-body couplings which tend to remove the extensive ground-state degeneracy.

The heat capacity $C_{V}(T)$, displayed in Fig. 5b, likewise behaves as expected for the SYK model with random independent couplings with small deviations becoming negligible in the large-$N$ limit. $C_{V}(T)$ is in principle measurable and we can see from Fig. 5b that its high-temperature behavior could be used gauge the relative strength of two- and four-fermion terms in the system.

As discussed in Refs. García-García and Verbaarschot (2016); You et al. (2016) many-body level statistics provides a sensitive diagnostic for the SYK physics encoded in the Hamiltonian (17). To apply this analysis to our results we arrange the many-body energy levels in ascending order $E_{1}<E_{2}<\dots$ and form ratios between the successive energy spacings

[TABLE]

According to Refs. García-García and Verbaarschot (2016); You et al. (2016) the SYK Hamiltonian can be constructed as a symmetric matrix in the Clifford algebra ${\cal C}\ell_{0,N-1}$ whose Bott periodicity gives rise to a Z8 classification with topological index $\nu=N\!\mod{8}$. As a result statistical distributions of the ratios $P(r)$ cycle through Wigner-Dyson random matrix ensembles with Z8 periodicity as a function of $N$. Specifically, Gaussian orthogonal (GOE), unitary (GUE) and symplectic (GSE) ensembles occur with distributions given by the “Wigner surmise”

[TABLE]

and parameters $Z$ and $\beta$ summarized in Table I for even $N$ relevant to our system.

As emphasized in Ref. You et al. (2016) the level spacing analysis must be performed separately in the two fermion parity sectors of the Hamiltonian (17).

Fig. 6 shows statistical distributions of the ratios $r_{n}$ computed for $N=24,26,28,30$ and 32 in our system. For the sake of clarity $P(\ln{r})$ is plotted along with the anticipated distributions for GOE, GUE and GSE given in Eq. (48). Unambiguous agreement with the pattern indicated in Table I is observed, lending further support to the notion that our proposed system realizes the SYK model. We checked that the Z8 periodic pattern persists for all $N$ down to 16. Additionally, the above results should be contrasted with the level statistics in the noninteracting case $J=0$, $K=1$ displayed in the bottom row of Fig. 6. In the absence of interactions Z8 periodicity is absent and the distribution of the ratio $r_{n}$ follows Poisson level statistics

[TABLE]

for all $N$. It is to be noted that no adjustable parameters are employed in the level-statistics analysis presented above.

V.2 Green’s function

Computing the Green’s function of the model is perhaps the most straightforward way of comparing the behavior of the system at finite $N$ to the large-$N$ limit solutions discussed in Sec. III. At the same time computation of propagators is numerically more costly because in addition to many-body energy levels one requires the corresponding eigenstates. We computed the on-site retarded Green’s function defined as

[TABLE]

Fourier transforming and using Lehmann representation in terms of the eigenstates $|n\rangle$ of the many-body Hamiltonian (17) one obtains, at $T=0$,

[TABLE]

where $\delta$ is a positive infinitesimal. From Eq. (51) the spectral function $A_{i}(\omega)={1\over\pi}{\rm Im}G_{i}^{R}(\omega)$ is readily extracted.

In Fig. 7 we display spectral function $A(\omega)={1\over N}\sum_{i}A_{i}(\omega)$ averaged over all Majorana zero modes. Physically this corresponds to a tunneling experiment with a large probe that allows for tunneling into all sites inside the hole. In agreement with the existing numerical results on the complex fermion version of the SYK model Fu and Sachdev (2016) we find that for system sizes we can numerically access (up to $N=30$) the conformal limit is approached only in a narrow interval of frequencies. In the low-frequency limit numerical results approach a constant value instead of the the $\sim 1/\sqrt{\omega}$ divergence expected in the large-$N$ limit. The dependence on $N$ is very weak, with the larger values showing reduced statistical fluctuations but otherwise qualitatively similar behavior. To convincingly demonstrate the conformal scaling of the Green’s function at the lowest frequencies numerical calculations using larger values of $N$ would be necessary. Unfortunately these are currently out of reach for the exact diagonalization method due to the exponential growth of the Hamiltonian matrix size with $N$. More sophisticated numerical techniques, such as the quantum Monte Carlo, could possibly reach larger system sizes.

Spectral functions calculated for the giant vortex setup exhibit larger statistical fluctuations compared to those computed with random Gaussian coupling constants $J_{ijkl}$ but are qualitatively similar when averaged over independent disorder realizations. Therefore, we conclude that the Green’s function behavior at finite $N$ supports the notion that our proposed system realizes the SYK model.

V.3 Out-of-time-order correlators and scrambling

Scrambling of quantum information – a process in which quantum information deposited into the system locally gets distributed among all its degrees of freedom – is central to the conjectured duality between the SYK model and AdS2 Einstein gravity. Black holes are thought to scramble with the maximum possible efficiency: they exhibit quantum chaos. For a quantum theory to be the holographic dual of a black hole its dynamics must exhibit similar fast scrambling behavior.

Out-of-time-order correlator (OTOC), defined in our system as

[TABLE]

allows to quantify the quantum chaotic behavior. For black holes in Einstein gravity scrambling occurs exponentially fast with $1-F(t)\sim e^{\lambda_{L}t}/N$ where the decay rate is given by the Lyapunov exponent $\lambda_{L}=2\pi T$ Maldacena et al. (2016). Similarly, for the SYK model in the large-$N$ limit one expects Kitaev (2015); Maldacena and Stanford (2016)

[TABLE]

Previous works Hosur et al. (2016); Fu and Sachdev (2016) gave numerical evaluations of $F(t)$ in the SYK model for $N$ up to 14 but found these system sizes to be too small to clearly show the expected $J$-independent Lyapunov exponent. Here we numerically evaluate OTOC for $N$ up to 22 and show that coupling constants obtained from the giant vortex geometry give qualitatively the same behavior as those for random independent coupling constants. Our results are summarized in Fig. 8 where we compute the on-site OTOC $F_{ii}(t)$ averaged over all sites.

For $J=1$ the OTOC is seen to rapidly decay to zero consistent with previous works on the SYK model Hosur et al. (2016); Fu and Sachdev (2016). The rate of decay is controlled by $J$: as in Refs. Hosur et al. (2016); Fu and Sachdev (2016) we find that $N=22$ is not large enough to observe the theoretically predicted $J$-independent Lyapunov exponent controlled by temperature, even when $J\gg T$. In addition we observe that adding sizable two-body tunneling term $K$ has only very modest effect on the behavior of $F(t)$ when the interaction strength is maintained. However, in the non-interacting case $(J=0,\ K=1)$ OTOC behavior changes qualitatively with the fast decay replaced by oscillations whose amplitude slowly increases.

VI Outlook: Towards the experimental realization and detection of the SYK model

Our theoretical results presented above indicate that low-energy fermionic degrees of freedom in a device with geometry depicted in Fig. 1 provide a physical realization of the SYK model. Additionally, all the ingredients are currently in place to begin experimental explorations of the proposed system. Superconducting order has been induced and observed at the surface of several TI compounds by multiple groups Koren et al. (2011); Sacépé et al. (2011); Qu et al. (2012); Williams et al. (2012); Cho et al. (2013); Xu et al. (2014); Zhao et al. (2015).

Importantly, Ref.  Cho et al. (2013) already demonstrated the ability to tune the chemical potential in (BixSb2-x)Se3 thin flakes through the neutrality point in the presence of superconductivity induced by Ti/Al contacts by a combination of chemical doping (tuning $x$) and back gate voltage. This is almost exactly what we require for the implementation of the SYK model. Well developed techniques (such as focused ion milling) exist to fabricate patterns, such as a hole with an irregular shape, in a SC film deposited on the TI surface. In the remainder of this Section we discuss in more detail the experimentally relevant constraints on the proposed device as well as possible ways to detect manifestations of the SYK physics in a realistic setting.

VI.1 Device geometry, length and energy scales

The key controllable design feature is the size of the hole, parametrized by its radius $R$. For simplicity in the estimates below we shall assume a circular hole but it should be understood that in a real experiment irregular shape is required to promote randomness of the zero mode wavefunctions. For the desired number $N$ of Majorana zero modes the hole must be large enough to pin $N$ vortices. Vortex pinning occurs because the SC order parameter $\Delta$ is suppressed to zero in the vortex core which costs condensation energy. Vortices therefore prefer to occupy regions where $\Delta$ has been locally suppressed by defects or in our case by an artificially fabricated hole. The optimal hole size $R_{N}$ for $N$ vortices in our setup can be thus estimated from the requirement that all the electronic states inside the hole that reside within the SC gap are transformed into zero modes,

[TABLE]

where $D(\varepsilon)=|\varepsilon|/2\pi v_{F}^{2}\hbar^{2}$ is the density of states of the TI surface. This gives

[TABLE]

with $\xi=\hbar v_{F}/\pi\Delta$ the BCS coherence length. In the absence of interactions a hole of this size will produce an energy spectrum similar to that depicted in Fig. 4b, with $N$ zero modes maximally separated from the rest of the spectrum.

In reality, if the SC film is in the type-II regime, a somewhat larger hole might be required to reliably pin $N$ vortices in a stable configuration and not create vortices nearby. The latter condition is that $B<B_{c1}$, where $B_{c1}$ is the lower critical field. Thus, the magnetic field to get the necessary flux is

[TABLE]

where $\lambda_{\rm eff}$ is the effective penetration depth of a thin SC film defined below Eq. (5). This gives

[TABLE]

Taking the standard expression for the lower critical field $B_{c1}=(\Phi_{0}/4\pi\lambda_{\rm eff}^{2})K_{0}(\kappa_{\rm eff}^{-1})$, where $\kappa_{\rm eff}=\lambda_{\rm eff}/\xi$, Eq. (57) becomes

[TABLE]

In type-II regime $\lambda_{\rm eff}>\xi$ and Eq. (58) will generally imply larger hole size than Eq. (55). A larger hole size would reduce the spectral gap to some extent but $N$ zero modes will remain robustly present. If the SC film remains in the type-I regime then there is no additional constraint on $R_{N}$ but the applied field must be kept below the thermodynamic critical field $B_{c}$ of the film.

These considerations impose some practical constraints on the material composition and thickness $d$ of the SC film. In general we want the film to be sufficiently thin so that scanning tunneling spectroscopy of the hole region can be performed. On the other hand we want it to be either in type-I or weakly type-II regime such that Eq. (58) does not enlarge the hole size significantly beyond the ideal radius given by Eq. (55). For Pb we have $(\xi,\lambda_{L})=(83,37)$nm. Taking $d=20$nm results in $\lambda_{\rm eff}\simeq 137$nm and Eq. (58) imposes only a mild increase in the hole size compared to the ideal, which should not adversely affect the zero modes. For Al we have $(\xi,\lambda_{L})=(1600,16)$nm and one can go down to very thin films and still remain in the type-I regime.

The TI film must be sufficiently thick so that it exhibits well developed gapless surface states. For Bi2Se3 family of materials this means thickness larger than 5 unit cells. TI films close to this critical thickness will also be easiest to bring to the neutrality point by back gating.

Using a hole close to the ideal size given by Eq. (55) will also promote the interaction strength. Intuitively it is clear that screened Coulomb interaction between electrons will have maximum effect on the zero modes if their wavefunctions are packed as closely together as possible. With this in mind one can give a crude estimate of the expected interaction strength $J$ as follows. Starting from Eq. (22) with $V_{0}=2\pi e^{2}\lambda_{\rm TF}/\epsilon$ and using Eq. (21) it is easy to show that

[TABLE]

where we identified the lengthscale $\zeta$ with the SC coherence length $\xi$. We can obtain a physically more transparent expression by introducing the Bohr radius $a_{0}=\hbar^{2}/m_{e}e^{2}\simeq 0.52$Å and the corresponding Rydberg energy $E_{0}=e^{2}/2a_{0}\simeq 13.6$eV,

[TABLE]

Several remarks are in order. Eq. (60) implies that for a fixed hole size $R$ the coupling strength grows as $J\sim N^{3/2}$. It is therefore advantageous to put as many flux quanta in the hole as can be stabilized. For the ‘ideal’ hole size $R=R_{N}$ given by Eq. (55) we have $M_{s}=2\pi^{3}N$ and the dependence on $N$ drops out. The amplitude of $J$ will then depend only on the coherence length $\xi$, screening length $\lambda_{\rm TF}$ and dielectric constant $\epsilon$ of the system. To get an idea about the possible size of $J$ we assume $\lambda_{\rm TF}\approx\xi$ and $\epsilon\approx 50$, appropriate for the surface of a TI such as Bi2Se3. Eq. (60) then gives $J\approx(1{\rm\AA}/\xi)17.8$meV. It is clear that using a superconductor with a large gap and short coherence length would aid the observation of the SYK physics in this system at reasonable energy and temperature scales. Taking Pb as a concrete example we have $\xi\simeq 52$nm, for $d>20$nm Eq. (58) does not impose additional restrictions on the hole diameter, and one obtains $J$ in the range of tens of $\mu$eV. This energy scale is accessible to scanning tunneling spectroscopy (STS) which, as we argue below, constitutes the most convenient experimental probe.

VI.2 Experimental detection

In our proposed setup the experimental detection of the signatures of the SYK state can be achieved using tunneling spectroscopy. Either planar tunneling measurement with a fixed probe weakly coupled to the TI surface or a scanning tunnel probe can be used. STS has the advantage of simultaneously being able to image the topography of the device with nanoscale resolution and measure the tunneling conductance $g(\omega)$ which is proportional to the spectral function of the system $A(\omega)$. A recently developed technique Ge et al. (2016) combines an STS tip with a miniature Hall probe which allows additional measurement of the local magnetic field $B$ at the sample surface. Such a probe is ideally suited for the proposed SYK model setup as it can be used to independently determine the magnetic flux and thus the number $N$ of Majorana fermions in the system.

In the large-$N$ limit of the SYK model $A(\omega)$ exhibits the characteristic $1/\sqrt{|\omega|}$ singularity (illustrated in Fig. 2a) which should be easy to distinguish from the semicircle distribution that prevails in a system dominated by random two-fermion tunneling terms. In the large-$N$ limit and at sufficiently low temperature $k_{B}T\ll J$ the detection of the SYK behavior via tunneling spectroscopy should therefore be relatively straightforward.

In a realistic setup the number of flux quanta $N$ will be finite and perhaps not too large. In this case our results in Fig. 7 show that the characteristic $1/\sqrt{|\omega|}$ singularity is cut off such that $A(0)$ is finite and grows with $N$ only very slowly. Additional results assembled in Fig. 9 indicate that even in this situation it is possible to distinguish the interaction-dominated SYK behavior from the behavior characteristic of the weakly interacting system with random two-fermion couplings. For $J\gtrsim K$ we observe a relatively smooth spectral density peaked at $\omega=0$ characteristic of the strongly interacting regime. In the opposite limit $J\lesssim K$ non-universal fluctuations that strongly depend on the specific disorder realization become increasingly prominent. Eventually, when $J\ll K$, the spectral function consists of a series on $N$ sharp peaks. These peaks occur at the eigenvalues of the $N\times N$ random hermitian matrix $iK_{ij}$ and represent the single-particle excitations of the non-interacting problem at $J=0$. At large $N$ these peaks merge to form a continuous distribution described by the semicircle law.

Full numerical diagonalization of the SYK Hamiltonian is feasible for $N$ up to 32 on a desktop computer and involves a matrix of size $2^{15}\times 2^{15}$ in each parity sector. By going to a supercomputer one can plausibly reach $N=42$, Boixo et al. (2016) but larger system sizes are out of reach due to the exponential growth of the Hamiltonian matrix with $N$. Experimental realization using the setup proposed here has no such limitation. Measurement of the spectral function in such a system could therefore help elucidate the approach to the large $N$ limit in which the SYK model becomes analytically tractable by field theory techniques. This has relevance to the spontaneous breaking of the emergent conformal symmetry at large $N$ and a host of other interesting issues extensively discussed in the recent literature Kitaev (2015); Maldacena and Stanford (2016); Sachdev (2015); Maldacena et al. (2016); Hosur et al. (2016); Polchinski and Rosenhaus (2016); García-García and Verbaarschot (2016); You et al. (2016). Measurement of the out-of-time-oder correlator $F(t)$ for $N$ larger than 32 could furthermore shed light on the emergence of the quantum chaotic behavior in the system, scrambling and the dual relation to the extremal black hole in AdS2. A protocol to measure $F(t)$ in a system of this type is currently unknown and this represents an interesting challenge and an opportunity for future study.

VII Conclusions

To conclude, we proposed a physical realization of the Sachdev-Ye-Kitaev model that utilizes available materials and experimental techniques. The proposal is to use the surface of a 3D TI at its global neutrality point proximitized by a conventional superconductor with an irregular-shaped hole and magnetic flux threaded through the hole. We demonstrated that the conventional screened Coulomb interaction between electrons in such a setup leads to a Majorana fermion Hamiltonian at low energies with requisite random four-fermion couplings. Detailed analysis indicates behavior consistent with that expected of the SYK model. We gave estimates for model parameters in the realistic systems and suggested experimental tests for the SYK behavior. This work thus provides connections between seemingly unrelated areas of research – mesoscopic physics, spin liquids, general relativity, and quantum chaos – and could lead to experimental insights into phemomena that are of great current interest.

Acknowledgements.

The authors are indebted to J. Alicea, O. Can, A. Kitaev, E. Lantagne-Hurtubise, M. Rozali, C.-M. Jian, I. Martin, and S. Sachdev for illuminating discussions. We thank NSERC, CIfAR and Max Planck - UBC Centre for Quantum Materials for support. DIP is grateful to KITP, where part of the research was conducted with support of the National Science Foundation Grant No. NSF PHY11-25915. Numerical simulations were performed in part with computer resources provided by WestGrid and Compute Canada Calcul.