Chaotic Correlation Functions with Complex Fermions

We study correlation functions in the complex fermion SYK model. We focus, specifically, on the h = 2 mode which explicitly breaks conformal invariance and exhibits the chaotic behaviour. We explicitly compute fermion six-point function and extract the corresponding six-point OTOC which exhibits an exponential growth with maximal chaos. Following the program of Gross-Rosenhaus, this correlator contains information of the bulk cubic coupling, at the conformal point as well as perturbatively away from it. Unlike the conformal modes with high values of h, the h = 2 mode has contact interaction dominating over the planar in the large q limit.


Introduction
A generic dynamical system is inherently chaotic [1]. For classical systems, chaos can be easily characterized by the sensitivity of trajectories with respect to initial conditions. For quantum systems, lacking in the concept of trajectories, the notion of chaos is more subtle. Often, quantum chaos can be characterized in terms of properties of the spectrum of the Hamiltonian. In the semi-classical approach, there is a relatively simple definition of chaotic behaviour, directly adopted from the sensitivity of classical trajectories with respect to initial conditions.
For classical dynamical systems, characterized by phase space coordinates {q(t), p(t)}, where q(t) and p(t) are generalized positions and generalized momenta. A particular trajectory is represented by q(t). High sensitivity of the late time trajectory with respect to the initial condition can be quantified as: where λ L is the so-called Lyapunov exponent and the right-most expression above is the Poisson bracket [1]. By virtue of the correspondence principle, we obtain a quantum mechanical characterization, by replacing the Poisson bracket with a commutator: {q(t), p(t)} → −i [q(t), p(t)] [2]. Instead of computing the commutator, one calculates the squared commutator, so that there is no spurious cancellation due to destructive phases. This argument, however, is limited and does not necessarily imply that allowing for such phases will always cancel the chaotic growth. In this article, we will calculate the cubic power of the commutator, which will explicitly display the exponential growth behaviour.
We define a generic function for the diagnostic of chaos: where n ∈ Z + , and V and W are two self-adjoint operators and the expectation value is defined with respect to a particular state of the system. Note that, in defining the chaos diagnostic in (1.2), we have recast the chaotic property as a feature of n-point correlation function of the system. In this article, we will explicitly discuss the case for n = 3 in a thermal state.
Before doing so, let us briefly look at the n = 2 case. Written explicitly, the commutator contains various four-point functions with no particular time-ordering, since t 1 and t 2 are defined without any ordering. For a thermal state expectation value, using the KMS conditions 1 , it is further possible to rearrange the various four-point functions in terms of two pieces: one timeordered four-point function and another out-of-time-ordered correlator (OTOC). These are given by V (0)V (0)W (t)W (t) and V (0)W (t)V (0)W (t) , respectively, choosing t 1 = 0 and t 2 = t. The time-ordered correlator does not display the exponential growth, it is contained in the four-point OTOC.
For n = 3, upon using the KMS condition, the chaos diagnostic in (1.2) has one time-ordered and two OTOC pieces. These are simply, V (0)V (0)V (0)W (t)W (t)W (t) (time-ordered) and V (0)W (t)V (0)W (t)V (0)W (t) , V (0)W (t)V (0)V (0)W (t)W (t) , W (t)V (0)V (0)W (t)V (0)W (t) , etc, which are OTOC. While a complete understanding of the behaviour of (1.2) for arbitrary n is desirable, we will explore an exact calculation for n = 3 in this article, with a particularly simple model. 2 Note that, the correlator of the form V W V W V W requires three time-folds in the Schwinger-Keldysh framework, while an OTOC of the form V W V V W W require only two time-folds. These OTOCs, generically, would not carry the same physical information. For the present model, however, we expect the same four-point OTOC physics in the second case. 3 The model we consider is a simple generalization of the so-called Sachdev-Ye-Kitaev (SYK) system [3][4][5][6], in which one considers fermionic degrees of freedom with an all-to-all interaction. The interaction coupling is drawn from a random Gaussian distribution with a zero mean value and a given width. In the large N limit, in which the number of fermionic degrees of freedom becomes infinite, the system becomes analytically tractable in the sense that the corresponding Schwinger-Dyson equations can be explicitly determined. The solution of this equation readily determines the two-point function, as a function of the coupling strength, in general. In particular, in the low energy limit, this Schwinger-Dyson equation is analytically solvable and yields a two-point function with a manifest SL(2, R) symmetry. In the infra-red (IR), this is described by a conformal field theory (CFT), and the two-point function breaks the conformal group into the SL(2, R) subgroup. In the large N limit, further, the four-point correlator can be explicitly calculated, which yields the corresponding Lyapunov exponent: λ L = 2πT , where T is the temperature of the thermal state. Here, we are working in natural units. This Lyapunov exponent saturates the so-called chaos bound [7]. Intriguingly, the chaos bound saturation also occurs for black holes, in which the local boost factor at the event horizon determines the corresponding Lyapunov exponent as well as the corresponding Hawking temperature. Only extremal black holes have an SL(2, R) global symmetry, due to the existence of an AdS 2 sector near the horizon. Correspondingly, the low energy conformal system coming from the SYK model can be shown to capture the essential physics of the AdS 2 [8].
The low energy effective action for the SYK model is simply given by a Schwarzian effective action, which can also be shown to arise from the two-dimensional Jackiw-Teitelboim theory in [8,9](In fact, exact derivation of the Schwarzian effective action has been done using the fluid/gravity correspondence [6]). However, in this context, the non-trivial statements of holography necessitates keeping a leading order correction away from the purely AdS 2 throat, as well as from the purely CFT 1 in the IR, hence it goes by the acronym of NAdS/NCFT. From the geometric perspective, AdS 2 appears in the following two cases: (i) in the extremal limit of 2 Note that, our explicit expressions of the OTOCs are written in an abuse of notation: KMS conditions will render some of the time-arguments to pick up imaginary parts, proportional to the inverse temperature β. We suppress these explicit factors here, for simplicity, but take those into account for the explicit calculation of the six-point correlator later. 3 We thank the Referee for raising this issue.
a black hole in asymptotically flat background, (ii) in the deep IR of an asymptotically AdS d+1background. Often, in the second case, the deep IR results from an RG-flow connecting a UV CFT d to an IR CFT 1 , as a result of a relevant density perturbation in the UV CFT. Holographically, such operators correspond to turning on a bulk U(1)-flux, in the simplest case. A standard example of this is the AdS-Reissner-Nordstrom black hole: It asymptotes to an AdS geometry and the extremal limit consists of an AdS 2 in the IR, which is supported by the flux. In the AdS 2 throat, the flux is a simple scalar field, and the boundary theory has a natural notion of a conserved charge and therefore a non-vanishing chemical potential.
The SYK model which is defined with N Majorana fermions, ψ i (i = 1, ..., N ) in (0 + 1) dimensions, cannot have a charge density and correspondingly a non-vanishing chemical potential. If we, instead, consider N complex fermions, ψ i and ψ † i , it is natural to introduce a ψ † i ψ i term in the Lagrangian, with chemical potential as the coupling. Similar to the SYK model, we consider a q-body interaction term, with q/2 number of ψ and q/2 number of ψ † , whose couplings are chosen from a Gaussian distribution, with a standard deviation denoted by J. In the limit q 1, this system becomes exactly solvable in an (1/q)-expansion [4], which is what we will use in evaluating the explicit correlation functions. Motivated from the previous paragraph, it is therefore natural to consider a particular generalization of the SYK model with a U(1) global symmetry.
The standard SYK theory, with Majorana fermions, has a particular operator at the conformal fixed point, whose four point OTOC displays chaotic behavior with Lyapunov exponent λ L = 2π β [3][4][5][6]. The generalization of this model to complex fermions was done in [38][39][40]. In [38], the low energy effective action of an SYK-model with complex fermions was discussed. It was shown that the presence of a non-vanishing chemical potential does not break the conformal symmetry in the deep IR, as long as one amends the conformal transformations with a gauge transformation. Thus, the resulting low energy effective action is simply a Schwarzian action along with a free bosonic theory with a standard kinetic term [14,17]. Therefore, from a strict IR perspective, the maximal chaos holds for any non-vanishing value of the chemical potential.
The non-triviality comes from the order of limits. The exponential growth of the four-point OTOC holds for a larger regime compared to the long-time (and therefore, deep IR) limit. For sufficiently large time, one recovers the maximal chaos. However, there exists an intermediate regime in which the exponential growth takes place with a different Lyapunov exponent. This is physically equivalent to staying in a medium energy scale and finding a chaotic behaviour of the correlation function at this energy scale. This associates naturally an RG-flow of the Lyapunov exponent itself.
In the standard SYK model, in the large q limit, the relevant scale in the system is provided by an effective coupling: qJ 2 2 q−1 , which has mass dimension one. The IR CFT resides in the J → ∞ limit, but the exponential growth of OTOC and subsequently the Lyapunov exponent can be obtained as a perturbation series in 1/J . This naturally gives an RG-flow of the Lyapunov exponent [4]. In [40] we studied the SYK model with complex fermions in the large q limit in the presence of a chemical potential µ. Here, in the UV Hamiltonian, we have two natural parameters: βµ and βJ and the effective coupling in the IR is given by, The strict IR is located at J eff → ∞ limit, and one can calculate systematically the RG-flow of the Lyapunov exponent in a perturbation series in 1/J eff . This RG-flow shows sensitive behaviour for the Lyapunov exponent as the UV parameter βµ is dialled up [40].
In keeping with the theme, in this article, we further compute higher point OTOC for complex fermion SYK-model, with a non-vanishing chemical potential. Our analyses follow closely the analyses in [41], in the large q limit. However, our analyses are performed in the complementary regime in that we completely focus on the operators that display chaotic nature and away from the conformal limit. In spirit of the NAdS/NCFT picture, this is rather natural regime to consider; in the context of chaotic properties of many body systems, this is an example of a tractable and explicit higher point OTOC which displays the expected exponential growth.
It is worth noting that, even though our calculations are performed with a non-vanishing chemical potential, all analyses are still valid near the conformal point. Therefore, it is not surprising that we find a maximal Lyapunov exponent. On the other hand, our numerical methods are expected to extract a Lyapunov behaviour only at early times. We demonstrate, with explicit examples, that such an early time calculation still captures a chaotic growth, with a (Schwinger-Keldysh) two-folded 6-point correlator. Furthermore, it captures the maximal Lyapunov with a (Schwinger-Keldysh) three-folded 6-point correlator. At present, we do not have a deeper understanding of an underlying, universal physics that may be responsible for this. However, we do certainly take note that, if generic, this provides us with a tremendous computational convenience in extracting similar chaotic growth, both qualitatively and quantitatively, in various other systems as well.
In this paper, after computing the fermion six point function with a non-vanishing chemical potential, we take the triple short time limit to estimate the the bulk three point correlator, away from the conformal limit. In this regard, we compute bulk three point function(triple short time limit of the fermion six point correlators, neglecting the Schwarzian mode) of the modes satisfying conformal invariance as well as the Schwarzian mode, using the techniques employed by Gross and Rosenhaus [41].
This paper is organized as follows. In Section 2, we briefly review the SYK model with complex fermions. In Section 3, we compute the six point fermion correlator in the triple short time limit. We then interpret it in terms of the bulk three point correlator in the IR limit of the conformal modes and check that we do indeed find them to be of the form of conformal three-point function, in the triple short time limit. We apply this technique in Section 4 to compute the six point function away from the conformal limit and we look at the enhanced contribution due to the non-conformal mode. In Section 5 we take the triple short time limit to determine the three point correlation function of fermion bilinears away from the conformal limit. We carry out this computation in the presence of a chemical potential µ. We also compute the relevant Eucledian correlators which on analytic continuation to real time, is supposed to give us the two fold OTOC as well as the three fold OTOC. From this data we extract the Lyapunov exponent. We conclude with the discussion of our results, and possible future directions.

