Creating and probing the Sachdev-Ye-Kitaev model with ultracold gases: Towards experimental studies of quantum gravity

We suggest that the holographic principle, combined with recent technological advances in atomic, molecular, and optical physics, can lead to experimental studies of quantum gravity. As a specific example, we consider the Sachdev-Ye-Kitaev (SYK) model, which consists of spin-polarized fermions with an all-to-all complex random two-body hopping and has been conjectured to be dual to a certain quantum gravitational system. Achieving low-temperature states of the SYK model is interpreted as a realization of a stringy black hole, provided that the holographic duality is true. We introduce a variant of the SYK model, in which the random two-body hopping is real. This model is equivalent to the origincal SYK model in the large-$N$ limit. We show that this model can be created in principle by confining ultracold fermionic atoms into optical lattices and coupling two atoms with molecular states via photo-association lasers. This development serves as an important first step towards an experimental realization of such systems dual to quantum black holes. We also show how to measure out-of-time-order correlation functions of the SYK model, which allow for identifying the maximally chaotic property of the black hole.


Introduction
The quantum nature of black holes is one of the most important subjects in theoretical physics, since the theoretical discovery of particle-emissions from a black hole due to quantum effects [1,2], which are often referred to as the Hawking radiation. Although there have been experimental searches for quantum black holes at the CERN LHC motivated by the predictions on the basis of theories of TeV-scale quantum gravity [3][4][5], no evidence of the black hole creation has been observed thus far [6][7][8][9]. In this paper, we present a completely different route to experimental studies of quantum gravity by exploiting both holographic principle and unprecedented controllability of optical-lattice systems loaded with ultracold gases [10].
In order to resolve paradoxes associated with the black hole evaporation that results from the Hawking radiation, the holographic principle [13,14] emerged, which claims that black holes, and more general quantum gravitational theories, are equivalent to non-gravitational theories in different spacetime dimensions. As a concrete example, the gauge/gravity duality conjecture [15], which claims the duality (i.e. the equivalence) between superstring/M-theory on certain spacetimes and quantum field theories, has been studied extensively. Although this conjecture has not been proven yet, it is believed to be correct at least in some simplest cases. For example, maximally supersymmetric matrix quantum mechanics (also known as the Matrix Model of M-theory [16,17]), which is conjectured to describe a black hole in type IIA superstring theory near the 't Hooft large-N limit [18], has been studied numerically starting in [19]. The agreement with the dual superstring theory prediction has been confirmed including the effect of virtual loops of string [20].
Thanks to their high controllability and cleanness, experiments with ultracold gases in optical lattices have succeeded in realizing various theoretical models, which were introduced in the contexts of condensed matter physics but did not have quantitative experimental counterparts. Examples include the Bose-Hubbard model [21], the Lieb-Liniger model [22,23], the Aubry-André model [24], the Harper Hamiltonian [25,26], and the topological Haldane model [27]. There have been theoretical proposals also for realizing lattice gauge models studied in high-energy physics [28][29][30]. These circumstances tempt one to expect that it may be possible as well to realize quantum field theories dual to quantum gravitational systems.
In this paper, we propose a possible way to create the Sachdev-Ye-Kitaev (SYK) model [31][32][33][34] experimentally with use of ultracold gases in optical lattices. The SYK model consists of spin-polarized fermions with an all-to-all random two-body hopping. Its thermal state is a non-Fermi liquid with nonzero entropy at vanishing temperature, which is called the with σ 2 g = (2N) −3 J 2 while ν s = + √ n ms J for even s and ν s = − √ n ms J for odd s.
model. Due to this, each Feynman diagram receives some correction after disorder average. However, such corrections are 1/N-suppressed in general, and hence this modified model agrees with the original model at large-N. In the following, we call the original SYK model with complex J ij,kl and the modified one with real J ij,kl 'complex-SYK' and 'real-SYK,' respectively. The 1/N-corrections to the real-and complex-SYK models are described by different sets of Feynman diagrams. In analogy to the duality between gauge theory and superstring, it is natural to expect that these two theories describes slightly different quantum gravitational systems whose classical limits coincide. In Appendix A, we indeed perform numerical comparisons between the two models to demonstrate that the computed physical quantities of the two models rapidly approaches each other as N increases. One of the severest bottlenecks for realizing the SYK model in optical-lattice experiments is the implementation of the all-to-all two-body hopping, because particles on lattice systems in general move the most dominantly via nearest-neighbor one-body hopping. In order to overcome this bottleneck, we consider a situation, in which two atoms are coupled with n ms molecular states, described by the following Hamiltoinan, Here, ν s , U s,s ′ , and g s,ij denote the detuning of molecular state s, the onsite interaction between two molecules in states s and s ′ , and the atom-molecule coupling constant. Using the degenerate perturbation theory up to the second order, we obtain the following effective See Appendix C for more detailed derivations of the effective Hamiltonian. A similar way of designing a kind of two-body hopping term, namely the ring exchange interaction, by means of intermediate two-particle states has been pointed out in previous work [42]. In the next section we elaborate how to prepare such a situation in optical-lattice systems while in this section we show that Eq. (11) serves as a quantitative approximation of the complex SYK model (1) when n ms is sufficiently large and ν s is appropriately tuned.
Let us suppose ν 1 = ν 2 = · · · = ν nms ∝ √ n ms . Then, if n ms is large enough, s g s,ij g s,kl νs should become Gaussian except for the diagonal elements (i, j) = (k, l) or (i, j) = (l, k) (note that g 2 s,ij is always positive). This happens because it is simply an n ms -step random walk for each set of indices (i, j, k, l). In order to improve the behavior of the diagonal elements, we take n ms to be even, and set ν s = + √ n ms σ ν for even s and ν s = − √ n ms σ ν for odd s. We assume that the distribution of the real g s,ij is Gaussian having the variance σ 2 = σ 2 g , with σ 2 g /σ ν = J/(2N) 3/2 . In this section we set σ ν = σ g = J/(2N) 3/2 for simplicity. As explained in Appendix B, if we identify s g s,ij g s,kl νs defined in this way with J ij,kl /(2N) 3/2 , the properties needed in the real-SYK model are satisfied at n ms = ∞. We collected samples by using independent real Gaussian random values of {g ij }. In Fig. 1, we plot the distribution of J ij,kl with {i, j} = {k, l} and 10 4 samples (a) and the diagonal elements J ij,ij with 10 5 samples (b). The distributions have different shapes for smaller values of n ms , but they quickly approach Gaussian distributions with corresponding variances for the real-SYK model as n ms increases. In Fig. 1(c), we plot the energy spectrum of this model using 10 4 samples with N = 10 and Q = N/2. The energy spectra become closer to that of the real-SYK model as n ms increases (for comparisons regarding other quantities, see Appendix A).

