A model of randomly-coupled Pauli spins

We construct a model of Pauli spin operators with all-to-all 4-local interactions by replacing Majorana fermions in the SYK model with spin operators. Equivalently, we replace fermions with hard-core bosons. We study this model numerically and compare the properties with those of the SYK model. We observe a striking quantitative coincidence between the spin model and the SYK model, which suggests that this spin model is strongly chaotic and, perhaps, can play some role in holography. We also discuss the path-integral approach with multi-local fields and the possibility of quantum simulations. This model may be an interesting target for quantum simulations because Pauli spins are easier to implement than fermions on qubit-based quantum devices.


Introduction
The Sachdev-Ye-Kitaev (SYK) model has been intensively studied for various motivations, ranging from condensed matter physics to quantum gravity via holography.Given the importance of the SYK model, it is natural to try quantum simulations.Indeed, there are a few attempts [1,2].Still, it is difficult to simulate the SYK model on quantum devices without some simplifications.One of the obstacles is that elementary degrees of freedom are fermions and fermions are non local when mapped to qubits.Specifically, via the Jordan-Wigner transform, Majorana fermions χa (a = 1, 2, • • • , N Maj ) satisfying { χa , χb } = 2δ ab are written in terms of Pauli strings (tensor products of Pauli matrices) acting on N spin = N Maj 2 spins as Here, we used Pauli matrices σx = 0 1 1 0 and the identity matrix which act on the local Hilbert space of each qubit.These long chains of Pauli matrices, which are as long as the number of degrees of freedom in the system, require a lot of resources (quantum operations) in digital quantum simulations.
In this paper, we will consider a spin model which is obtained by replacing all Pauli σz operators in (1) with the identity Îs.Such a theory contains only SU(2) spin variables (Pauli matrices σx and σy ) on N spin sites.For brevity, we will denote this model as 'SpinXY4' in this paper.XY refers to σx and σy and 4 refers to the number of Pauli operators in the interaction.
There are several reasons we are interested in such a model.First of all, this model can be studied more easily on quantum computers.Therefore, if this model inherits some interesting features of the SYK model, it will be an interesting target for quantum simulation in the near future and, hopefully, serve as a good starting point for the experimental study of quantum gravity via holography [3,4,5].We could hope that much of the physics is preserved by the replacement of fermions with spins, given that the Sachdev-Ye model (SY model) [6], which is closely related to the SYK model, is a model consisting of SU(M ) spin variables.A potential advantage of this model over the SY model is that there is only one limit (N spin → ∞) while the SY model requires the large-spin limit (M → ∞) and the many-spin limit (N spin → ∞).The simple structure in terms of spin-1/2 variables makes the simulation on qubit-based quantum devices straightforward. 1 Note that an important motivation for the large-M limit in the SY model is to avoid the spin-glass phase, and hence, we would like to know if a spin-glass phase appears in SpinXY4.
Our findings are the following: • We studied the density of states (DoS) up to N Maj = 2N spin = 34 by exact numerical diagonalization, collecting many samples with different random couplings.For small N spin , the DoS is almost indistinguishable from the one for the SYK model.As N spin increases, we see a small discrepancy near the edge, although the bulk of the spectrum looks very similar to SYK. (Sec.3) • Statistical properties of the energy spectrum are consistent with those of Random Matrix Theory (RMT), suggesting the absence of the spin-glass phase except for a few low-energy modes.(Sec.4) • The spectral form factor (SFF) has a long ramp that suggests a strongly-chaotic nature, similar to the SYK model.(Sec.4) • For some values of N spin , some correlation functions are quantitatively close to the counterparts in SYK at any time scale.(Sec.4.3.2 and Sec. 6) • While the Edwards-Anderson (EA) parameter defined using the σz operators as well as a generalized version of the EA parameter decrease monotonically as a function of the system size for a majority of the energy spectrum, their increase suggests that a small number of low-energy states behave as in the spin-glass state.(Sec.5) We believe these findings provide us with good motivation for further investigations.This paper is organized as follows.In Sec. 2 we give the precise definition of the model.We also provide an incomplete list of potential generalizations.In Sec. 3 we study the density of states.We make a quantitative comparison with the SYK model and find an intriguing resemblance, except for the edges.In Sec. 4 we study the correlation of energy eigenvalues and compare it with that of Random Matrix Theory.We observe striking similarities with the SYK model: agreement with RMT is observed except for a small number of low-lying modes, and the agreement extends to a wide energy band (equivalently, a long ramp is observed in the spectral form factor).In Sec. 6 we study two-point functions.The late-time behavior is consistent with RMT and similarities with the SYK model are observed.For certain choices of operators and values of N , we observe a striking quantitative coincidence at all time scales.In Sec. 5, we introduce a generalized version of the EA parameter, defined between eigenstates belonging to the two parity sectors, and study it along with the EA parameter.In Sec. 6, we study the two-point functions as a function of time and find strong similarities with the SYK model.In Sec. 7 we introduce a pathintegral formulation for the description of the large-N limit of these spin models based on collective multi-local fields following a closed set of Schwinger-Dyson equations.We end the paper by commenting on the implementation of a Trotterized Hamiltonian evolution of the SpinXY4 model on a quantum device in Sec. 8 and then we conclude with an outlook.
Note added: We note that SYK-like behavior has been observed in the spectral function of random Heisenberg magnets for low and finite frequencies [7], where the ground state is spin glass [8].Also see [9] for saddle-point equation study and numerical study for the SpinXYZq model (see below for definition), which was initially introduced as the quantum p-spin glass model [10].Possibility of studying black hole spacetimes by a spin system with another type of four-spin random couplings has been previously discussed in [11,12].