SYK model with complex fermions
The SYK model with complex fermion in 0 + 1 dimensions is defined by the Hamiltonian with all to all random interaction between q fermions, An exhaustive study of this model is done in [39], we will mention some of the essential features that will be necessary for our analysis. In addition to the higher dimensional operators of t ψ i which behave in a manner similar to those found in the SYK model with Majorana fermions; we also have the operators of the formÕ n = 1 The lowest lying mode of these operators give the Schwarzian mode and the U (1) charge respectively. In absence of a mass like term in the action the two point function of the particle and anti-particle are the same in the free case as well as the low energy limit of the interacting theory.
where, G c (τ ) is the propagator in the conformal limit. In the low energy, i.e., IR limit it is possible to obtain the four point function of the fermions using the expansion in the eigenbasis of the quadratic Casimir operator. We will skip the details and only state the results here. Since we have complex fermions, i.e., ψ i = ξ i + iη i , in case of the correlation functions we have contributions of two different kind, While the first piece, namely, the real part is anti-symmetric under the exchange of t 1 and t 2 , the second piece is symmetric.
In case of the four point function if we consider the time reversal invariant contribution this leads to two different contributions namely F A (τ 1 , τ 2 , τ 3 , τ 4 ) and F S (τ 1 , τ 2 , τ 3 , τ 4 ) which are respectively anti-symmetric and symmetric under t 1 ↔ t 2 and t 3 ↔ t 4 . The first term, i.e., F A (τ 1 , τ 2 , τ 3 , τ 4 ) is identical to the SYK with Majorana fermions but the second term is new and occurs in the complex fermion model. From [39] we have, is the conformal cross ratio and Ψ A and Ψ S are linear combinations of the eigen-functions of the quadratic Casimir. Ψ A is antisymmetric, while Ψ S is symmetric under the transformation, which effectively exchanges the first two or last two arguments of four point function. Finally k A and k S are eigenvalues of the retarded kernels (for antisymmetric and symmetric) which commute with the Casimir.