Creating the model
In this section, we explain how to create the model (11), a simplified version of the SYK model, in a system of optical lattices loaded with ultracold gases. We consider a two-dimensional gas of spin-polarized fermionic atoms confined in an optical lattice. In the proposed scheme, we utilize the PA process that coherently converts two atoms into a bosonic molecule in a certain electronic (or hyperfine), vibrational, and rotational state [37].
We assume that molecules are confined also by the optical-lattice lasers confining atoms.  However, since in general the lattice depth for molecules may be controlled independently from that for atoms, we assume that the former has the sign opposite to the latter. In this situation, the potential minima of the molecular optical lattice sit right next to those of the atomic optical lattice, as illustrated in Fig. 2(a), such that we do not have to take into account the effects of the onsite interactions between an atom and a molecule, which would otherwise complicate the levels of the atomic and molecular bands. We assume that the optical lattices are so deep that atoms and molecules in each lattice site are completely isolated. To make the manipulation of the system easier, we remove all the atoms in the lattice sites neighboring to occupied sites. We also assume that each occupied atomic lattice site contains Q atoms. We regard the band degrees of freedom in the atomic site as the physical site index of the SYK model. More specifically, the first, second, third, . . ., N-th bands correspond to i = 1, 2, 3, . . . , N sites. We write the energy of the lowest molecular band and that of the i-th atomic band as E m and E a,i .
Let us introduce a PA laser, which couples atomic bands i(≤ N) and j(≤ N) with the lowest molecular band. The frequency of the PA laser is chosen as where E (2) i,j = E a,i + E a,j , and ν denotes the detuning. We consider a situation in which all the combinations of the two atomic bands (i, j) are coupled via independent PA lasers as shown in Fig. 2(b). For such a situation to be possible, |ν| has to be larger than the linewidth of the PA lasers Γ PA and that of the molecular state Γ ms . In addition, the condition |ν| ≪ ∆ min has to be satisfied, where ∆ min denotes the minimum level spacing in E (2) i,j ≤ E where the atom-molecule coupling constant is given by Ω i,j (r) denotes the intensity of the PA laser while w m (r) and w a,i (r) represent the Wannier function of the 1st molecular band and the i-th atomic band. The absolute value of the detuning |ν| is assumed to be much smaller than the onsite interaction U between two molecules in order to avoid double occupancy of the molecules. For the same reason, U has to be much smaller than the minimum level spacing ∆ min in E i,j or sufficiently larger than the maximum level spacing ∆ max in E (2) i,j . Moreover, the level spacing between the first and second molecular bands ∆ MB is assumed to be larger than ∆ max in order to avoid accidental couplings between higher molecular bands and the atomic bands via the PA lasers.
A PA molecule has many vibrational and rotational states. When ∆ max <∆, we can extend the scheme described above to include couplings of two atoms with multiple molecular states, where∆ denotes the minimum level spacing of the involved molecular states. The extended system is now described by Eq. (10). When |ν s | ≫ |g s,ij |, one can integrate out the molecular degrees of freedom through the second-order perturbation theory with respect to the atom-molecule couplings, leading to the effective Hamiltonian of Eq. (11). Notice that the precise condition for the second-order perturbation theory to be valid is shown in Appendix C. We emphasize that the coupling constant g s,ij can be controlled independently with respect to indices s, i, and j because each coupling is created via an independent PA laser. Setting ν s to be bimodal, √ n ms σ ν for even s and − √ n ms σ ν for odd s, and the distribution function of g s,ij to be Gaussian with the variance σ 2 = σ 2 g , the coupling J ij,kl ≡ (2N) 3/2 s g s,ij g s,kl /ν s becomes Gaussian random for sufficiently large n ms as shown in the previous section.
It is useful to summarize the necessary conditions for this scheme to be valid in terms of the several relevant energy scales, Here, t i and τ exp denote the intersite hopping for atoms in the i-th band and the lifetime of the experimental system, which is approximately a few seconds in typical experiments of ultracold gases in optical lattices. Notice that the total number of necessary PA lasers in this scheme is n ms × N(N − 1)/2. In Appendix D, we discuss the feasibility of this proposed scheme by taking 6 Li and a double-well optical lattice [43,44] (see also Fig. 3), as specific choices of atomic species and lattice configuration.