Definition of the model
spins instead of N Maj Majorana fermions.Let Ôa be the counterpart of χa , i.e., Ôa is obtained by replacing σz with Î in χa .Specifically, Ô2j−1 = σj,x and Ô2j = σj,y , where The Hamiltonian of the model is the following: in which the couplings J abcd are chosen from the standard normal distribution and η abcd is the number of spins whose both x and y components appear in (a, b, c, d), e.g., η 1357 = 0, η 1235 = 1, η 1234 = 2.We need i η abcd for the Hermiticity of the Hamiltonian. 2We will compare this model with the SYK model with q = 4 which we rename for conciseness as 'SYK4': We chose the normalization of the random couplings J abcd in such a way that the large-N Maj limit of the SYK model simplifies.Specifically, the energy E and entropy S scale as N 1 Maj when the temperature T is fixed to be an order-N 0 Maj value and characteristic time scales, such as the decay rate of a two-point function, are of order N 0 Maj .Despite an apparent similarity at a formal level, the Hamiltonians (5) and (7) are clearly different because we are using different building blocks: Pauli spins Ô in the former and fermions χ in the latter.We could interpret Ô2a−1 ± i Ô2a as the creation and annihilation operators of a hard-core boson3 rather than a fermion.

Parity
A convenient basis of the Hilbert space is { s 1 , s 2 , • • • , s N spin }, where s a = ±1 (a = 1, 2, • • • , N spin ) represents a spin up or spin down at each site.Because σ a,x and σ a,y change s a to −s a (up to down), and because the Hamiltonian is a sum of products of four of them, the product N spin a=1 s a is conserved.We can see this also by noticing that Ĥ and Γ ≡ σ1,z ⊗ σ2,z ⊗ • • • ⊗ σN spin ,z commute.Therefore, the Hamiltonian can be written in a block-diagonal form with two blocks consisting of γ ≡ We will call γ = +1 the parity even sector and γ = −1 the parity odd sector.They correspond to the parity-even and odd sectors in the SYK model.

Possible variants
Similarly to the case of the SYK model, we can consider many variants of the SpinXY4 model.

q-local models (SpinXYq)
We can take the number of spins in each interaction term to be a generic number q: where the standard choice of the normalization factor is N = q!(N Maj − q)!/N Maj !. 4Note that q can be odd, in which case parity is not conserved. 5

Adding or removing σ z
In the SpinXY4 model defined above we allowed σ x and σ y on the same site to appear in the same interaction term.We can forbid this to happen and this amounts to setting η = 0.Such a modification should not change the theory in the limit of N spin → ∞.However, there are some differences that are not captured in the large-N spin limit.For example, the universality class (from RMT) is the Gaussian unitary ensemble (GUE) for any N spin when σ x and σ y are allowed at the same site, while it is the Gaussian orthogonal ensemble (GOE) for even N spin and GUE for odd N spin when σ x and σ y are not allowed at the same site. 6 We could also consider the random q-local coupling of σ a,x , σ a,y , and σ a,z with a = 1, 2, • • • , N spin .Such a model could be called SpinXYZq.The density of states for this model has been studied in Ref. [10].

Complex model
The analog of complex fermions are Ô2a−1 ± i Ô2a .By using them, we can define the analog of the complex SYK model.

Coupled SYK-like models
We can also prepare multiple copies of the SpinXY model and couple them.A particularly interesting model of this kind would be the analog of the coupled SYK model [18] which could be used to study the traversable wormhole [19].Note that the traversable wormhole is a promising target of experimental quantum gravity via holography [20], and there has been an attempt to study the SYK model on a quantum device in this context [2].

Qudit models
We can define a model replacing Pauli operators with spin-s representations with s > 1 2 and correspondingly replace qubits with qudits as the fundamental quantum registers for quantum simulations.