Large q limit
We now augment this q-point interaction with a quadratic coupling term by introducing a chemical potential µ which couples to the conserved charge i ψ † i ψ i . The fermion propagator in the Fourier space derived from the quadratic part of the action is, Once we take the interaction terms in the Hamiltonian into account it gives the dressed propagator. In the large N limit, the contribution of melonic diagrams dominates. If we also take large q limit (q < N ) then the loop corrected propagator is amenable to analytic computations, where, G 0 (µ, τ ) is the free propagator in real space. To compute the two-point function in the interacting theory one sets up the Schwinger-Dyson equation and seeks a solution in the large q limit. This Schwinger-Dyson equation at finite inverse temperature β can be cast in the form of a differential equation, The solution to this equation in given by [40], where J eff is defined in (1.4).

Correlation Functions
Let us begin with the short time, i.e., τ 1 − τ 2 = τ 12 → 0 limit of the four point function both for the symmetric and anti-symmetric case, When we calculate the the six point function of the complex fermions we go to different short time limits, where the correlation functions take some effective form. In the triple short time limit we calculate it as an effective three point function of the fermion bi-linear operators. This way one can compute the correlation function near points where different arguments approach each other yielding poles and by the property of being analytic everywhere else we get the full contribution.
In the remaining part of this article we calculate the O(1/N 2 ) coefficient of the six point function with respect to the 1/N expansion. To this order there are contributions from the contact diagrams as well as from the planar diagrams (see Fig.1). We will now write down the corresponding expressions: where the contributions of S 1 (contact) and S 2 (planar) are exactly same as those in [41], i.e., they agree with corresponding results for the Majorana fermions. In case of the SYK model with complex fermions, if we demand time reversal invariance, (since the Hamiltonian is itself time reversal invariant) we have only two other contributions. Now, 3) is the contact diagram contribution. Here, we have written only one particular assignment of the arguments; there are other possible assignments whose contributions account for the factor of 1/15 on the left hand side. There are total 15 possible independent configurations 4 .
We will use same symbol h to denote the conformal weight of the bi-linear operators both for F A and F S , although the values are different for the two: for The planar contribution is given by, where, Using the Selberg integrals in its special and general forms, one obtains: Using the short time expansion of four point amplitudes, we get: where explicit expressions of the constants c, ξ,c,ξ are given in appendix A. The integrals I (1) and I (2) are given by, (3.9) 4 Instead of taking τ 1 , τ 2 for the first four point function argument we could take any other two, for example τ 1 , τ 3 . Also interchanging the vertices among themselves This way the total number of possibilities is given by 6 The integral (3.8) can be simplified by the change of variables, τ a → τ 1 − (1/τ a ), and τ b → . The simplification is done by first decomposing the integral into sums of integrals. Namely the integration from −∞ to ∞ will be written as a sum of two, an integral from −∞ to 0 and an integral from 0 to ∞. We implement the change of variables on each fragment separately, simplify each of them before recombining them back. At the end of this exercise we get, (3.10) These change of variables are followed up by another pair of change of variables which are carried out in a sequential manner. We will first implement τ a → τ a − (1/τ 31 ), and τ b → τ b − (1/τ 31 ) and then we will rescale the integration variables τ a → (τ 53 τ a )/(τ 31 τ 51 ) and where, α = −h n + 1, β = −h k + 1, and γ = hn+hm+h k 2 − 1. To summarize, the aim of the above exercise was to obtain a conformal three point correlation function for fermion bi-linear operators denoted above by I (1) nmk . As in [41], we divide the Selberg integral,S f ull 2,2 , into different parts. This is achieved by decomposing the integral into three pieces [−∞, 0], [0, 1] and [1, ∞] for each integration variable. This results in six Selberg integrals with appropriately modified arguments. Carefully keeping track of the signs gives, (3.12) The generalized Selberg integrals and some important results which are used above are given in [41], but for completeness we give the relevant definitions here, (3.13) In a similar fashion one can manipulate I (2) nmk to bring it in a form of the conformal three point function. This computation, however, is considerably more involved, we instead do the analysis in the large q limit. The I (2) nmk in our case differs from that obtained in [41] by only the sgn functions while the rest of the integrand has exactly the same form. So for us also at large q, I (2) nmk takes the form, (3.14) In our case of courses (2) nmk is different from s nmk obtained by Gross and Rosenhaus [41].