Measuring observables
Once the SYK model is realized in optical-lattice systems, various observables can be measured. One of the most interesting signatures of a black hole formation is the fast scrambling quantified by OTOC functions [38], with which the chaotic nature of the system can be  the OTOC functions [45], We follow the protocol to explain how to measure the OTOC functions in our system in the specific case thatV =ĉ i andŴ =ĉ j , namely Since the SYK model is homogeneously random and has no meaningful distance, this correlation function takes only two different cases, namely the onsite case (i = j) and the offsite case (i = j). Hence, it is sufficient to show the cases that i, j ∈ {1, N}.
The protocol requires a control qubit interacting with the probed system [45]. We assume that a double well occupied by a single atom plays the role of a control qubit and regard the state in which the atom occupies the right (left) well as the |0 C (|1 C ) state of the qubit as shown in Fig. 4(a). We also assume that the species of the qubit atoms is different from that of the SYK atoms and that the optical lattice potential for the former can be controlled independently of that for the latter. We locate the qubit double well in such a way that its left well is well overlapped with the site for the SYK atoms [see Fig. 4(b)]. In this situation, the qubit atom has the onsite interactionsŨ i with the SYK atom in band i when the qubit state is |1 C . Specifically, supposing that the optical lattice for the SYK atoms is given by Eq. (D1), that for the qubit atoms may be formed by the following double-well optical lattice [43,44], whose spatial distribution for V ′ 0 < 0 and R ′ = 0.3 is depicted in Fig. 5. The protocol to measure C i,j (t) is summarized as follows: whereÎ denotes the identity matrix.X C andŶ C denote the x and y components of the Pauli matrices for the control qubit. Taking the offsite case that i = 1 and j = N, let us elaborate this protocol item by item. The onsite case can be treated in a very similar way. First, since apply a π-pulse with frequency ω = E a ′ − E a,1 −Ũ 1 , which resonantly couples the E a,1 and E a ′ states if the qubit state is |1 C . Here E a ′ denotes the energy of the lowest band of the atoms in another state. Hence, the application of the π-pulse leads to the operation ofĉ 1 for the |1 C state and that of the identity matrix for the |0 C state, thus creating process (ii).
We next turn off the atomic interaction and perform the unitary time evolution of the system until τ = t. This corresponds to process (iii). At τ = t, we apply a π-pulse with frequency ω = E a ′ − E a,N , which resonantly couples the E a,N and E a ′ states. This corresponds to the operation ofĉ N , namely process (iv).
Process (v) requires the sign of the Hamiltonian to be inverted. Such a manipulation can be made for our SYK model simply by inverting the sign of the detuning for all the PA lasers. We perform the unitary time evolution of the inverted Hamiltonian until τ = 2t. At τ = 2t, we set the atomic interaction to be repulsive (Ũ i > 0) and apply a π-pulse with frequency ω = E a ′ − E a,1 , which leads to the operation ofĉ 1 for the |0 C state and no operation for |1 C . This is nothing but process (vi).
In order to obtain X C , we need to measure the population of the bonding state in the qubit double well. Such a measurement can be made by means of the band-mapping techniques used in Ref. [44]. On the other hand, in order to obtain Ŷ C , we need to measure the current between the two wells, which is possible using the optical lattice microscope techniques [46][47][48][49][50]. Thus, process (vii) is feasible.
It has been also suggested that the degeneracy of the ground state in the SYK model can be read off from the low-temperature behavior of the single-particle Green's functions [31], We suggest that G i,j (t) can be measured in a way very similar to the one described above.
Specifically, skipping processes (iv) and (v) corresponds to a protocol to measure G i,i (t).
The offsite case (i = j) is also possible simply by replacingĉ i operation in process (vi) witĥ c j .

