A new class of SYK-like models with maximal chaos

We investigate a model closely related to both the original Sachdev-Ye-Kitaev (SYK) model and the $\mathcal{N}=1$ supersymmetric SYK model. It consists of $N$ real Majorana fermions and $M$ auxiliary bosons with Yukawa interactions. We consider the large $N$ and $M$ limit and keep the ratio $M/N$ fixed. The model has two branches characterized by the conformal dimensions of fields, which we compute as a function of the ratio $M/N$. One of the branches contains the supersymmetric saddle for $M=N$. Furthermore, we determine the Lyapunov exponent of the model and find maximal chaos independent of $M/N$.

In particular the SYK model is a (nearly) conformal field theory (CFT) in the infrared, and is assumed to have a nearly anti de sitter (AdS) dual in this regime [10][11][12]. For these low energies the model can be described by a Schwarzian, which also appears on the bulk side in the AdS 2 dilaton gravity. The model has been intensively studied the past few years. There exists many generalizations including higher dimensions [13][14][15][16][17], flavours [18,19], tunable chaos [20] and supersymmetry [16,21].
In this paper we consider a particular model closely related to the N = 1 supersymmetric extension of SYK. Instead of having an equal number, N , of fermions and bosons we consider the case where we have M bosons and N fermions and study its behaviour as a function of the ratio M/N . In section 2 we will introduce the model and discuss in more detail the relation to the (supersymmetric) SYK model. Afterwards, in section 3, we consider the effective action. We derive the equations of motion and consider the solutions at strong coupling. We find two families of solutions that we label by their conformal dimensions at M = N (rational or irrational). Comparing the entropy of both solutions we determine that the rational solution is the dominant saddle for M = N . In section 4 we compute the Lyapunov exponent and find that is independent of M/N due to a subtle cancellation.

The Model
The model consists out of N Majorana fermions obeying {ψ i , ψ j } = δ ij and M (auxiliary) bosons. We will use indices a, b to denote the bosons and i, j, k for the fermions (no ambiguity will arise). The Lagrangian is given as follows: where ψ denote the Majorana fermions and φ the bosons. The coupling C aij is defined to be antisymmetric in the last two indices, which are contracted with the Majorana fermions. In [22] a similar term was studied as a perturbation upon "normal" SYK. The fermions are dimensionless, whereas the bosons φ and couplings C have dimension of E 1/2 . Notice that we have two parameters M and N . We are interested in taking the limits of both M and N going to infinity but keeping M/N fixed. In other words we have that M = αN for some fixed α. From now on we will always assume that two a indices are summed up to M whilst the other i, j, k, .. are summed up to N . Lastly, we let the coupling be disordered averaged by the following distribution: Here J has the dimension of energy and is larger than zero. We can now compute some basic one-loop diagrams for both the fermions and the bosons. We show the one-loop corrections to the two point functions in figure 1, which are proportional to some power of M/N (that can easily be checked). In fact one can check that any boson loop adds a factor of M N and each fermion loop

Relation to SYK
Let us first examine the relation to the original SYK model [1,23] with Hamiltonian To check the similarity we start by plugging in the algebraic equation of motion for φ a back into the Lagrangian. The equation of motion is found to be φ a = i 2 C a ij ψ i ψ j . After plugging it into eq. (2.1) we obtain the Hamiltonian: This is also the presentation that one can see in [22]. We can then use the antisymmetry in the last two indices of C aij and the commutation relations of the Majorana fermions to rewrite this to: where we defined the constant E 0 = − 1 16 C 2 aij (recall that a is summed to M and i, j up to N ). Comparing now to the standard SYK Hamiltonian, eq. (2.4), we find: (2.7) The notation indicates that the asymmetry on the right hand side is only in i, j, k and l, which in turn of course means that J ijkl is completely asymmetric. The above expression for the J coupling shows us that the model is essentially obtained by performing a Hubbard-Stratonovich (HS) transformation on SYK. Of course apart from this HS transformation we have also chosen a different distribution (see eq. (2.2)) compared to SYK. This means that J ijkl are no longer the independent Gaussian variables and this is the cause of the differences between the models.

