Chaos in the Fishnet

We consider the computation of out-of-time-ordered correlators (OTOCs) in the fishnet theories, with a mass term added. These fields theories are not unitary. We compute the growth exponent, in the planar limit, at any value of the coupling and show that the model exhibits chaos. At strong coupling the growth exponent violates the Maldacena-Shenker-Stanford bound. We also consider the mass deformed versions of the six dimensional honeycomb theories, which can also be solved in the planar limit. The honeycomb theory shows a very similar behavior to that exhibited by the fishnet theory.


Introduction
Chaos in many-body quantum systems can be probed using non-time-ordered four point functions [1,2,3,4,5,6,7]. In this language, the butterfly effect is the statement that for rather general operators W and V , the thermal expectation value of the square of the commutator becomes large at late time. For simple operators in large N systems, there is a long period of exponential growth [5] and generically we expect The rate λ L defined by this exponential quantifies the strength of chaos and remarkably, it is governed by a bound λ L ≤ 2π β [8]. The exquisite paper [7] evaluated λ L in a weakly coupled large N quantum field theory, at finite temperature. Our study is heavily influenced by [7] and we are repeating the same analysis for a much simpler class of theories. For this reason it is useful to review salient parts of [7]. The model considered is of a single hermitian matrix Φ ab ( x, t) of mass m interacting through a Tr(Φ 4 ) self coupling. The observable considered is a color and spatially-averaged version of the squared commutator (we indicate the color sums explicitly and use Tr for the thermal average computed with thermal density matrix ρ at inverse temperature β) By splitting ρ into two square root factors, we place the two commutators on opposite sides of the thermal circle, which nicely regulates some divergences without affecting the value of λ L . The computation simplifies at large N since only the planar diagrams contribute. Expanding, the two commutators gives four terms with each computed by a particular analytic continuation of the Euclidean correlator. Each term can also be represented with path integral contours in complex time, where some real-time folds are appended to the Euclidean thermal circle. For each term we generate the perturbation series by expanding the interaction vertex, integrating each vertex along the contour and then applying Wick's Theorem with contour-ordered propagators [7]. In [7] the region of integration for the interaction vertices is restricted to the real-time folds. The rationale for this restriction is that the integral over the thermal circle corrects the thermal state. Although these corrections are important for getting the exact C(t), the claim is that they do not affect the spectrum of growth exponents. This is intuitively convincing and further, these two simplifications are valid at any coupling. We will employ the same simplifications in our study.
To proceed [7] makes one more approximation. To compute λ L to leading order in the coupling, only the fastest-growing function of time at each order is retained. This simplification is dramatic: it restricts the diagrams that are to be summed to the dressed ladder diagrams. The rate of growth of the sum of ladder diagrams becomes the problem of finding the largest eigenvalue of a one-dimensional integral equation. This can be diagonalized numerically. The resulting λ L can be determined very explicitly. One point that we have taken note of, is that λ L is proportional to m −1 indicating that the important degrees of freedom are the highly populated, frequently colliding low energy quanta with E ∼ m. The dominant effect does not come from the thermal scale quanta as one might have thought. For the systems we study we find a divergent λ L as m → 0 suggesting a similar behavior.
Our goal is to carry out the same analysis, but for the fishnet model deformed with a mass term. We will also ultimately need to resort to numerical methods and the mass deformation avoids difficulties from the m −1 dependence of λ L . The fishnet family of models, proposed in [9] are non-unitary, non-supersymmetric CFTs obtained by taking a special double scaling limit of γ-deformed N = 4 super Yang-Mills theory. This limit combines a weak coupling limit (the 't Hooft coupling is sent to zero) with an infinite imaginary twist, γ j → i∞ leaving the theory with three finite effective couplings ξ j = λe iγ j /2 . The limit decouples the gauge fields and the gaugino so that the model has three complex scalars and three complex fermions with a certain chiral interaction. This chiral interaction significantly reduces the number of Feynman diagrams that can be drawn, so that at large N there is often only a single diagram at a given order in the 't Hooft coupling! In this way the double scaling limit significantly simplifies computations of interesting physical properties [10,11,12,13]. The fishnet was further studied in [14] using asymptotic Bethe ansatz methods. In all of these results the conformal symmetry plays a crucial role. We have spoiled this symmetry with our mass deformation. Nonetheless, the deformed fishnet continues to provide a beautiful model because there are still very few Feynman diagrams that contribute. Indeed, we will again find that the computation of the OTOC reduces to summing ladders which occurs because of the specific structure of the fishnet vertex and not because we assume weak coupling. In this way, the fishnet provides a simple solvable toy model that exhibits chaotic dynamics. This is discussed in detail in section 2. Our numerical results for λ L , developed in section 3, show the behavior of λ L as a function of mass for 0.001 < mβ < 0.05. For example, we have at mβ = 0.01. Consequently, as the 't Hooft coupling is increased λ L will violate the bound of [8]. Recall that the fishnet theory is not unitary. The authors of [8] suggest that unitarity, analyticity and causality are the crucial assumptions necessary to prove the bound. They used Hamiltonian systems for which unitarity and causality are manifest. Our result supports the idea that unitarity is an important ingredient in the proof of [8].
The fishnet theory has close cousins in both 3 and 6 dimensions. We consider the six dimensional theory in section 4, again reducing the problem of determining λ L to finding the largest eigenvalue of a one-dimensional integral equation. We do this numerically in section 5. Our numerical results again show that the chaos bound is violated at strong coupling.

Fishnet OTOC
The model that we study is described by the Lagrangian density The fishnet CFT would have 3 complex adjoint scalars. For our purposes it is enough to study a model with two complex scalars so, for simplicity, we will do this. The model is not unitary because the interaction term is not hermitian. Our first observation is that there are no planar diagrams that correct the two point function. For example, the one loop mass correction and two loop wave function correction are both non-planar graphs, thanks to the specific form of the quartic fishnet interaction. See Figure 1 where the one loop mass correction graph is shown. Thus, at large N and at any value of the coupling, the two point functions are given by their free field values. As explained in [7] both the retarded Green's function and the Wightman functions are needed. These are given by Next consider the large N limit of the four point function 1 . The only diagrams that contribute are the ladder diagrams and, as we have just explained, none of the rungs or side rails are dressed at large N . Consequently, the sum of ladder diagrams satisfies the integral equation where the kernel of the integral equation is 2 and where in (2.3) the function f (ω, p) is defined by the Fourier transform of the double commutator The integral equation (2.3) is exact at large N , not restricted to weak coupling. This integral equation is very similar to the equation obtained in [7]. The only difference is that in [7] the two point function is dressed with a two loop self energy. The imaginary part of this self energy leads to an exponential decay of the two point function reflecting the fact that collisions with thermal excitations can knock the particle into a different momentum state so that it has a finite lifetime. As we explained above, this process vanishes at large N , thanks to the form of the fishnet interaction. This correction drives thermalization, so that the planar fishnet theory does not thermalize. For the exponential growth that we are interested in, we need the late time behavior of the correlator. This is determined by a double pole whose contribution is isolated by replacing [7] Now, making the ansatz and (after Fourier transforming ω) f (t, p) = e λ L t f (p) we find We want to do the integral over d 3 k in (2.8). From the explicit form of R(·) we know that R(k ± ) only depends on | p|, | k| and | p − k|. Changing to spherical coordinates, and using the symmetry of the problem, which implies that f (p) = f (| p|), we have where the kernel in the above equation is This defines an eigenvalue problem for the exponent λ L and is the equation that we study numerically in the next section.
There is a final point that deserves discussion. The Lagrangian (2.1) should be supplemented with extra double-trace interactions that are induced by quantum corrections [17,18]. This double-trace renormalization is a leading effect at large N and so can contribute in the planar limit. The possible double trace terms are For the four point correlator computed by C(t) we only need to consider ξ 1 . The contribution of similar double trace interactions, to four point functions, was considered in [15]. In this work integral equations obeyed by the exact four point functions were obtained by introducing suitable graph-building operators. This analysis ignores the contribution of the double trace interactions. A perturbative analysis which adds the double-trace vertices was then considered. These double trace interactions are needed to remove UV divergences introduced by some of the graphs that contribute. The conclusion [15] is that in the end the exact finite coupling solutions are not modified by the double trace interactions 3 .
We would expect these double trace contributions also play an important role in removing UV divergences in the mass deformed theory. However, in the deep UV we expect the mass deformed theory and the massless theory behave in the same way, so it is reasonable to expect the double trace terms would simply remove UV divergences and nothing more. As described in the introduction, by splitting ρ into two square root factors we place the two commutators on opposite sides of the thermal circle, which nicely regulates the potential UV divergences. In this case there are no UV divergences to regulate and we have simply ignored the double trace terms.

Numerical Results
To solve (2.9) we will perform the integral in (2.10) numerically and discretize p and k to obtain a matrix equation. To perform the integral in (2.10) it is helpful to change integration variable from θ to y = k 2 + p 2 − 2kp cos θ. It is convenient to change integration variables before we discretize (2.9). We set so that we map the range 0 ≤ k < ∞ to the range 0 ≤ u < 1. The eigenvalue problem for the exponent then becomes where, after a suitable redefinition of the eigenfunctions f (u), the kernelK(u, v) is symmetric so that we are guaranteed that it can be diagonalized. The largest eigenvalue of the matrix obtained fromK(u, v) after we discretize is the growth exponent λ L .
Numerical convergence is rapid and the numerics converges after discretizing u and v uniformly with about 200 points. Discretizing with a finer lattice does not change the value obtained for λ L . All eigenvalues are positive, in contrast to the results of [7] where a continuum of negative eigenvalues were found. The difference is that for the model of [7] the retarded propagators receive a self energy correction, with an imaginary part at two loops. This imaginary part leads to exponential decay of the two point correlation functions indicating a finite lifetime of single particle states. This is because at finite temperature a particle can be knocked into a different momentum state through collision with a thermal excitation. Our numerical results are given in Figure 2.

Honeycomb OTOC
We study a cubic theory in d = 6 dimensions. The coupling constant of the cubic theory is dimensionless so the massless theory is classically conformal. Following our analysis of the fishnet theory, we will again add a mass term 4 . The model has vertices as shown in Figure  3. The corresponding Lagrangian density is (green for X, red for Y and yellow for Z) The Lagrangian is not hermitian. We have used different couplings for the two different vertices. The theory appears to be non-perturbatively unstable, but we will study it in perturbation theory around the point X = Y = Z = 0. In the massless limit it defines a solvable conformal field theory, realizing Zamalodchikov's integrable honeycomb [20]. Specifically, we study the planar limit obtained by taking N → ∞ holding λ 1 = g 1 √ N and λ 2 = g 2 √ N fixed.
In the planar limit the two point functions Tr(X(x 1 )X † (x 2 )) , Tr(Y (x 1 )Y † (x 2 )) and Tr(Z(x 1 )Z † (x 2 )) are not corrected, so we can use the free field two point functions. Further the three point functions Tr(X(x 1 )Y † (x 2 )Z † (x 3 )) and Tr(X † (x 1 )Y (x 2 )Z(x 3 )) are not corrected, so that the β functions for these couplings vanish at large N . Finally, a simplification The planar four point function is again given by summing ladder diagrams. The sum of ladders (again denoted f (ω, p)) obeys the following equation For the honeycomb theory the rung function is Notice that, in contrast to the fishnet problem, the rung function here does depend on both ω and p. The diagrammatic version of equation (4.2) is given in Figure 4.
The late time behavior is then given by (make use of (2.6)) which can be transformed to time space to give Notice that we have evaluated the rung function at ω = 0 which can be done since we only want the leading late time behavior. Now, as above, making the ansatz and Fourier transforming from ω to t, we find where k ± = (±ω k , k). Here all vectors are 5 dimensional. Making the ansatz f (t, p) = e λ L t f (p) we now find The largest eigenvalue of this eigenproblem tells us about chaos in the honeycomb model.

Numerical Results
It is possible to do 5 of the 6 integrals over l appearing in the rung function (4.3) analytically. The final integral must be dealt with numerically. First, eachG propagator comes with a delta function. The integral over l 0 can be dome immediately using one of the delta functions. We then choose a coordinate system for the spatial components of l such that p = (p, 0, 0, 0, 0) k = (k cos θ, k sin θ, 0, 0, 0) (5.1) with θ the angle between p and k. Using the variables (0 ≤ L ⊥ < ∞, 0 ≤ φ 1 ≤ π, 0 ≤ φ 2 < 2π) l = (l 1 , l 2 , L ⊥ sin(φ 1 ), L ⊥ cos(φ 1 ) sin(φ 2 ), L ⊥ cos(φ 1 ) cos(φ 2 )) (5.2) we can immediately, thanks to rotational invariance, do the integrals over φ 1 and φ 2 . The remaining delta function can be used to do the integral over L ⊥ . Next, changing variables from l 1 and l 2 to v 1 and v 2 where we are able to do the integral over v 1 with the residue theorem. Thus, in the end we are left with a single integral to do The final integral is evaluated approximately. First, rescale the v 2 coordinate as follows Then, in the integrand make the replacement This replacement is accurate for small β and has the advantage that now the final integral of V 2 can be performed analytically. Further, due to the sinh's in the denominator of (5.4) the large momentum contribution is exponentially suppressed. The cosh on the RHS of (5.8) ensures this suppression, which becomes an issue when the eigenvalue problem for λ L is solved. We have verified that this approximation is accurate at small β by comparing with a direct numerical integration of (5.4). Denote this small β rung function by R β (p, k, θ) The eigenvalue problem we need to solve is Thanks to rotational invariance, three of the integrals can be done analytically, leaving an integral over the magnitude | k| = k and an integral over the angle θ between k and p After discretizing, the term in brackets defines a matrix whose eigenvalues give the growth exponent λ L . We consider eigenfunction p 2 f (p) so that discretization gives a symmetric matrix that can be diagonalized. As before we make the transformation (3.1), which maps the infinite k interval onto the interval [0, 1]. The transformation has a Jacobian of 3 (1−u) 2 . Discretize the u and v intervals uniformly with spacing δu to obtain the matrix (k n = k(u n ) and p m = p(v m )) whose largest eigenvalue gives the growth exponent λ L . We have again absorbed a factor into the eigenfunction to ensure a symmetric R mn .
Similar to the case of the fishnet model, we find that R mn only has positive eigenvalues. This is presumably again because there is no correction to the two point function, implying that even at finite temperature our single particle states have an infinite lifetime. The results for the growth exponent, λ L , at small β are given in Figure 5.

Discussion
The fishnet CFT is a fascinating model. It is an interacting non-supersymmetric CFT in four dimensions, which is rare and already makes it interesting. Further, it exhibits integrability which it inherits from its N = 4 super Yang-Mills theory parent. For integrability to work out, conformal invariance is crucial. In this study we have used the same model, but with a mass deformation which obviously ruins conformal invariance. The resulting model is still a beautiful model with enough simplicity that we can solve for correlation functions exactly in the planar limit. In this way, the model has the potential to teach us about phenomena that may be hard to reach at weak coupling.
We have studied certain out of time order four point functions and have found that the model exhibits chaos. Further, it violates the bound of [8]. Our model is not unitary and we interpret the violation of the chaos bound as evidence that unitarity is an important ingredient in the bound of [8]. The fact that we have managed to carry this computation out in the fishnet model suggests it is a good model to keep in mind for other computations. An obvious question is to ask if the spectral form factor can be computed in this model.
One could also ask how the planar results will get corrected once non-planar corrections are included. It is clear that there will be non-planar corrections to the two point functions, and that at two loops, these will have an imaginary part. Thus, a finite lifetime for particles is a non-planar effect in the fishnet (and honeycomb) models. Once included, we expect an exponential decay of the two point functions and hence thermalization.
Finally, we could repeat the computations we have carried out in this study for the three dimensional sextic model too. Based on the results of this study, we expect the three dimensional model exhibits chaos and that again the chaos bound will be violated.