Bottlenecks and possible solutions
While we have shown in the previous sections that the SYK model can be realized in principle by using ultracold fermionic atoms in optical lattices coupled via PA lasers, there remain difficulties in practice, which have to be resolved in order to realize actual experiments. In this section, we describe such difficulties and discuss possible solutions to them.
12 First, the severest bottleneck is the number of necessary PA lasers, n ms × N(N − 1)/2. For instance, when n ms = 36 and N = 16, 4320 PA lasers are required. Implementation of lasers with as many frequencies as O(10 3 ) in a single experiment is difficult with current experimental technology. A possible way to circumvent this difficulty is as follows. The Gaussian randomness of the coupling J ij,kl in the SYK model, which requires the use of multiple molecular states leading to the factor of n ms in the number of necessary PA lasers, might be needed only to make the theory analytically solvable. In other words, the modified SYK model of Eq. (11) with only one or a few molecular states may exhibit the SY state at low temperatures. Indeed, in a supersymmetric generalization of the SYK model, which has been proposed very recently, the system exhibits the SY state even though the coupling J ij,kl is not Gaussian random [51]. In future theoretical studies, it will be important to examine the robustness of the SY state in the absence of the Gaussian randomness in more general situations without supersymmetry. If the SY state can survive at n ms = O(1), the number of necessary PA lasers for N 10 will be reduced to O(10 2 ). Even in this case, preparing such a number of PA lasers remains as a challenge.
The second bottleneck is that all of the multiple PA lasers have to have ultranarrow linewidth. This requirement stems from the condition (16), meaning that the linewidth has to be much smaller than the detuning |ν s |. As shown in Appendix D, if we choose the double-well optical lattice of Eq. (D1) and set N = 16, V 0 = −60E R , R = 0.59, and θ = π/6, then the minimum level spacing is given by ∆ min = h × 66.7 Hz such that Γ PA 2π × 1 Hz is required. State-of-art experiments have successfully stabilized a laser with a single frequency to the extent that Γ PA ∼ 2π × 0.1 Hz, aiming to its application to optical-lattice atomic clocks [52,55,56]. However, achieving such narrow linewidth for lasers with multiple frequencies is challenging for current experiments. An alternative route to circumvent this difficulty may be to design a new configuration of optical lattice optimized for enlarging ∆ min significantly compared to the case of the double-well optical lattice.
Third, it is unclear which molecular states are suited for our purpose because some information regarding molecular properties is currently unknown. Specifically, while it is more desirable to have stronger coupling between atomic and molecular states, i.e., larger Franck-Condon factor, we do not know which molecular states have relatively stronger coupling. Information regarding linewidths of electronically ground-state molecular states is insufficient as well. Moreover, in order to satisfy the conditions (18) and (19), we need to confirm that |U s,s ′ | is not too small by accident but currently we do not know the values of the swave scattering lengths determining the interaction between two molecules. These unknown properties can be revealed in a step-by-step manner with current experimental technology while it requires a lot of efforts.
Fourth, while we assumed that the optical-lattice depth for molecules can be controlled independently of that for atoms, such a situation has not been realized in experiments thus far. However, assuming that a PA molecule consists of two electronically ground-state atoms with different hyperfine states, at least one of the two atoms forming the molecule has a hyperfine state different from atoms in the SYK system. An optical lattice whose depth can be controlled independently of two hyperfine states has been already realized in experiments [57,58]. The optical lattice of this type should also allow for independent control of the lattice depths for the atoms and molecules.
Finally, as for the measurement scheme of the time-dependent correlation functions, the preparation of qubit atoms interacting with the probed system has never been realized thus far, while some other theoretical works recently proposed similar measurement schemes using control qubits [45,59,60]