Figure 2
This is the leading diagram that contributes to interactions between replicas.
In the figure the top three lines would be associated with a different replica index than the other three below. One can check that this figure is proportional to 1/N .

Relation to supersymmetric SYK
The N = 1 supersymmetric SYK model was introduced in [21], the Lagrangian density is given by: There are two important differences compared to the model described in eq. (2.1). The first important aspect is that there are N bosons, which is the same as the number of fermions (which has to be true for supersymmetry). Secondly the coupling C ijk in the supersymmetric case has to be completely antisymmetric. Note that the equal number of bosons and fermions is also necessary for the antisymmetry in the coupling.
In other words, starting from eq. (2.1) we can obtain the supersymmetric model by setting M = N and making the coupling completely antisymmetric. It is precisely when the coupling is completely antisymmetric (and hence M = N ) that the Lagrangian is invariant under supersymmetry transformations.

Effective action and saddles
To find the effective action we will follow the standard procedure of averaging over the disorder in C aij by using the replica trick (see appendices in [18,24]). As in the usual SYK case we will assume replica diagonal matrices. To justify this we have to compare log Z and log Z since assuming replica diagonal matrices corresponds to evaluating the latter instead of the former. The usual argument (see e.g. appendices in [24]) is to consider diagrams that are in log Z but not in log Z. The leading diagram belonging to the former but not the latter is shown in figure 2 and as can be verified it is suppressed by 1 N M . Thus in the large N limit these contributions will be subdominant. Using replica symmetry, the result of the disorder average becomes: For the last term we introduced brackets below the sum to indicate that j, k sum up to N whilst a sums up to M . We now introduce bilocal fields for both the fermions and bosons as follows: We insert them into the partition function by Lagrange multipliers. Afterwards we are only left with Gaussian integrals for both the fermions and bosons. Completing these leads to: Where on the left hand side we divided out a factor of N , but could just as well have taken out M (recall that M/N is fixed). Let us now vary with respect to G φ and G ψ to obtain the self energies: These equations can also be obtained using the melonic structure of the Feynman diagrams at large N and M , just as in ordinary SYK. The Schwinger-Dyson equations are obtained by varying with respect to the Σ (we assume time translation symmetry and go to Fourier space):

Two saddle points
In order to solve the above equations we have to assume the strong coupling limit βJ 1. This implies that in eq. (3.6) we can ignore the first terms on the right hand side. Hence we can write the equations as follows (we have Fourier transformed back to time): We then use the following (conformal) form for the two point functions: To obtain conditions on the conformal dimensions we plug these into the saddle point equations above, eq. (3.7) and eq. (3.8). Afterwards we Fourier transform using [21,23]: Some other useful relations for Γ functions are . (3.14) After plugging this al in we obtain the following relations: These relations (for M = N ) have also been derived in [21]. By comparing the frequency dependent parts we obtain the first condition on the conformal dimensions: As a side note, under this condition the saddle point equations have the conformal symmetry, very analogous to the original SYK model: Where f (τ ) a smooth function (in one dimension Conf(R) ∼ = Diff(R)). To obtain results for finite temperature we use this symmetry with f being the exponential map for example.
Coming back to eq. (3.15), we can obtain another constraint by taking the quotient, which yields the (transcendental) equation: This result, for M = N , is also obtained in [21], although it contains some typos. In [22] it is also shown for M = N , albeit in a different form. The second condition, eq. (3.18), can also be recast to an equation for ∆ φ using eq. (3.16): First we solve eq. (3.18) for M being equal to N . This case overlaps with supersymmetric SYK (as commented upon in the introduction) and we find the same solutions as in [21]. The first solution is given by:

20)
We label this solution as the "rational" solution. In the supersymmetric model this solution is the one that preserves supersymmetry. In that case, the supersymmetric Ward identity G φ = ∂ τ G ψ , together with eq. (3.9) implies ∆ φ = ∆ ψ + 1 2 [21], obviously obeyed by eq. (3.20).  There is another solution with positive conformal dimensions, it is however irrational:

21)
As one can easily check this does not satisfy ∆ φ = ∆ ψ + 1 2 and hence would break supersymmetry. A similar situation arises in [25] where there are also two solutions, one preserving and one breaking the supersymmetry.

Arbitrary M and N
Let us now vary the ratio M/N and find the conformal dimensions as a function of this ratio. We solve eq. (3.18) numerically and show the results in figure 3. There are two "families" of solutions, labelled by their behaviour at M = N . The rational solution was also found in [22].
When M/N becomes large the rational and irrational flow to the same point. This can be understood by considering the defining equations eq. (3.18) and eq. (3.19). When one takes the limit of M/N going to infinity there is only one solution left: Similarly, the behaviour for small M/N can be understood by taking the appropriate limits in the defining equations. This is equivalent to considering the limit N/M going to infinity in eq. (3.19). Consider the following two limits: Applying these in eq. (3.19) shows us that for small M/N we either need to consider the case where ∆ φ is very small or the case where it goes to one. The latter corresponds to the rational family. For the irrational ∆ φ 1 case we find from eq. (3.19) that it behaves exactly as: The corresponding dimensions for the fermion can of course be found by eq. (3.16). Lastly let us investigate the rational ∆ ψ for small M/N . Observing figure 3, we assume that ∆ ψ is small and consider eq. (3.15). We can then solve exactly as follows:

Dominant saddle
In this section we will determine what is the dominant saddle by comparing the entropies of both solutions. In particular we consider the case M/N = 1, since we know here exactly the behaviour of the rational solution as a function of q (see below). For the computation we will follow [21] and use the model for a q-interaction (meaning a vertex with one boson and q − 1 fermions, with q odd), see appendix A for an overview of the changes. The free energy becomes: Now we derive with respect to q (we continue the values of q to the reals) such that we don't have to evaluate the first terms. We take the fields to be on-shell such that we only need to explicitly take the partial derivative of the last term where the Gs are now the finite temperature versions, obtained by conformal symmetry mentioned before (eq. (3.17)): The integral can then be computed straightforwardly (using the periodicity in the τ variables): Where C is a constant independent of β. The constant term is a diverging quantity independent of q contributing to the ground state energy but will not contribute to the entropy, similar to the scenario in [21]. It is important to note that apart from the overall factor, the M/N dependence is also in A q−1 B (eq. (A.3)) and the conformal dimension ∆ ψ ( figure 3). For now we will consider the case M = N .
The entropies S R and S I are labelled by their rational of irrational origin, see eq. (3.20) and eq. (3.21) respectively, and given by: where C R and C I are integration constants. Furthermore we have called the q dependent parts R(q) and I(q) (these do not contain constants).

The rational branch
Let us first consider the rational branch. In this case we can always solve the exact dependence of the conformal dimension on q (see appendix A): The above expression means that the integrand in eq. (3.30), r(q), can be computed using eq. (3.29): and hence we can compute the entropy for the rational case: To fix the integration constant we will consider the limit q → ∞. For the case M/N = 1 we can follow exactly [21], section II.C. There the results in a large q expansion are obtained: Where v is an integration constant related to βJ [21]. The expansion in large q can still be made in a similar manner (although the functions g x (τ ) may now contain factors of M/N ). From these expressions above it becomes clear that the q → ∞ limit reduces to free fermions. It also allows us to fix the constant C R since: Where in the last line we used the above observation that it should reduce to a free fermion entropy.