Density of states
In this section, we define the density of states (DoS) by taking the average over many samples with different random couplings J abcd .Practically, we introduce a binning separation of the energy spectrum and count the number of energy levels in each bin.When we combine many samples, we can take a very fine binning width (due to the large statistics of counts).In Fig. 1 we show the DoS obtained in this way.We can see similar shapes across different system sizes.In Fig. 2, we compared SpinXY4 and SYK4 in the same panel.The two densities are almost indistinguishable, except for a tiny discrepancy near the edges.
Note that we did not separate the two parity sectors to obtain these results.Whether we separate or not the two parity sectors, we see almost identical densities. 6To see this, it is convenient to perform a unitary transformation that maps σa,y and σa,z to σa,z and −σ a,y The normalized density of states for SpinXY4 (left) and SYK4 (right).The contributions of the two parity sectors are not separated.The number of samples is 2 ρ(E)

Edge of the energy spectrum
Let us look closely at the edge of the spectrum, where small deviations between the two models are apparent.Fig. 3 is a zoomed-in view of the lower edge of the DoS from N Maj = 2N spin = 16 to 34.The horizontal axis is E/|⟨E 0 ⟩ SYK |.As N Maj increases we see a small but clear discrepancy between SpinXY4 and SYK4.
For the SYK model, the DoS behaves as ρ(E) ≃ A sinh(B (E − C)) near the lower edge, where C = ⟨E 0 ⟩, and A, B, C were estimated analytically [21,22,23].A natural question is whether the SpinXY4 shows a similar pattern.A nontrivial technical issue here is that the smallest eigenvalue tends to have a large fluctuation at finite N .To deal with this issue, we consider the distribution of [23] for each parity sector.Note that E 0 is subtracted in a sample-by-sample fashion.This option could remove sample-by-sample fluctuations of E 0 .The distribution of E ′ is more relevant than that of E when we consider the low-temperature region with the quenched averaging.In Fig. 4, we plotted the density of E ′ for SpinXY4 and SYK4.While we can see sharp edges for both models, the discrepancy grows as N increases.In Fig. 5, we tried to fit the density of Although this fit ansatz is not bad for N spin = 12 and 13, we do not find a nice fit for N spin = 14 and 15.N spin = 15 SpinXY4 0.00539588 sinh(6.8447E'

Level correlations
In this section, we compare the correlation in the energy spectrum with that of Random Matrix Theory (RMT).We will study two sectors corresponding to γ ≡ N spin a=1 s a = ±1 separately.Unlike SYK4, we do not find eigenvalue degeneracy within each sector nor between the two sectors.We do observe agreement with RMT except for a small number of low-lying eigenvalues.Such an agreement suggests that this model is ergodic rather than in a spin-glass phase.(See e.g., Refs.[24,25] for the spectral analysis of spin glass.)As we will see in Sec.4.3, the spectral form factor of our model resembles that of SYK4.This implies a very strongly chaotic nature of the model.

Nearest-neighbor level spacing
To compute the eigenvalues, we can utilize the block-diagonal structure of the Hamiltonian, i.e., we can diagonalize 2 N spin −1 ×2 N spin −1 blocks corresponding to γ = ±1 separately.In each sector, we sorted energy eigenvalues in increasing order as The nearest-neighbor level spacing is defined by s i ≡ E i+1 − E i .To compare it with RMT, we need to unfold the spectrum.Here we use the fixed-i unfolding [26], i.e., we define the unfolded spacing si by si = s i /⟨s i ⟩ J for each i.
In Fig. 6, the distribution of the unfolded level spacing P (s i ) is plotted for several values of i.For N spin .For N spin = 11, although a significant difference from RMT can be seen only for i = 1, we see almost no difference from RMT at i ≥ 2. For N spin = 15, we see a larger deviation from RMT at small i.However, the agreement with RMT is not bad already at i = 4 and it is hard to see a difference from RMT at i ≥ 10.For the SYK4, a good agreement with RMT is observed even for i = 1 [27].

Neighboring gap ratio
By using the unfolded level spacing si , we define the neighboring gap ratio r i as In the left panel of Fig. 7, we plotted Good agreement with RMT (the GUE universality class) [28] is observed at i ≥ 4. In the right panel, the same quantities for SYK4 are plotted for the values of N corresponding to GUE.Again, a good agreement with RMT is observed at i ≥ 4.
Note that the gap ratio can be sensitive to the unfolding near the edges of the energy spectrum.By using the unfolded level spacings, we can see good agreement even near the edge of the spectrum.Only the parity-even sector was used.In SpinXY4, the eigenvalues in the parity-even and odd sectors are not correlated but the plots for the two sectors are indistinguishable.In SYK4, parity-even and parity-odd sectors have the same eigenvalues when N spin = N Maj /2 is odd.SpinXY4 SYK4

Spectral form factor
A convenient quantity to see the correlation of energy eigenvalues in a wider energy band is the spectral form factor (SFF).The SFF can be defined for each parity sector as where Here the sum over states j is taken in γ = +1 or γ = −1 sector.The SFF starts with 1 at t = 0 and shows the slope, dip, ramp, and plateau.The ramp and plateau are universal among chaotic systems.If the ramp is longer (equivalently, if the onset of the ramp is earlier), the energy spectrum agrees with RMT in the wider energy band.We plotted g(t, β = 0) for our model in Fig. 8.We can compare it with the same quantity for the SYK model.We can see similar long ramps.
The onset of the ramp can be hidden by the slope.To see the onset of the ramp more accurately, a modified spectral form factor h(α, t, β) defined by [29,30] where is useful.By tuning a parameter α appropriately, the slope can fall much more quickly and the hidden part of the ramp can be revealed.h γ=±1 (α = 1, t, β = 0) is plotted in Fig. 9.

Spectral Form Factors with and without separating parity sectors
In the above, we defined the SFF for each parity sector.We could also combine two sectors.Specifically, by using In Fig. 10, we plotted g full (t, β) and g γ=±1 (t, β).We can see that at an early time and at late time.It can be understood as follows.[Left] g full (t, β) = g γ=±1 (t, β) can be seen at an early time.
• To see the late-time behavior, we use . The first two terms on the right-hand side are described by RMT with matrix size 2 N spin −1 × 2 N spin −1 .The third term vanishes because there is no correlation between two parity sectors and hence Z γ=±1 fluctuates around zero independently, averaging to zero: ⟨Z γ=+1 Z * γ=−1 ⟩ J = 0. Therefore, the late-time behavior of ⟨|Z full |⟩ J coincides with that of RMT.The factor-2 difference in (15) is explained by the difference in the dimension of the full and parity-fixed Hilbert spaces.

Comparison with SYK at all time scales
As we have seen above, the late-time features of the SFF capture the fine-grained energylevel correlations.On the other hand, at the early time, the SFF is sensitive to the density of states.Therefore, the observation so far indicates that the SFF of SpinXY4 resembles that of SYK4 closely both at an early time and at a late time.Now we would like to ask if the similarity can be observed at all time scales.
For quantitative agreement at the late time, we choose N such that N spin = N Maj /2 is odd because then both SpinXY4 and SYK4 are in the GUE universality class and hence we can expect the precise agreement at a late time.With such a choice of N , there is a two-fold degeneracy in the energy eigenvalues in SYK4.Therefore, there are 2 N spin −1 independent eigenvalues used in the SFF.If we keep only one of the parity sectors in the spin model, the numbers of eigenvalues match.Therefore, we compare g γ=±1 (t, β) in SpinXY4 and g(t, β) in SYK4.
In Fig. 11, we plot the spectral form factor for N Maj = 2N spin = 26 and 30, β = 0, 1,  and 2. In addition to g(t, β), we plot the 'connected part' defined by The agreement is strikingly good, although a small discrepancy is visible around the dip.

Edwards-Anderson parameter
The Edwards-Anderson parameter [31] is a standard tool to see if a given system has a spinglass phase or not.In this section, we study the generalized Edwards-Anderson parameter q gEA (j), here defined for the j-th lowest energy normalized eigenstates |E j ⟩ (E),(O) as [9] q gEA (j) = 1 Note that we do not include α = z in the sum because ψ (O) j σi,z ψ (E) j = 0 due to the parity conservation.We also study q zEA defined by Numerically, we observed that q zEA takes a nonzero value only for odd N spin .In Fig. 12 we plot the value of q gEA (j) as a function of the eigenstate index j ∈ [1, 2 N spin −1 ].For clarity we only plot the results for even N spin , however the results for odd N spin qualitatively agree with the even N spin case as we see below.Due to the symmetry concerning the overall sign of the Hamiltonian, the distributions of q gEA (j) and q gEA (2 N spin −1 + 1 − j) are identical.At N spin ≤ 14 we observed the following pattern: • q gEA (j) at small and fixed j increases as a function of N spin , indicating that the lowest energy eigenstates behave as spin glass states.
• For j > O(10 1 ), q gEA (j) decreases as a function of N spin .However, it is possible that q gEA eventually increase with N spin at any fixed j, if N spin becomes sufficiently large.
• At fixed N spin , q gEA (j) shows a power-law decay as a function of the eigenstate index, until j reaches ≈ 2 N spin −2 .
• The smallest value of q gEA decreases exponentially as N spin is increased.
Removal of terms with η abcd > 0, where multiple operators acting on the same spin are chosen, affects the value of q gEA (j) slightly but does not qualitatively change the pattern above.In Fig. 13, we plot q gEA and q zEA for odd N spin .We observe that q gEA behaves similarly to the case of even values of N spin .While q zEA (j) shows a similar behavior with a smaller number of j showing increase with N spin between N spin = 11 and N spin = 13, it decreases exponentially as N spin is increased for all j when terms with η abcd > 0 are removed from the model (5).
In summary, q gEA (j) suggests some low-energy states are in the spin-glass phase, although q zEA (j) suggests the opposite may be the case if terms with η abcd > 0 are suppressed.More studies will be needed to have a conclusive statement.In this context, we note that Ref. [9] studied the Edwards-Anderson parameter for the ground state of the SpinXYZq model.They observed a slow decline of the Edwards-Anderson parameter that is consistent with the absence of the glassiness, although the signal could mean that q = 4 is sitting at the border between spin-glass and ergodic cases in the sense that q = 3 and q = 5 respectively exhibit clear growth and decline of the Edwards-Anderson parameter.Also, see Ref. [32] that analyzed an analogue of the EA parameter to study the phase structure of the SYK model.

Two-point function
We consider the two-point function Note that we take the sum over all the energy eigenstates from both parity γ = ±1 sectors.We will take the average over random couplings separately for the numerator and denominator.Furthermore, we take the average over a = 1, • • • , N Maj .Here we consider the annealed average: We will also consider If the system is chaotic, late-time behaviors of such correlation functions should be understood based on RMT.We can repeat the argument for the SYK model [21] without a substantial change.
As for G x,y (t), the operators Ôa connect states with different parity, and hence, ⟨E| Ôa |E ′ ⟩ 2 is nonzero when |E⟩ and |E ′ ⟩ are in different parity sectors.Other than that, it can be approximated by a smooth function of E −E ′ , as suggested by the eigenstate thermalization hypothesis (ETH).As far as the late-time behaviors are concerned, we can approximate it with a constant.Then, contributions from two sectors with different parity, which are not correlated, cancel out and we do not see the ramp and plateau.As we can see in Fig. 14, this is indeed the case.We can see a close similarity with two-point function of ψ i (t) and ψ i (0) in SYK4 with N Maj = 2N spin = 16, 20 and 24 [21].
As for G z (t), the operators σj,z do not change the parity.Therefore, the late-time behavior resembles the sum of SFF in two parity sectors, and hence, we expect the ramp and plateau.In Fig. 15, we do see such a pattern.

Comparison with SYK model at all time scales
In Sec.4.3.2,we observed that the spectral form factors from SYK4 and SpinXY4 can be close at all time scales.Let us see if a similar coincidence can be seen for the two-point functions.
We take N spin odd so that both models are in the GUE universality class.The eigenvalues in two parity sectors in SpinXY4 are not correlated while the eigenvalues in the two  Note that Im G z (t) = 0 holds for β = 0, which is numerically confirmed.1024 samples are used, and the average over all operators and samples is taken before the absolute value is computed.10,14,18,22,26).Note that Im G z (t) = 0 holds for β = 0, which is numerically confirmed.1024 samples are used, and the average over all operators and samples is taken before the absolute value is computed.
parity sectors of the SYK model are paired.Therefore, we compare fixed-parity sectors in SpinXY4 and SYK4, and we choose the operators that do not mix different parity sectors.Specifically, we study the two-point function of σa,z = −i Ô2a−1 Ô2a = −i χ2a−1 χ2a , which is G z (t) defined by (21).
The results are shown in Fig. 16, for β = 0, 2. Overall, we find them remarkably similar to each other.For β = 0, we observe good agreement at early and late times, although some discrepancy is visible in between.For β = 2, we can see a small difference at late time as well.