Away from the Conformal Limit
In this section we carry out the calculation of correlation functions away from the conformal IR fixed point. In our earlier work [40], we studied the effect of introducing a chemical potential µ, in the SYK-model with complex fermions. We found that a non-zero µ takes us away from the conformal limit since it explicitly introduces a scale in the problem. The effect of introduction of this scale parameter is reflected in the chaotic behavior of the model, namely, it brings down the value of the Lyapunov exponent. We computed the required quantities and studied the maximally chaotic mode(in the large q limit where things could be handled analytically).
We write below the relevant expressions in the large q limit. The two point function(to the leading order in large N ) in the region τ > 0 is given by, where, The above relation can be written in a compact manner by writing, We now aim at calculating the enhanced contribution to the four point function slightly away from the conformal limit with the chemical potential µ. Note that since we want to be slightly away from the IR, we will keep µβ to be small and expand all functions in this variable. Then it can be interpreted that we move slightly away from the IR by turning on a small chemical potential.
The solution to this equation with appropriate boundary condition is well known. In fact this is the same equation as obtained in [4]. The solution is given by, (withñ = n/ν) The quantization condition on h is obtained by demanding that the wave function vanishes at x = 0, i.e.,x = π(1 − ν). As we approach the conformal limit ν → 1 this solution actually diverges for generic values of h near 2 (we are interested in the h = 2 eigenfunctions). But we want values of h such that the solutions are finite or vanishing, so the first or second argument of the hypergeometric function has to be a negative integer. This gives the quantization of h near 2 to be, This gives the shift in the eigenvalue k = 2 h(h−1) to be, This result is identical to the shift obtained in [4], only difference being that ν now depends on the effective coupling βJ eff instead of βJ .