Conclusion
We have suggested that ultracold gases in optical lattices can be applied to experimental studies of quantum gravity under the assumption of the holographic principle. As a specific example, we have proposed that creating the SYK model, whose low-temperature state has been conjectured to be dual to AdS 2 black holes [31], is in principle possible. We have shown how to measure the OTOC functions and the single-particle Green's function, which characterize the properties of the black hole, with use of a control qubit consisting of an atom in a double well. Moreover, we have discussed practical difficulties in realizing our proposal with current experimental technology, and how they might be circumvented. We emphasize that while our proposal to realize the SYK model in experiment is incomplete in a practical sense because of the remaining difficulties, it has made a first step towards the experimental realization of the SYK model. In this sense, the present work has significantly advanced our original idea that quantum gravity can be studied in optical-lattice systems loaded with ultracold gases with the help of the holographic principle.
We chose our specific example because it looked the simplest among the currently available models with holography. However, the SYK model is still rather complicated in the sense that it has an unnatural two-body hopping that has to be Gaussian random. Hence, it will be useful to explore its simplified variants or other quantum gauge models with holography that can be created more easily in optical-lattice experiments. We finally note that while the Hawking radiation is one of the most important issues regarding quantum black holes, whether the Hawking radiation can be seen in the SYK model is not clear at this stage.
Answering to this question will be an imperative task for future theoretical studies. At very least, it should be possible to study a variant of the information puzzle, associated to the way that the information is encoded in a black hole [67].