Path-integral approach
We discuss how large-N spin systems can be studied systematically with path integral methods.To develop a systematic large-N (and 1/N ) expansion, the following features are needed: i) An invariant (collective) set of variables Φ(a) needs to be identified, generally as singlets under a U(N ), O(N ), Sp(2N ) or S N group operating on the system.
ii) A closed set of Schwinger-Dyson (SD) equations needs to be deduced and/or iii) a collective action describing the 1/N dynamics of collective variables Φ(a) needs to be established.
We note that items ii) and iii) should have an equivalent description in terms of Feynman diagrams, e.g., planar diagrams in matrix models, and bubble diagrams in vector models.
For theories of spin degrees of freedom, however, none of the required features items i) to iii) were obvious so far.We will see shortly that for spin systems the relevant symmetry group is S N , which induces an infinite set of collective variables.This S N symmetry is featured by expanding the kinetic term in the Lagrangian, thus our formalism described below applies to all spin systems.Let S a i = 1 2 σ i,a denote the spin operator, where a = x, y, z, and i = 1, . . ., N spin denotes the site index.From here on, we use the letter N instead of N spin , i.e., N = N spin .Furthermore, we take the variance of the random coupling to be J 2 .(Previously, we took J = 1.)Up to the 1/N -suppressed terms, the real-time path integral after the disorder average is given by where s = 1 2 and S ± i ≡ S x i ± iS y i .In the above formula we only consider the η = 0 sector in the Hamiltonian (5) since in the large N limit the η > 0 sector is of lower order in 1/N and hence can be dropped.To be more specific, one can see from above that the η = 0 sector gives a potential of order O(N ).On the other hand, one can show that the η = 1 and η = 2 sectors are of order O(1) and O(1/N ) respectively, and thus are suppressed in the large N limit.The kinetic term in the Lagrangian is This term does not have the O(N ) symmetry as opposed to the SYK model [33], while the S N symmetry is manifest.We can use the constraint to write S z = √ s 2 − S + S − and expand the denominator into a Taylor series as This expression motivates us to use the S N -singlet multi-time collective variables with L being the length of the sequence (number of pairs of S + i S − i ).The multi-time labels of these collective variables are themselves identical under the S L -exchange and we can therefore consider time-ordered sequences The set of collective fields extends the bi-locals operational in the O(N ) symmetry case.
The infinite sequence of multi-time collective variables will be shown to close under SD equations, giving a basis for the large-N limit of these spin-chain models.We note that these S N invariants are due to the kinetic terms in the Lagrangian, as such their appearance is a universal feature for spin systems.
In the strong coupling limit (1/J → 0), the potential term is dominant in the action (22).Then, O(N ) symmetry emerges and the bi-local description applies.Since it is analogous to the bosonic SYK model, and has been shown by [34] that the replica non-diagonal configuration is of lower energy, we will consider the quenched averaging which involves n replica fields.To compare with the well-known results in the SYK model, we will work in the Euclidean time τ = it.In the 1/J → 0 limit, we see that after rescaling the kinetic term in the action drops out, such that the replica representation of the partition function with the disordered average is where a is the replica index, and the Euclidean action is with an emerging O(N ) symmetry S i → O ij S j .This allows a bi-local as the invariant collective field: where we use the variable X to package the time variable t and the replica index a [33].
The bi-local field is symmetric under the exchange X ↔ Y : By contrast, the bi-local field in the SYK model is anti-symmetric [33].Thus, the partition function can be written as The Jacobian J [ϕ] is A standard way to deal with this is to introduce an auxiliary field, integrate out S ± i , and then eliminate the auxiliary field by solving the saddle-point equations.The end result is With this Jacobian, we have the collective action in the strong coupling limit In contrast to the SYK model [33], the coefficient in front of the Jacobian term is minus instead of plus.As mentioned before, for the low temperature (large β) limit, replica indices should be added with possible replica non-diagonal solutions [34].The SD equation is δA col [ϕ]/δϕ = 0, giving the relation We now consider finite coupling with only the S N symmetry and show the explicit form of the associated collective and SD equations.The general collective scheme for specifying the Jacobian J [Φ] applies [35].It represents a change to the invariants Φ({t} L , {t ′ } L ) defined in (25), with L = 1, 2, . . . .On the right-hand side of (25), the sum over i is analogous to the trace in the matrix model.This sequence of 'single-trace' fields is analogous to the 'loop' or 'word' variables of matrix models.Hence the basic building blocks will be the 'splitting' and 'joining' of 'single-trace' fields.The 'splitting' operation is resulting in a sum of variables of length L − 1.By using a L and b L−1 to mean ({t} L , {t ′ } L ) and ({t} L−1 , {t ′ } L−1 ), this operation can be written schematically as Note that the counterpart of this operation in the matrix model splits a loop into two loops.The 'joining' is where the l-th and k-th spins are taken out of the sequence.This is then a linear combination of traces of length L + K − 1, or schematically Often, it is hard to obtain the Jacobian explicitly, while it is not hard to determine ω and Ω as illustrated above.Still, we can write the saddle-point equation explicitly without knowing the Jacobian [35]: This is the large-N SD equation written explicitly in terms of the collective variables.As demonstrated in [35], this general formula applies to the O(N ) vector model, U(N ) Yang-Mills gauge theory, etc.It applies to the large-N spin systems as well, whose relevant collective variables are S N singlets, and one needs to substitute (38) and ( 40) into this formula.This set of equations represents a natural multi-time generalization of bi-local SD equations.It offers a possibility to search for more general ground state configuration of relevance at small temperatures.As a concrete example, we may apply these equations explicitly to the strong coupling limit.For simplicity let us assume that in this case we can have the replica-diagonal solutions such that we can ignore the replica indices.The action (29) can be written in terms of the Φ 1 : Since the action A[Φ] only depends on Φ 1 in the strong coupling limit, we see that the SD equations ( 41) reduce to giving where ϕ(τ, τ ′ ) = Φ 1 (τ ; τ ′ ) + Φ 1 (τ ′ ; τ ).Explicitly, for L = 1, 2, we have We see that the equation for L = 1 (45) is consistent with the saddle point equation of the collective action (36) we derived before.These equations have the recursive pattern that Φ L is determined by the Φ L−1 and Φ 1 .Let Φ L and ϕ 0 be the solution of the L = 1 part.Then, the following ansatz solves the above Schwinger-Dyson equations: Thus, all multi-local fields are determined solely by Φ 1 , consistent with that in the strong coupling limit the only degree of freedom is the bi-local field.