The enhanced four point contribution
Let us now look at the four point function and use the above result to figure out the enhanced contribution for the Schwarzian mode slightly away from the conformal limit. We begin with the expansion of the four point function in the basis of eigenfunctions of the Kernel (using the variable θ = 2πτ β on the thermal circle and the period becomes 2π), To find the enhanced contribution of the Schwarzian or h = 2 mode we use the eigenfunction of the Casimir for h = 2 and the shifted eigenvalue in the denominator. In the numerator we just use the eigenvalue with h = 2 in the IR. This is done to ensure that we are only slightly away from the conformal limit driven by introducing a small chemical potential. Here we will use all results for the large q limit, where, and for large βJ eff we have used 1 − ν ∼ 2 βJ eff ≈ 2 βJ + q 2 − 1 (µβ) 2 4βJ . We will now carry out the sum over n. The final expression after all simplifications can be written in a compact form [4]. To proceed with the six point function using these results we use numerical computation, in carrying out the integration. Subsequently, we extract the chaotic behavior of the OTO six point correlator, even though we do not have an analytic result. This will be discussed in a later section.

The Contact and the Planar diagrams
There are two types of topologies of Feynman diagrams contributing to six point function. The contact diagrams and the planar diagrams (see Fig. 1). We will now claim that among these diagrams, which contribute to the six point function at the leading order in ∼ 1/N , the contribution of the contact diagrams dominates over that of the planar ones by an order q 4 , for the enhanced non-conformal mode contribution to the four point function. So at large q, the contact diagrams dominate over the planar ones, and hence we will look at only the former. But before that we will demonstrate this dominance of the contact diagrams below by studying the scaling behaviour in the large q limit of the contact and planar contributions.
The contact contribution is given by, (4.16) Whereas the planar contribution is, where, is the amputated four point function. We can use the SD equations to write, Since we are working at finite temperature, we have to do the Matsubara sum. However, notice that the iω + µ term has no poles, so when we evaluate the sum using the contour integration prescription, contribution of this part vanishes and we are left with, Now we can convert the τ integrals to θ integrals via appropriate scaling and we get the contribution from the contact diagram to be, where, In terms of the θ variable we have, (4.22) and in the large q limit, for large but finite βJ eff , (4.23) Here we have put ν = 1 inside the sine function which is consistent to the leading order with ν → 1 as βJ ef f → ∞. As a consequence the (βJ eff ) 2 coefficient of the contact diagram as well as the amputated four point function cancels out due to the (βJ eff ) 2 appearing in the denominator of (4.23).
The planar six-point diagram, on the other hand, is given by the product of three amputated four point functions hence, when we take the product and convert the τ integrals to θ integrals in (4.17), we get, Taking the ratio of S c with S p we see that, So in the large q limit as one can easily see that the contact diagram is far more dominant compared to the planar ones and hence it is justified to consider the contribution of the contact diagrams only.