The irrational branch
Unfortunately we can't solve analytically the q dependence of the irrational solution. We did manage to find a good fit by − 1 2q + 1 q−1 , which matches the numerical results very well for large q. In fact, we only find small deviations for low values of q. In figure 4 we plot the numerical results, the best-fit solution and the rational solution.
To conclude which of the entropies is bigger (i.e. which is the dominant saddle) we will investigate the integrands as a function of q, see figure 5. From this plot we can see that the irrational integrand is bigger than the rational one and as q increases their difference decreases.
We now investigate again the behaviour for large q. By using the approximate solution (see figure 4) we can obtain an expression for i(q) at large q. This is done by using the approximate solution in eq. (3.29): From this we can guess the leading order behaviour of I(q) in eq. (3.30): We can see that the irrational integrand (i(q)) is larger than the rational one (r(q)).
Where we also gave the rational behaviour for large q, which can be obtained from eq. (3.33). It can be seen that for these large q we have I(q) < R(q). Now to fix the integration constant for the irrational case (see eq. (3.30)) we can take the strict q → ∞ limit. In this limit the two solutions (irrational and rational) will coincide since their conformal dimensions will equal. Hence we also have that their entropies must coincide: (3.41) Using eq. (3.39) it becomes clear that the first term is zero and hence C I = C R = 1 2 log 2. We can then conclude that for large q: S I < S R . Further more since the slope of R(q) is always smaller than that of I(q) (since i(q) > r(q)) we find that this conclusion holds for any q.
To conclude this analysis, we found that the rational solution is the dominant saddle in this model at M = N . To extend the analysis for arbitrary M/N we would need to find a best-fit of the conformal dimensions as a function of q similar to above. Afterwards we can follow the same procedure to figure out the dominant one. Since this numerical analysis is quite tedious, we have left it for future work.