Toward quantum simulation
We have already discussed in the introduction how the SYK4 model requires long chains of Pauli matrices when embedding the Majorana fermions on qubit degrees of freedom (e.g. using a Jordan-Wigner transformation).Those Pauli strings have a length that grows linearly with the size of the system N spin , making it prohibitively challenging to approach the many-spin (N spin → ∞) limit.On the other hand, the advantage of SpinXY4 over SYK4 is that each term in the Hamiltonian involves at most only four qubits, regardless of the size of the system.A review of the computational resources for the quantum digital simulation of the SYK4 model can be found in Ref. [36,37] and a recent experimental trial for N = 6 Majorana fermions on a superconducting qubit device has been reported in Ref. [38].
As an example of what building blocks are required for the digital quantum simulation of the dynamics of SpinXY4, we focus on a first-order Suzuki-Trotter decomposition and reduce the simulation to a product of 4-qubit unitary operations.We can think of considering only spin operators acting on 4 different spins.Practically, if Û ≡ e −iJδtσ 1,x σ2,x σ3,x σ4,x can be realized for Jδt ≪ 1, the Hamiltonian time evolution can be coded into a circuit using native single-qubit and two-qubit quantum gates.We restrict to this exponential of a Pauli string because site indices can be handled by swapping qubit labels, and it is straightforward to replace σj,x with σj,y (or σj,z ) by a change of basis realized with singlequbit gates.Let us note that having Pauli strings with several terms that are exponentiated is a very common occurrence in quantum chemistry applications [39], such as in the Unitary Coupled-Cluster ansatz, and there exist numerous techniques to synthesize the corresponding quantum circuits, such as those based on phase gadgets and ZX-calculus [40].
As an example, the unitary operator Û above can be applied on 4 qubits using 6 CNOT gates in a staircase pattern, sandwiching a single-qubit Z rotation R z (α) = e − 1 2 iασz with angle α = J • δt, while the Hadamard gate H is used at the beginning and at the end of the circuit: The Hamiltonian (5) can contain terms that are acting on 2 qubits, 3 qubits, or 4 qubits at most.These terms will involve in general all qubits in the system, and all qubits will eventually be connected to all other qubits.For the purpose of Trotterized digital quantum simulations, a system of qubits arranged with a local geometry will require a large number of SWAP gates to implement all the interactions.On the other hand, we can expect that trapped-ion devices, such as Quantinuum H-series systems [41], can tame the non-local nature of the interaction.In the case of the quantum charged-coupled device architecture of H-series [42], qubits are realized by ions that can physically move on the device, effectively implementing all-to-all connections with no additional gate overhead [43].One additional feature of the Quantinuum H-series systems is the native 2-qubit gate ZZPhase(α), which implements directly the operator e − 1 2 iα(σ i,z ⊗σ j,z ) between any pair of qubits i and j with an infidelity that is proportional to the angle α, and around 0.5 − 2.0 × 10 −3 [44].Using such arbitrary-angle two-qubit gate we can express the circuit for Û above with one less 2-qubit gate, replacing the Z rotation and the neighboring CNOT gates by a single ZZPhase: R z (α) Overall, when taking into account the large number of terms in the Hamiltonian (5), this results in a great reduction of the total circuit depth, making the circuit for the Trotterized simulation less susceptible to noise [44].In recent demonstrations of quantum algorithms on Quantinuum H-series devices, circuits with a number of 2-qubit gates between 600 and 1000 were run without significant loss of signal [45] making use of tailored error detection techniques [46].Moreover, in the application of quantum optimization algorithms, a recent paper has implemented circuits with e −iθσ 1,z σ2,z σ3,z σ4,z Hamiltonian terms on Quantinuum H-series devices with up to 1000 2-qubit gates [47], using an optimization algorithm to reduce the number of gates by arranging Hamiltonian terms.The possibility of exploring the SpinXY4 variants described in Sec.2.2, such as introducing σz or reducing the number of terms to sparsify the interactions, using digital quantum simulations on real hardware is therefore a near-term challenge we would like to pursue in the future.