The six point function
Although one can get an analytic answer for the enhanced four point contribution slightly away from the conformal limit, calculation of the full six point function becomes somewhat messy to carry out analytically. We therefore compute the six point function using numerical methods.
Let us first summarize the results, we will then we state all the relevant values used in carrying out these computations.
• We first compute the six point contribution with three h = 2 modes, for a fixed µβ, keeping all the time arguments separate and then we take the short time limit θ 2 → θ 1 , θ 4 → θ 3 , etc. We then change the value of µβ and compute the correlator again. We see that the six point function decreases in this limit as µβ is increased, while keeping βJ fixed at some large value.
• We can instead first take the triple short time limit and then carry out the integrals numerically for all the possible nonconformal contributions, i.e., We find that among the three terms in (4.26), the first contribution seems to vanish to all orders in µβ expansion, on the other hand the remaining two terms, although are small, have contribution at the same order, namely O ((µβ) 2 ). These results will get corrected as we go to higher orders.
• To benchmark the code we compute λ (1) 11k (as was done in [41] for all three conformal modes) 5 for the contact diagrams and plot it against k, where for large k, h k = 2k + 1 + 2∆ + O(1/k). We find the similar fall off behavior at large k.
Let us now look at some details of the analysis. One of the things that we have to keep in mind is that we are slightly away from the conformal limit because we have turned on a small µβ. We need to be careful while working with the conformal modes. Due to explicit scale in the theory, the modes may not be conformal anymore. In other words, the normalized four point contributions of these modes may not be a function of only the cross ratio χ anymore. However, for small βµ, the conformal perturbation theory makes sense and within this limit using the conformal basis is justified.
If we recall the eigenvectors of the Kernel then we see that, and to obtain the conformal modes one has to go to the IR, do the sum over n in the four point function to obtain the sum over integer values of h as well as the integral over the principle continuous series. One then deforms the contour to pick up the poles at k c = 1, eigenvalue of the kernel in the conformal limit. In the IR limit the exponential µβ factor becomes equal to 1, but since it has no n dependence it plays no role when we carry out the sum over n. Therefore, slightly away from the conformal limit we will have (small µβ), The above expression breaks conformal invariance and this is the four point function we will be working with away from the conformal limit.
Since βJ eff appears as an overall factor we strip off this factor and look at the integrals only. For small µβ this factor is large but finite. We will keep µβ fixed at µβ ∼ 7.4 × 10 −4 .

The Short time and OTO behavior of the Six point function
In this section we will study short time limit of the six point function as well as its out-oftime ordered behaviour. In both cases the expressions are quite involved and we had to take recourse to numerical methods. We will first look at the short time limit. Let us first look at the short time limit of the six point function. Since the terms we will be looking at are a part of the non-conformal piece, we will first compute it by keeping all six times different. We will then take the limit in which one of the temporal variable is approaching another, for example θ 2 → θ 1 , while holding rest of the time variables fixed. In this limit we will study the behaviour of the six point function. We do not have analytic control over this computation and hence we take a recourse to the numerical methods. Figure  2 contains the plot of numerical evaluation of the behavior of the six point, involving only the h = 2 modes, as a function of θ 1 in the short time limit, θ 2 → θ 1 while keeping θ 3 = θ 4 , and θ 5 = θ 6 fixed. In the triple short time limit, it is easy to see that the six point function vanishes.
While the contribution to the non-conformal piece coming from the product of h = 2 modes is shown in Fig. 2, computing the contribution for one or two h = 2 modes requires some care because in this case the planar diagrams may have enhanced contribution but in the large q limit, it will still be subdominant. In a similar way it can be checked that the planar contribution is finite by analysing behavior of F amp (θ 1 , · · · , θ 4 ) as a function of its arguments. However, we do not explicitly compute this because we do not need it in our analysis.