Appendices:
A Comparison of the SYK model with its variants In this appendix, we compare the three versions of Sachdev-Ye-Kitaev (SYK) model, namely the original (complex) one, the real one, and the modified one (11) to demonstrate that the last one is a quantitative approximation of the original one even at finite N and n ms . We set the reduced Planck's constant and the Boltzmann constant to be = k B = 1 We first compare the complex and real SYK models. In the complex and real-SYK models, normalized by dividing by N, are shown. The energy E is calculated as The disorder average E and Q are taken by using random couplings. From the plots, we can see a clear agreement at large N. As expected from the standard 1/N-counting, two theories give the same values of E and Q up to the sub-leading corrections of order N 0 . We have also calculated the entropy S defined by The entropy approaches S = N log 2 as the temperature T is increased, which is expected from the fact that all 2 N states can equally contribute in the high-temperature limit.  The real-time, same-site density-density correlation function i n i (t)n i (0) − n i (t) · n i (0) /N calculated for (a) N = 6 using 10 3 samples, and (b) N = 10 using 10 2 samples. The data for T = 1 and 0.1 for the complex and real SYK models are plotted together with the data for the model (11).
Note that log Z < log Z in general. The result is shown in Fig. A3. For T /J 1, S/N is already almost converged at N = 6, while for smaller T , S/N is an increasing function of N. Notice that the entropy of the complex SYK model at finite N has been computed also in Ref. [68].
As an example of a two-point function, in Fig. A4 we present the same-site density-density correlation function C nn (t), which is defined using the number operatorn i =ĉ † iĉ i by For higher T the obtained normalized entropy is linear in 1/N, however for T /J ≪ 1 the curve is more convex, which suggests that the actual N → ∞ limit may be significantly larger than the value plotted and may converge to a finite value as T → 0.
Here, the operatorÔ(t) at time t isÔ(t) = e iĤtÔ e −iĤt and the expectation values are calculated as in the cases of the charge and the energy: with |ψ (Q) j being the many-body wavefunction corresponding to the eigenenergy E (Q) j . Having corroborated the quantitative agreement between the complex and real SYK models, we next compare the modified SYK model of Eq. (11) with the real SYK model. As shown in Fig. A5 for Q and E , we observe that the results are already similar for n ms = 8 for N = 6 and 10. In Fig. A6 (a) we plot the entropy S as a function of the temperature for N = 6, 8, 12 for n ms = 16, along with the result of linear extrapolation to 1/N → 0. As in Fig. A6 (b), for T 0.1J the obtained entropy is almost linear in 1/N, however for lower temperatures the dependence of S on 1/N is more convex, indicating that the linear fit from N = 6, 8, 12 may be underestimating the value of S(N → ∞) at T → 0 and that S may converge to a finite value.
It is worth stressing the importance of large n ms . Even when n ms = 1, the model (11) looks like the real-SYK model if we identify J ij,kl /(2N) 3/2 with g ij g kl /ν, where g ij ≡ g 1,ij and ν ≡ ν 1 . However, the distribution of the latter is not Gaussian in general for a given distribution of {g ij } and even worse, i.e., the randomness is not strong enough; for example, when g 12 g 12 /ν and g 34 g 34 /ν are large, g 12 g 34 /ν is also large. Note also that g ij g ij is always positive.
B Properties of J ij,kl = (2N ) 3/2 n ms s=1 g s,ij g s,kl /ν s In this appendix, we explain basic properties of J ij,kl ≡ (2N) 3/2 nms s=1 g s,ij g s,kl νs . We take ν s = + √ n ms σ ν for even s and ν s = − √ n ms σ ν for odd s, and the Gaussian weight of g s,ij is chosen to be e −g 2 s,ij /(2σ 2 g ) √ 2πσg , with σ 2 g /σ ν = (2N) −3/2 J. It will turn out that this corresponds to the Gaussian random coupling J ij,kl needed for the real-SYK model, with J = 1. Generic values of J can be realized by rescaling g s,ij and/or ν s . Firstly let us show that the distribution of x ≡ J ij,kl converges to e −x 2 / √ π for (i, j) = (k, l) and e −x 2 /2 / √ 2π for (i, j) = (k, l). Then, we should show that, when real numbers which means the width is 1. We can also show J ij,kl J pq,rs ∝ (δ ip δ jq − δ iq δ jp )(δ kr δ ls − δ ks δ lr ) + (ij ↔ kl). Let us note that we only have to show J ij,kl J pq,rs ∝ δ ip δ jq δ kr δ ls + (ij ↔ kl) for i < j, k < l, p < q and r < s. It is equivalent to show that ( s g (s) Q /ν s ′ ) = 0 unless I = P, J = Q or I = Q, J = P , where indices I, J, P, Q represents (i, j), (k, l), (p, q) and (r, s). With this notation, If I = J or P = Q, we can rewrite it as It is a bit tricky to show JJJ = 0; actually it holds when n ms is infinity. For simplicity, let us suppose I = J, P = Q, V = W . Then, For the same reason, we have and so on.
We further note that, if we can introduce complex g s,ij , we may identify where the non-perturbative partĤ 0 and perturbative partV are given bŷ We see from Eq. (C2) that the non-perturbative energy depends only on the number of particles in each molecular states. This means that all different atomic configurations with no molecule are degenerate inĤ 0 . The non-perturbative energy of these degenerate states is given by E 0 = 0. We define the Hilbert subspace spanned by all these degenerate states with no molecule as D. We note that states with a molecule or two molecules appear as virtual states in the second-or fourth-order perturbation.
In order to derive the effective Hamiltonian, we perform the Schrieffer-Wolff transforma- whereP 0 is the projection operator on D.
To determine the transformation matrixŜ we require that in eŜĤ m e −Ŝ all the matrix elements connecting states in D with those outside of D are zero, i.e., that eŜĤ m e −Ŝ is block-diagonal [70]. While our main purpose is to derive Eq. (11) of the main text, which corresponds to the effective Hamiltonian up to the second-order perturbation, we here describe the terms up to the fourth order, in order to discuss the validity condition of the second-order approximation. Notice that the odd order terms do not exist in the effective Hamiltonian because inV all the matrix elements connecting two states with the same number of molecules are zero.
The second-and fourth-order terms can be formally written aŝ Substituting Eqs. (C2) and (C3) into Eqs. (C6) and (C7), we obtain (C13) It is obvious thatĤ (2) eff is equivalent to Eq. (11) of the main text. Since we have assumed that |ν s | ≪ |U ss ′ |, the first term in the right-hand side of Eq. (C10) is much smaller than the second term. Hence, we neglect the first term and compareĤ (2) eff with the second term inĤ (4) eff in the following discussions. We recall that g s,ij is assumed to be Gaussian random with the standard deviation σ = σ g [see Eq. (21) of the main text] and that ν s is assumed to be (−1) s √ n ms σ ν . Combining these assumption with Eq. (C11), we obtain The scale of the eigenenergies ofĤ (2) eff is set by 4N 3 K 2 ij,kl /3! while that ofĤ (4) eff is set by 8N 7 L 2 ii ′ ,jj ′ ,kk ′ ,ll ′ /7! [33]. In order for the second-order approximation to be valid, the former must be much larger than the latter. This condition implies that In Appendix D, we show that the condition of Eq. (C16) can be safely satisfied in a realistic situation.  In this appendix, in order to discuss the feasibility of our scheme for creating the SYK model, we consider the following optical lattice, which is Eq. (20) of the main text. This optical lattice consists of two square optical lattices and a represents the lattice spacing of the one with the shorter period. We assume that V 0 < 0 for atoms while V 0 > 0 for molecules. In Fig. 3 of the main text, we show the spatial profile of this potential for V 0 < 0, R = 0.59, and θ = π/6. Such an optical lattice is often used to create a double-well optical lattice [43,44], whose unit cell is a double well potential, and is advantageous for the proposed scheme in the sense that the band levels of the atomic site have no degeneracy as shown in Fig. D1, where the eigenenergies of a single atom in the optical lattice at zero quasi-momentum are plotted for V 0 = −60E R , R = 0.59, and θ = π/6. E R ≡ 2 π 2 2ma 2 denotes the recoil energy. If there are any degenerate levels, then ∆ min = 0 such that the condition (16) can not be satisfied.
To evaluate the energy scales appearing in the necessary conditions (15)- (19), let us specifically choose 6 Li as the fermionic atoms confined in our system. Use of this species in cold-atom experiments is rather standard. Setting the lattice spacing to be a standard value, namely d = 532 nm, leads to the recoil energy E R = h × 29.2 kHz. Taking the values of the parameters used in Fig. D1 and setting N = 16 immediately give max(t i ) ∼ h × 0.5 Hz, ∆ min = h × 66.7 Hz, and ∆ max = h × 1.96 MHz. If we set |ν s | = h × 10 Hz, n ms = 36, and σ g /σ ν = 0.3, then J = h × 27.2 Hz. Hence, the first condition (15) is safely satisfied. We note that with these values of N, |ν s |, n ms , and σ g /σ ν the condition (C16) is safely satisfied as well.
Since the second condition (16) requires the information of the linewidths, let us first estimate Γ ms . If a PA molecule consists of one electronically ground-state alkali atom and one electronically excited alkali atom, a typical scale of the linewidth is Γ ms ∼ 2π × 10 MHz [37], which is much larger than |ν s | and can not be used for the present scheme. In contrast, if a PA molecule consists of two electronically ground-state alkali atoms, Γ ms is much smaller in general. A coherent coupling to such an electronically ground-state molecule can be created with use of the two-photon Raman PA techniques [37]. For instance, in the case of 87 Rb atoms confined in an optical lattice, molecular states with linewidths as narrow as Γ ms ∼ 2π × 1 kHz have been observed [71]. In the case of 6 Li atoms, detailed experimental searches for linewidths of electronically ground-state molecules in optical lattices have not been performed. However, it is known at least that Feshbach molecules of 6 Li, which correspond to an electronically ground-state molecular state, can have lifetime as long as 10 s in the absence of an optical-lattice potential [72,73], meaning that its linewidth can be as narrow as Γ ms ∼ 2π × 0.1 Hz.
As for Γ PA , state-of-art experiments have developed lasers with ultra-narrow linewidth for application to optical-lattice atom clocks such that the linewidth can be as low as Γ PA ∼ 2π × 0.1 Hz [52,53,56]. Using the electronically ground-state molecules and the state-of-art lasers, the condition (16) can be overcome in principle. Notice, however, that implementation of such narrow linewidths for all PA lasers with many different frequencies has never been realized thus far.
Furthermore, we assume that V 0 = 2 × 10 5 E R for the molecular optical lattice so that ∆ MB = h × 10.9 MHz. Since the level spacing of the rotational states of a 6 Li 2 molecule is typically∆ ∼ h × 100 MHz, the third condition (17) is also satisfied. Finally, we estimate that |U s,s ′ | ∼ 3 MHz under the assumption that the s-wave scattering lengths between two molecules take a typical value |a s | ∼ 100a B , where a B denotes the Bohr radius. With this estimation of |U s,s ′ |, the fourth and fifth conditions (18) and (19) are satisfied. This means that it is in principle possible to create the modified SYK model (11) at least up to N = 16 by means of the proposed scheme with the specific choices of the optical lattice potential of Eq. (D1) and the atomic species of 6 Li.