Conclusion and discussion
In this paper, we defined and studied the randomly coupled spin model (SpinXY4) by replacing Majorana fermions in the SYK model (SYK4) with Pauli spin operators.We found striking similarities between this model and the SYK model.We conclude that this is an interesting model of quantum chaos that can be simulated more easily on quantum computers.
There are many directions to be explored.It would be nice if we could solve this model or some variants analytically.For the SYK model, the effective action in terms of bi-local fields provided us with a better understanding of the model itself and its relation to gravity.
Hence, if we could do a similar analysis in terms of multi-local fields, we could understand if this model has a connection to gravity via holography.It would also be interesting to study various variants of the model including those suggested in Sec.2.2.We might be able to find an even simpler target for quantum simulation, or we might be able to find good models for holography or condensed matter physics.
As a final remark, we point out the similarity between the SpinXY4 Hamiltonian and the interactions in the matrix model for quantum black hole (see Refs. [4,5]).The matrix model contains several N × N matrices consisting of N 2 bosonic degrees of freedom.The interaction part of the Hamiltonian consists of O(N 4 ) 4-local terms of these bosons.In the coordinate basis truncation, each bosonic operator can be written as a sum of σz , and hence the entire interaction consists of the sum of 4-local interactions of σz s.For this reason, the quantum simulation of SpinXY4 may be a good starting point for the simulation of the matrix model.Furthermore, the Yang-Mills theory can be embedded into the matrix model [48,49], and hence, the same technique can be used to study Yang-Mills theory, and probably, the standard model of particle physics.

Figure 2 :
Figure2: Density of states from N Maj = 2N spin = 16 to 34.We can see that SpinXY4 and SYK4 have almost the same distribution except for a small discrepancy near the edges.See Fig.3for the zoom-in picture near the lower edge.The contributions of the two parity sectors are not separated.

Figure 3 :
Figure 3: Density of states from N Maj = 2N spin = 16 to 34 near the edge.The horizontal axis is E/|⟨E 0 ⟩ SYK |.The contributions of the two parity sectors are not separated.

Figure 5 :
Figure 5: The density of E ′ = E − E 0 (sample-by-sample subtraction) for SpinXY4 for N spin = 12, 13, 14, 15 and fit by ρ(E ′ ) = A sinh(B √ E) near the edge.The fit is sensitive to the fit range.

Figure 6 :
Figure 6: Distribution of the unfolded level spacing P (s i ) for i = 1, 2, 3, • • • .SpinXY4, N spin = 11 (top, left), SYK4, N Maj = 22 (top, right), SpinXY4, N spin = 15 (bottom, left), SYK4, N Maj = 30 (bottom, right).Only the parity-even sector was used.In SpinXY4, the eigenvalues in the parity-even and odd sectors are not correlated but the plots for the two sectors are indistinguishable.In SYK4, parity-even and parity-odd sectors have the same eigenvalues when N spin = N Maj /2 is odd.