Maximal Lyapunov from six-point OTOC
Before discussing the details, let us emphasize that even though we have analytic control over the ingredients, thanks to the large q-expansion, the calculation of the Euclidean correlator involves performing multi-variable integrations over the Euclidean time. The resulting integral, for the 6-point correlator, is a function of six Euclidean times. Ideally, one needs to first obtain this function of six variables, and subsequently carry out the analytic continuation corresponding to the Lorentzian correlator that one wants to compute. However, this is a challenging task in general, and it is even more so when one lacks an analytic expression for the Euclidean 6-point correlator. In view of this, we will discuss below a simple-minded numerical approach that appears to capture the correct Lyapunov behaviour.
Before considering the 6-point function, let us concentrate on the 4-point function evaluated in [4]. It is known to produce to well-established maximal value λ L = 1, in units of T = 1/(2π). The assignment of the time arguments in the alternating configuration, in [4], was set to be: This implies that both θ 1 and θ 2 are analytically continued to real time and from their assignment on the thermal circle we find figure 3 as the corresponding Schwinger-Keldysh contours. Clearly, in this case, one computes a two Lorentzian time-folded correlator which corresponds to the true 4-point OTOC. The ordering of the imaginary time co-ordinates do not matter when analytically continuing to real time. There the ordering is fixed by epsilon prescription. The fact that the above chosen configuration yields the expected value of the maximal Lyapunov exponent, suggests that the information of the epsilon prescription can be translated to the configuration we choose. Furthermore, note that, even though the general 4-point OTOC is a function of four real-time variables, the choice in (5.1) simplifies this dependence to only one real-time variable, denoted by t. A similar statement holds for the Euclidean correlator.
Given the above choice of time-assignments, we perform the following tasks: (i) We numerically generate Euclidean data points, this yields a function F 4 (θ), where F 4 is the Euclidean four-point correlator and θ is the only Euclidean time variable, with our proposed assignment in (5.1). (ii) We numerically fit the data with a guess-function. By trial and error, and motivated by simplicity, we have used a cosine function as the fitting function, in which the frequency and the amplitude of the fitting function are found out as a result of the numerical fitting. (iii) We read off the frequency to be the corresponding Lyapunov exponent, since where () 1,2 are overall constants. In the large time limit, such a term is dominated by the e λ L t , and therefore the Lyapunov is obtained from the frequency of the Euclidean periodic function.
We have explicitly checked that, the above numerical implementation does reproduce λ L = 1, within acceptable numerical accuracy, for the OTOC considered in [4]. This is not a watertight argument that such an analysis will work for arbitrary cases, but, in absence of a better method, worth exploring. Furthermore, note that one does not need the complete information about the frequency spectrum of the Euclidean correlator, but only the dominant one, to extract the maximum Lyapunov exponent. This may provide a justification of why this method should work, however, we do not have any further reasons to offer at this point.
We will subsequently use a similar numerical method to compute the 6-point OTOCs, both with two time-folds and three time-folds. Finally, we will fit the data and extract the maximum Lyapunov exponent.

2-fold Schwinger-Keldysh contour
First we compute the quantity: which corresponds to the OTOC with two time-folds and τ denotes the Euclidean time. Thus, this is not a truly 6-point OTOC, since the out-of-time physics is captured by a 4-point OTOC that lies inside the 6-point correlator. Motivated by our discussion above, in this case, we assign: and we choose τ = −0.5, −0.4, ...., 0.5. Thus we have, Note that, within the regime of the τ -variable: τ ∈ [−0.5, +0.5], we ensure the ordering: θ 1 < θ 3 < θ 5 and therefore θ 2 < θ 4 < θ 6 , by choosing the configuration in (5.3). This ordering ensures that there are no operator crossing and the corresponding correlator remains a twofolded OTO correlator, for the regime of τ considered above. We note that this is only a representative configuration, however, the qualitative physics does not appear to significantly change.
We set the following values for different quantities: The 2-fold contour, corresponding to the assignment in (5.3), is shown in figure 4. Now we compute the numerical data, as a function of τ , the result of which is shown in figure 5. To this data, we fit a function of the form: Clearly, λ L is not far from the maximal value, but it is not sufficiently close.