Chaos
In this section we will investigate the chaos or Lyapunov exponent of the model as a function of the ratio M/N . We will first review shortly the basics of such a computation and then move on to our model. The main tool for quantifying quantum chaos are so called Out of Time Order Correlators (OTOC) [5,[26][27][28][29]. For a more elaborate review of chaos and calculating these correlators see chapter 8 in [16], the first section of [8] and a discussion in [24].
From a quantum mechanical point of view we can take two arbitrary Hermitian operators V and W and consider the commutator [W (it), V (0)] (with real time t ∈ R). The argument of the operator is imaginary since we consider it to be Euclidean time, as will be the case for our operators later on. The commutator describes the influence of small changes of V on later measurements of W (or the other way around). One particular indicator of these effects of chaos, which we will also use, puts the operators on the thermal circle [8]: (4.1) Where the brackets denote the thermal trace, the precise factors of β will not be important for us. For late enough times t (to be precise, between the dissipation and scrambling time [8]), quantum chaos dictates that this correlator will grow exponentially. By considering all the terms that arise in the above correlator one can show [8] that the exponential growth of the correlator arises due to the exponential behaviour of the related correlator: These out of time order correlators F (t) are usually studied in the context of quantum chaos, and we will use these as well. Schematically the OTOC eq. (4.2) behaves as [8,16,28]: The exponent λ L is called the Lyapunov exponent and it quantifies the chaos of the system. In the coming section our goal is to extract this Lyapunov exponent from the OTOCs. In general one can follow two approaches. The most obvious one is to compute the full four point function and continue these Euclidean correlators to real time. An easier option, however, is to consider the so called retarded kernel and its eigenfunctions [1,23,28]. In the context of ladder diagrams, kernels are the operators that add one more ladder to the diagram. For the OTOCs it has to be the retarded kernel due to the complex time contours specified by OTOCs similar to eq. (4.1). For a review of this procedure including the complex time contours, ladder diagrams and the application to ordinary SYK see [16].
The key idea of this procedure is to consider an exponentially growing OTOC on which the kernel(s) are acting. Under the assumption of this exponential growth one can find that it is precisely the eigenfunctions of the kernel with eigenvalue one that govern the chaotic behaviour. More intuitively, the growth rate of OTOC is determined by the demand that adding another ladder should not change the total sum. In the rest of this section we explain this procedure in more detail.

Retarded kernels
Let us now turn to our model and consider the four point functions (or OTOCs) that we want to compute. In [30,31] the chaos is also calculated for similar circumstances and we will comment upon this method at the end. We will consider the four point functions ψ ψ ψ ψ , ψ ψ φ φ , φ φ ψ ψ and φ φ φ φ . This is because acting with kernels on these diagrams will result in mixing between them and hence we can not consider them separately. The explicit OTOCs we will consider are of the form: F ψψ (t 1 , t 2 ) = tr [y ψ(t 1 ) y ψ(0) y ψ(t 2 ) y ψ(0)] , (4.4) Where y is defined as y 4 = ρ(β). Diagrammatically these OTOCs are four point functions (ladder diagrams) with an arbitrary large amount of rungs. The other combinations of ψ and φ listed above have similar expressions and are denoted by F ψφ , F φψ and F φφ . The two subscripts of F ij denote the two incoming and two outgoing species, respectively.
Let us now consider all the (retarded) kernels necessary for our model, which we draw in figure 6. Note that there is no kernel K 22 since there is no such interaction in the Lagrangian. It then becomes clear that acting with the K 12 and K 21 kernels causes mixing between the four point functions. To get expressions for them we need first the necessary propagators in the diagrams. For the horizontal propagators we need the retarded ones (due to the complex time contours, (c) K 12 kernel Figure 6 Here we show the relevant kernels for the chaos computation. The subscripts of the kernels denote that they are elements of a matrix. The total matrix acts on a vector consisting of diagrams which starts with either two ψ or two φ lines.
see [16]): Recall that the arguments are imaginary since we consider complex Euclidean time. We can then use the finite temperature two point functions from eq. (3.28) to find: And similarly for φ: Lastly we need the ladder rung (lr) propagator, 1 which is obtained by simply continuing the Euclidean propagator τ → it + β 2 : Here x denotes ψ or φ and b x denotes A or B respectively. The form of this propagator is the same for fermions and scalars since we only need to consider τ > 0 here. We can then write down the expressions for the kernels. Note that each vertex gets a factor i from inserting it on a Lorentzian time fold in the contour and apart from this we also give K 11 and K 21 an additional minus sign due to the ordering of the contour (see also Here we show the matrix integral equation that the OTOCs obey. The black boxes indicate the arbitrary large amount of rungs in the ladder diagrams. The very left vector consists out of all the OTOCs, the first vector on the right hand side denotes the zeroth order contributions to these and the last term is the kernel acting upon the vector of the four point functions. The matrix product also includes a convolution between the kernels and the OTOCs. Notice that the 4 × 4 kernel matrix has a 2 × 2 block diagonal structure. [30,31]). The resulting form of the kernels is: 2 Where the times t ij = t i − t j are shown in figure 6.

Integral matrix equation
Now that we have obtained the retarded kernels we go back to our four out of time order correlators. All together they obey an integral matrix equation as shown in figure 7, this is a generalization of the one particle version seen for example in [16]. In the figure we have put all the OTOCs in a four component vector seen on the very left (and right) side. These are exactly the OTOCs we named F ψψ , F ψφ , F φψ and F φφ before. Our (drawing) conventions are such that for the very left vector the times t 1 and t 2 are on the top left and bottom left of each four point function, respectively.
The first vector on the right hand side denotes the free contributions to the four point functions. Clearly F ψφ and F φψ don't have these since there is no such free propagator. The matrix consists out of the retarded kernels discussed above and depicted in figure 6. Note that the matrix product in the last term also has an implicit convolution (which we will explicitly compute later on).
Quantum chaos implies that for late enough times these OTOCs will show exponentially growing behaviour, as discussed shortly in the introduction of this section. So let us make the following exponential growth ansatz: (4.10) Where i, j can denote ψ or φ and f ij denote functions of the time difference. Under the assumption of exponential growth the matrix equation figure 7 will simplify due to suppression of the zeroth order contributions. In fact, as one can easily check, the free diagrams will exponentially vanish compared to the exponential growth of the other terms. We are then left with the following equation: Where the F s now obey the ansatz eq. (4.10) and the matrix multiplication still involves the convolutions. However, we see that the problem can in fact be reduced to two identical problems because of the block diagonal structure. Hence we don't have to refer to the outgoing lines of the OTOCs (the second subscript of the F ij ) and consider simply: This leads to the following equations: (4.14) Where we have now explicitly written out the convolutions. The two equations are mixed and can be combined to give: F ψ (t 1 , t 2 ) = (K 11 * F ψ )(t 1 , t 2 ) + (K 12 * (K 21 * F ψ ))(t 1 , t 2 ) . (4.15)

Lyapunov exponents
To actually solve the integrals we need to find the functions f i (t 12 ) in eq. (4.10) such that eq. (4.12) is satisfied. We take the following form of the functions, similar to [16,23,30,31]: Where the C i denote non-zero real constants, and we have h as the free exponential growth parameter. The Lyapunov exponent can be found by λ L = − 2πh β . The crucial integral for the computations is as follows 3 Using the above identity we can then calculate the following integrals, reminiscent of eigenvalue equations: We use the following notation: Where the last line follows immediately from the relation between the conformal dimensions eq. (3.16). Then after some calculations we obtain the values of the k ij : , (4.21) (4.23) Note that here A and B are the coefficients of the (retarded) propagators, for which we only have an expression for A 2 B, see eq. (3.15). We can now use the above k ij along with eq. (4.18) in the integral equation eq. (4.15), to get: (4.24) Of course, one could also have used the eigenfunctions (eq. (4.18)) first in eq. (4.13) and afterwards solved the mixing. Either way the equation resulting from the chaos regime is: (4.25) We pick then some fixed M/N (which fixes ∆ ψ and A 2 B) and numerically solve this equation for the Lyapunov exponent λ L = − 2 π h β , which yields the solution h = −1. As it turns out, for h = −1 all the M/N dependence drops out and we find in fact maximal chaos for all values of M/N : (4.26) Motivated by these numerical results we analytically checked whether k 11 + k 21 k 12 = 1 for h = −1. To do so we use the identities from eq. (3.13) and also the following: . (4.27) Using these identities and the expressions for A 2 B, we obtain the following simplified expressions, valid at h = −1: Hence even though k 11 and k 12 k 21 individually depend on M/N (since ∆ ψ does), the combined result exactly cancels.
Lastly let us shortly mention another method of obtaining the chaos, outlined in [30,31]. In these articles the approach is to take the matrix of kernel eigenvalues, the k ij , and diagonalize it. Afterwards one of the eigenvalues is set to one. Let us consider this matrix: for which the resulting eigenvalues are k ± = 1 2 k 11 ± k 2 11 + 4 k 21 k 12 . The growing behaviour is found when k + = 1, which amounts to k 11 + k 21 k 12 = 1, consistent with our method.

Discussion
In this article we have investigated new SYK-like models with M bosons and N fermions. The parameter M/N determines the behaviour of the model and for M = N our model is related to the supersymmetric SYK model. We have found that there are two families of solutions in the model, distinguished by their conformal dimensions which we plotted as a function of M/N in Figure 3. For M = N the rational solution coincides with the supersymmetric solution found in [21]. We have shown that this branch is the dominant saddle for M = N . It would be interesting to see if this remains true for arbitrary M/N , it might be that at some point the two branches switch from being (sub)dominant.
Apart from this we investigated the Lyapunov exponent λ L = − 2πh β , which shows exactly maximal chaos independent of M/N . This is due to some non trivial cancellations in the M/N dependences at h = −1. Due to the maximal chaos the model has a holographic interpretation as a black hole. It would be interesting to understand the role of M/N in this holographic description. Concretely it would be interesting to find the Schwarzian for this model, in particular the coefficient in front of the Schwarzian action, related to the heat capacity, and its dependence on M/N . We leave this for future research.

A The model for a q-point interaction
In this appendix we will shortly show how the model and some results change when we consider an interaction vertex of degree q, so an interaction with one boson and q − 1 fermions. The integer q is supposed to be odd, but later we will continue it to arbitrary real values. The model we consider in the main text has q = 3. When we apply this, the coupling i 2 C aij φ a ψ i ψ j goes to i (q−1)! C ai 1 i 2 ...i q−1 φ a ψ i 1 ...ψ i q−1 . Integrating out the bosons would lead to a Hamiltonian with a vertex containing 2q − 2 fermions. We further assume that q N .