Figure 8 :
Figure8: The spectral form factor for SpinXY4, N spin = 8, 9, . . ., 16 (left) and that for SYK4, N Maj = 2N spin = 16, 18, . . ., 32 (right).Only the parity-even sector is used.Note that SYK4 has a two-fold degeneracy in eigenvalues in each parity sector when N Maj ≡ 4 mod 8 and such a degeneracy shifts the height of the plateau by factor 2. The number of samples is 2 28−(N spin ) for both SpinXY4 and SYK4.

Figure 11 :
Figure 11: g(t) (left) and g c (t) (right) for SpinXY4 and SYK4 are compared.The results for β = 0, 1, 2 are plotted from bottom to top.[Top] N spin = 13.32768 samples are used for both models.[Bottom] N spin = 15.8192 samples are used for both models.The parity-even sector is used for both SpinXY4 and SYK4.

5 N 2 NFigure 14 :
Figure 14: |G xy (t)| plotted for the SpinXY4 model at β = 0 (left) and β = 2 (right) from top to bottom.The vertical axis is linear (logarithmic) in the upper (lower) plots.1024 samples are used, and the average over all operators and samples is taken before the absolute value is computed.

Figure 15 :
Figure 15: |G z (t)|, Re G z (t), and Im G z (t) for β = 0, 2 plotted for the SpinXY4 model.Note that Im G z (t) = 0 holds for β = 0, which is numerically confirmed.1024 samples are used, and the average over all operators and samples is taken before the absolute value is computed.
Nspin σ1,y ⊗ σ2,y ⊗ • • • ⊗ σNspin,y , which is real and symmetric if N spin is even.Therefore, if N spin Ĥ′ are real and symmetric, and hence, they are in the GOE universality class.When N spin is odd, there is no specific structure, and hence, we observe the GUE universality class.
28−N spin except for SpinXY4 with N spin = 17 and for SYK4 with N Maj = 34.