3-fold Schwinger Keldysh contour
Now let us consider the quantity: which is a true 6-point OTOC. The corresponding assignment of the Euclidean coordinates are given by This implies: . This configuration corresponds to the following 3-fold contour: As before, it is straightforward to check that, for the range of τ ∈ [−0.5, +0.5], the configuration in (5.6) respects the ordering: θ 1 < θ 3 < θ 5 and therefore θ 2 < θ 4 < θ 6 . Thus, the 3-folded OTO correlator remains a 3-folded one, for the entire variation of τ ∈ [−0.5, +0.5].
The numerical result is plotted in figure 7, along with the fitting curve, of the form: From the fit we obtain the following values of the coefficients and frequency: Once again, the value of the corresponding Lyapunov is close to the maximal one.
Before concluding this section, let us note that our calculations are done near the conformal point and therefore it is not surprising that the corresponding Lyapunov exponent is numerically close to its maximal value. However, a few comments are in order.  Figure 7. Six point OTOC with 3-fold Schwinger-Keldysh contour.
Step-size= 0.05, λ L = 1.01853 Firstly, note that all our calculations are performed in the Euclidean description. By construction, the Euclidean time τ ≤ β, where β is the inverse temperature. At best, we can ascertain τ ∼ O(β). Naively, within this regime, the analytic continuation to the Lorentzian time, via τ → it, should capture physics at time-scales t ∼ 1/T ∼ t d , where t d is the dissipation time-scale. Therefore, our method is expected to capture the early-time physics that is still close to the physics at the dissipation time-scale.
It is therefore quite interesting that, at a pragmatic level, both the two-fold and the threefold OTO correlators seem to capture an exponential growth even at this early time-scale. Moreover, the three-folded OTO seems to suggest a maximal Lyapunov exponent already at these time-scales, while the 2-fold OTOC (which is a four-point OTOC, fused with a two-point correlator; and not a truly 6-point OTOC) suggests an exponential growth, slightly below the maximal behaviour.
Naively, such a behaviour is not completely unexpected, since a higher point time-ordered correlator can decay faster than a lower point time-ordered correlator. Stated simple, slightly after the dissipation time-scale, a 4-point function φφφφ ∼ φφ φφ ∼ e −2t/t d , while φφφφφφ ∼ φφ φφ φφ ∼ e −3t/t d , assuming that the two-point function decays as φφ ∼ e −t/t d . Therefore, as compared to the 4-point OTOC, the 6-point OTOC can already begin showing a chaotic growth that is closer to the large time behaviour of the corresponding correlator. In general, therefore, an n-point function, for sufficiently large n, can display the maximal Lyapunov growth at very early times. Of course, these statements are only plausible ones, based on the limited evidence that we have gathered for this particular example. To gain some confidence in this plausibility argument, we studied behaviour of the six point OTOC with two fold and three fold contour with same ordering but for a couple of slightly differing choices of the argument. We find that the Lyapunov index is reasonably robust under such variations which further justifies our argument. It will be very interesting to explore this in more detail by considering different arguments of the OTOC as well as considering higher point OTOC correlators.

Conclusion
We have computed the fermion six point function in the SYK model with complex fermions in the presence of a non-vanishing chemical potential. We then took triple short time limit of this correlation function so that it appears as a three point function of fermion bilinears. We show that the three point function of fermion bilinears, for h = 2 modes, have the scaling property of conformal field theory three point function, as is expected as a generalisation of the results of [41] to the complex fermion case. Like in [41], we find that the contribution of the contact three point graphs in the large q limit is subleading compared to that of the planar graphs.
We also compute three point function of fermion bilinears for the h = 2 mode. This mode is known to break the conformal invariance of the SYK model, both spontaneously as well as explicitly. This mode is known to exhibit chaotic behaviour with the Lyapunov exponent λ L that saturates the chaos bound. The three point function of bilinears in this case has a behaviour different from those of the conformal, i.e., h = 2 modes. In this case we find that in the large q limit, the contribution of the planar graphs is subleading compared to the contact graphs.
We have also explored the out-of-time-order physics of the associated six-point correlation function, upon taking the analytic continuation of the Euclidean result. As we have emphasized before, our approach is expected to extract an early-time Lyapunov growth in the corresponding OTOC. It is rather interesting that, for a higher point OTOC (e.g. the 6-point one), this early-time Lyapunov growth already features closely a maximal growth, while the 4-point OTOC already hints at the chaotic growth behaviour with a close to maximal Lyapunov behaviour. While there is no apparent inconsistency with the "maximally braided, k-OTO" correlation function studied in [42], it raises the possibility that the late-time Lyapunov growth may already be pragmatically distillable at early-time physics, provided one explores an n-point OTOC, for sufficiently large n. Evidently, an explicit calculation of an n-point OTOC, for larger and larger values of n becomes more and more difficult. So, one is not able to bypass the difficulty of understanding large-time physics, by simply computing a sufficiently high-point OTOC at early-times. These points deserve a better and more thorough scrutiny, and perhaps a natural point to explore is to consider the formal limit of n → ∞. We hope to come back to this aspect in future.
As a further future direction to explore, coming back to the model at hand, since the couplings of the SYK model are chosen from random gaussian distributions, it is also tempting to ask if one can apply techniques of stochastic quantisation to reconstruct the bulk description. We hope to report on this soon.
their constructive remarks which helped improve the results and presentation considerably.
with the following definitions: