Quantum Chaos, Thermodynamics and Black Hole Microstates in the mass deformed SYK model

We study various aspects of the mass deformed SYK model which can escape the interiors of pure boundary state black holes. SYK boundary states are given by a simple local boundary condition on the Majorana fermions and then evolved in Euclidean time in the SYK Hamiltonian. We study the ground state of this mass deformed SYK model in detail. We also use SYK boundary states as a variational approximation to the ground state of the mass deformed SYK model. We compare variational approximation with the exact ground state results and they showed a good agreement. We also study the time evolution of the mass deformed ground state under the SYK Hamiltonian. We give a gravity interpretation of the mass deformed ground state and its time evolutions. In gravity side, mass deformation gives a way to prepare black hole microstates that are similar to pure boundary state black holes. Escaping protocol on these ground states simply gives a global AdS2 with an IR end of the world brane. We also study the thermodynamics and quantum chaotic properties of this mass deformed SYK model. Interestingly, we do not observe the Hawking Page like phase transition in this model in spite of similarity of the Hamiltonian with eternal traversable wormhole model where we have the phase transition.


Introduction and Summary
Recently there are many progress in understanding the spacetime dynamics like the chaotic behavior in black holes [1,2] or making traversable wormholes [3,4,5]. The Sachdev-Ye-Kitaev (SYK) model [6,7] plays an important role to study these behaviors. The SYK model is an UV complete model that shares the same low energy dynamics with the Jackiw-Teitelboim gravity [8,9,10]. We can solve the SYK model analytically in the large N , low energy limit [11]. We can also study the large N Schwinger-Dyson equation directly numerically [11] and further we can study the model numerically at finite N [12,13,14] by exact diagonalization, which serves a complemental way to analyze the system. Traversable wormholes are made by introducing a direct coupling between two asymptotic boundaries [3,5]. After turning off the coupling between two sides, we can make a two sided black holes [5,15]. Therefore, traversable wormhole protocol also gives a way to prepare the thermofield double states. We can consider the similar question for a single sided system whether we can prepare a "thermofield double" state in a single copy of CFT from some massive deformation. An analog of the thermofield double states for single sided case is the boundary state [16,17], which have a simple entanglement structure. In gravity side, the dual spacetimes are black hole microstates [18,19,20], and by state dependent mass deformation we can reveal the interior of these single sided black holes [21]. Employing the proposal for holographic duals of boundary CFT (BCFT) [22,23], these microstates are modeled by geometry with end of the world (EOW) branes. Now the problem is whether we can prepare the black hole microstates from massive deformations. In the field theory side, the relation between boundary states and gapped ground states are discussed in the literature. Boundary states and the time evolution of them are used to model so called quantum quench [24,25]. The quantum quench 1 is the time dependent process where we suddenly turn off the mass term and evolve the state by the Hamiltonian at critical points. In this context, boundary states are used to approximate the initial gapped ground states. This variational approximation is also used to study the massive deformations of conformal field theory [26] and gives a qualitative picture of the phase diagram.
In the SYK model, an analog of boundary states and their gravity interpretation are proposed in [27]. We can also reveal the interior by the mass deformation [27,28]. Motivated by the above observations on the boundary states, gapped deformation of CFT and their connection to black hole microstates, in this paper we study the ground state of a mass deformed SYK model that also first appeared in [27] and its relation to the SYK boundary state in detail. A merit to consider the SYK model is that we can analyze directly the mass deformed theory itself. This enables us to compare the exact ground state and the variational approximation and check the validity of the variational approximation, which is usually difficult.
Another motivation to study this model is to gain an insight to the gravity counterpart of the chaotic/integrable transition. The relevance of the quantum chaotic properties to the black hole physics was pointed out in [29,30,31]. One way to characterize the quantum chaos is the exponentially growing behavior of the out-of-time-ordered correlation functions (OTOC) [32] in time which is quantified by the quantum Lyapunov exponent. The OTOC was also studied in the SYK model [11] in the large N limit, where the quantum Lyapunov exponent was found to saturate the bound proposed in [2] in the strong coupling limit. This also supports the relation between black holes and the quantum chaos.
Recently the quantum chaotic property was also studied in a variety of deformations of the SYK model [13,14,33,34,35,36,37,38,39,40]. Among these developments in [13] the authors considered a deformation of the SYK Hamiltonian by a random mass term and found that the quantum Lyapunov exponent decreases as the mass parameter is increased. More interestingly, the authors also found that the Lyapunov exponent vanishes at some finite value of mass parameter, where the system exhibits the "chaotic/integrable transition" [13]. The chaotic/integrable transition was also captured by another characterization of the quantum chaos through the level statistics [41,42,43,44] as a sharp transition.
With the relation between quantum chaos and the black holes in one's mind, it may be natural to speculate the gravitational phenomenon dual to the chaotic/integrable transition to be the Hawking-Page transition [45]. Note, however, that these two phenomena do not necessarily happen at the same time. The Hawking-Page like transition can be captured as a first order phase transition through the entropy, or equivalently, the free energy F = − 1 β log Z(β) = − 1 β dE ρ(E) e −βE where ρ(E) is the number density of the states. On the other hand, it is known that the level statistics diagnoses the chaoticity of a quantum system correctly only after the energy spectrum of the system is unfolded [44,43], which essentially subtract the information of ρ(E) from the spectrum. Also, the connection between chaos and the thermalization property of the system implies that one cannot see whether the system is chaotic or not just by looking thermodinamic quantities such as the free energy [46]. Nevertheless, there are some examples [14,34] where we have multiple evidences that they are indeed correlated. It would be interesting to ask what kind of additional properties of a model can relate the chaotic/integrable transition and the Hawking-Page like transition indirectly, which will define an interesting class of theories.

Summary of the paper
We have studied various properties of the mass deformed Hamiltonian which is first proposed in [27]. The model consists from N Majorana fermions ψ i and the Hamiltonian is given by with mean J i 1 ···iq = 0 and variance J 2 We gives an effective action in terms of the collective variables and derive the large N Schwinger-Dyson equation for the mass deformed Hamiltonian. We study this Schwinger-Dyson equation numerically. In the zero temperature case, we found the analytical solution for this Schwinger-Dyson equation in the small mass parameter limit. The diagonal correlation function is related to the SYK correlation function by the conformal transformation. The off diagonal correlation function is also determined in the conformal limit. Using these Euclidean correlation functions in the ground state, we study the several physical observables, which show the non trivial scaling with respects to the mass parameter µ in (1.1) for small µ limit.
We have also used the SYK boundary state, which is also first proposed in [27] and interpreted as a black hole microstate, as a variational approximation for the ground state of the mass deformed Hamiltonian (1.1). We studied this variational approximation both numerically and analytically. We compare the numerical results in exact ground state and in variational approximation and two results show good agreement in the entire mass parameter regime. In the small µ limit, we compare the analytical results in exact ground state and in variational approximation. We found that the scaling with respect to the mass parameter coincide but the proportional constants are different. However the coefficients themselves are also very close. Therefore, the variational approximation is a good approximation but is not a perfect and has an order N difference from the exact ground state.
In section 3.3 we have also computed the large N free energy by solving the Schwinger-Dyson equations at finite temperature numerically. In contrast to the results in the models with a similar Hamiltonian [5,33,34], in our model we have not found a phase transition.
We have also studied our model for finite N in section 4. We have computed the overlap between the boundary state and the true ground state, and have found these two states are close to each other. This result supports the validity of the variational approximation for the ground state in section 3.2. We have also diagnosed the quantum chaoticity of the system by computing the adjacent gap ratio [47,48,49,50] (we explain more detail in section 4.2). As a result we found that the system is chaotic for any value of the deformation parameter and at all energy scale.
In section 5 we have studied the mass deformed model in the large q limit where we can analyze the mass deformed theory analytically beyond the conformal limit. The large q results are consistent with the conformal limit results when the mass parameter µ is small. In the large q limit, the variational approximation actually coincides with the exact ground state in all mass parameter regime. We checked this agreement from the calculation of the ground state energy and the other observables also perfectly coincide. Furthermore, in the large q limit we compute the overlap between exact ground state and the SYK boundary state for the variational approximation. We found the saddle point solution that gives the maximal overlap 1. Even at finite temperature, we can solve the system analytically. We checked analytically that there are no phase transition in the model (1.1). At the order of β ∼ q, the chaos exponents grows from 0 to 2π β that is maximal [2]. Finally, we give a gravity interpretation of the mass deformed ground state. We studied the time evolution of the mass deformed ground state under the SYK Hamiltonian. This is a setup of quantum quench where we suddenly turn off the mass term. We solve this quench problem analytically in the conformal limit, determined the time evolution of the reparametrization mode and found that the system thermalizes. We interpret this dynamics of the reparametrization mode in gravity. The original geometry is interpreted as the global AdS 2 with EOW brane which is static under the global time translation. The quench corresponds to the black hole generation. Therefore, we interpret the mass deformation as a protocol to obtain atypical black hole microstates that are similar to pure boundary state black holes. This is a single sided analog of the preparation of the thermofield double from the two coupling mass deformation [5,51]. We have also applied our gravity interpretation to the escaping interior protocol starting from the mass deformed ground state. The insertion of the mass deformed Hamiltonian before the quench just delays the generation of the black hole and shifts the position of the black hole horizon. As a special case, we can apply the deformed Hamiltonian eternally and we can escape the interior eternally. Actually, this is exactly the global AdS 2 with EOW brane. This is a single sided counterpart of the identification of eternal traversable wormholes and global AdS 2 . We also point out that the mismatch of the sign of spins in the ground state and the mass term leads to the excitation. Too match excitation leads to the huge excitation and we expect that they finally lead to the black hole generation because the system shows the chaotic behavior at high energy as we studied in this paper. Therefore, the escaping protocol is successful only when we choose the mass term in a correct state dependent way.
The organization of this paper is as follows. In section 2, we review the mass deformed SYK model and boundary states those are introduced in [27]. We also review their gravity interpretation as black hole microstates and their time evolution. In section 3, we study the mass deformed SYK model in the large N limit. We derive the large N Schwinger-Dyson equation for the mass deformed SYK model and solve them both analytically and numerically.
In section 4, we study the overlap and the level statistics at finite N and compare with the large N results. In section 5, we study the large q limit of the mass deformed SYK model. We calculate various quantities analytically and checked the consistency with the results in large N , finite q analysis. In section 6, we describe the gravity interpretation of the mass deformed SYK model. In section 7, we discuss some implications of our results and possible future works. In appendix A, we present the derivation of the large N effective action for collective variables. In appendix B, we discuss the numerical method to solve the large N equation of the mass deformed SYK model. In appendix C, we show the detail of large q analysis. In appendix D, we discuss the detail of the large N , q = 4 case.

The model
The main purpose of this paper is to study the following Hamiltonian with mean J i 1 ···iq = 0 and variance J 2 We also defined S k = −2iψ 2k−1 ψ 2k . This Hamiltonian was first introduced in [27]. First we review some properties of this Hamiltonian with its connection to particular pure states in the SYK model. We also review their gravity interpretation of pure states and the evolution under mass deformed Hamiltonian.

A review of the SYK model
The Sachdev-Ye-Kitaev (SYK) model [6,7] is the system of N Majorana fermions ψ i , which obey the anti commutation relation {ψ i , ψ j } = δ ij , with the Hamiltonian with mean J i 1 ···iq = 0 and variance J 2 This Schwinger-Dyson equation comes from the Euclidean action We can solve the equation (2.3) numerically. In the large q limit, we can solve the equation (2.3) analytically. It is also possible to solve at long time (1 J τ 12 N ) regime by ignoring the derivative term ∂ τ in (2.3) as where the scaling dimension ∆ and the coefficient c ∆ is given by When we ignore the term ∂ τ , the equation of motion have a reparametriztion symmetry ). In the low temperature (1 J β N ) and long time (1 J τ 12 ), the thermal correlation function is obtained from the ground state answer (2.5) with the reparametrization f (τ ) = tan π β τ . The low energy reparametriztion symmetry is actually broken by the UV effect. This leading breaking term is given by the Schwarzian action [11] (2.7) The constant α S can be determined numerically. For example, α S ≈ 0.00709 for q = 4, α S ≈ 0.00403 for q = 6 and α S ≈ 0.00257 for q = 8 [11]. At large q, α S goes as α S ∼ 1 4q 2 .

A review of pure states and mass deformation of the SYK model
In the SYK model we can also study the real time evolution of particular pure states. To define such states, first we define a set of spin operators from Majorana fermion operators. They are defined as They satisfy S 2 k = 1, which means that the eigenvalues of S k 's are ±1. Moreover, they are mutually commuting with each others [S k , S k ] = 0. Therefore, we can consider the simultaneous eigenstates of the spin operators S k 's: where s = (s 1 , · · · , sN 2 ) is a set of eigenvalues. These 2 N 2 states span the SYK Hilbert spaces. This condition can also be written as (2.10) We can produce lower energy states by including the Euclidean evolution |B s (β) = e − β 2 H SY K |B s . We can interpret the (2.10) as a transparent boundary condition between fermion field ψ 2k−1 and ψ 2k in the path integral language. Because the states |B s form a basis, the average of the correlators over all choices of s k reproduces the thermal ensemble exactly: This is a true statement for any operator O.
In the large N limit, the model possesses an emergent O(N ) symmetry. This O(N ) symmetry includes an element f 1 that flips the sign of ψ 2 . Similarly, there are elements f k that flips the sign of ψ 2k . Each of f k 's also flips the sign of the spin operator S k . This flip element f k maps the |B s to other state |B s where s is given by the flip of s k from s. Therefore, the norms B s (β)|B s (β) have the same value in all |B s states because of this emergent symmetry. On the other hand, we saw in (2.11) that the average over all the |B s is equivalent to the thermal one. Therefore, the norms of |B s is equal to the thermal partition function in the leading of the 1/N expansion: Similarly, the two point functions such as ψ 1 (τ )ψ 1 (τ ) or ψ 2 (τ )ψ 2 (τ ) are also individually invariant under the flip groups. They are called diagonal corerlators [27]. They also become the same with the thermal correaltors in the large N limit: we obtain Now, the off diagonal correlator times the boundary condition s k becomes a 4 point function with 2 ψ 1 's and 2 ψ 2 's. Therefore, in the large N limit these are flip group invariant correlation function. Then, We also know that in the large N limit the normalization factor becomes B s (β)|B s (β) = 2 − N 2 Tr(e −βH SY K ), and therefore the off diagonal correlator becomes Tr(e −βH SY K ) .
(2.17) Further because we are taking the large N limit, four point function factorizes to the product of 2 point functions and the off diagonal correlator becomes Here we obtain another minus sign because we need contract ψ 2k 's and to do that we need to exchange the order of ψ 2k−1 and ψ 2k .
We can think of the state |B s (β) as a state after projection measurement of thermofield double state [52,53]: (2.19)

conformal limit and the symmetry of correlation functions
In the conformal limit (βJ 1, J |τ − τ | 1), the correlation function becomes Especially, by analytically continuing to real time τ → it, we obtain Therefore, for example the spin operator expectation value S k (t) = −2iG off (t, t) is In this coordinate, it is manifest that the translation τ P → τ P + c is a symmetry of both of the diagonal correlator and the off diagonal correlator. This is the same symmetry of the Poincare patch in AdS 2 . Later we consider the gravity setup with similar symmetry.

evolution under the mass deformed Hamiltonian
Now we consider the evolution under the mass deformed Hamiltonian: In the low energy limit and small µ limit, we can treat this deformation as The Schwartzian action is The term H(t) is evaluated as .
(2.30) Therefore, the coupling to the reparametrization mode becomes Then, we obtain a Lagrangian for f (t): We can write the Schwartzian term using a Lagrange multiplier λ(t) as When we integrate over λ(t), this impose the condition φ(t) = log f and the action reduces to the original one. By introducingη = µ(c ∆ ) 2 J α S , the low energy action becomes On the other hand, the initial condition we consider is set by the Euclidean evolution with the Euclidean action Here we putη = 0 because the Euclidean evolution is given by the Hamiltonian without mass deformation. The solution we are interested in is The equation of motion for λ(t) gives (2.37) The EOM for f gives λ (τ ) = 0. Now we consider the time evolution with the initial condition φ (0) = 0, e φ(t=0) = π 2 J 2 β 2 . Because the equation of motion for f implies that λ(t) is constant, we can set λ = −4J . The Lagrangian is now Therefore, the evolution is simply given by the motion of a particle with a potential This potential crosses 0 at φ = φ × that is given by which gives e φm = (η∆) The Lorentzian dynamics is simply described by the particle motion under this potential with the initial condition A schematic form of the potential is described in Fig. 1. When the deformation parameterη V ( ) ⇥ m Figure 1: A schematic form of the potential for a particle φ(t).
satisfies φ 0 < φ × , the motion of the particle is confined in a finite region. In the SYK model, this especially means that the expectation value S k (t) does not decay and the system does not thermalize. Especially, when φ 0 = φ m , the particle sits on the bottom of the potential and does not oscillate. In this case, using f = e φ the expectation value of the spin operator becomes (2.46) Therefore, the time evolution keeps the expectation value of the spin operator and prevents thermalization. φ 0 = φ × gives a relation between β and µ, which becomes π β(µ)J 2 = (η∆)

the large q limit
For later purpose, we also consider the large q limit of the pure states |B s (β) . The correlator is approximated as The correlation function in the large q limit becomes . This solution can also be written as (2.53) 2 In the notation of [11], we can write the correlation functions as where v ∈ [0, 1] and v satisfies πv The relation with that in our paper is given byα = πv β anď γ = π 2 − πv 2 .

Gravity interpretation of pure states
According to [27] here we consider the gravity configuration that have features in common with the SYK setup. Currently we do not know the precise dual gravity theory of the SYK model. However, the Nearly-AdS 2 gravity has some features in common with the low energy limit of the SYK model. Especially, they share the same low energy theory that is described by the Schwarzian action [10,11]. Therefore, we consider the gravity setup that is similar to the SYK pure states.
In Euclidean signature, the diagonal correlator is the same with the thermal correlator. This is interpreted as the Euclidean black hole or hyperbolic disc H 2 and we imagine that there is a boundary at some finite but very large circle [10]. The difference is the existence of the special point P that corresponds to the insertion of projection operator |B s B s |. Imagining the existence of N bulk fields, this is interpreted as the boundary condition that relates the bulk fields in pairs like ψ 2k−1 = is k ψ 2k at the point P . Except P we impose the same, standard boundary conditions with the thermal case. Other property is the symmetry of the correlation function. We saw that both of diagonal and off diagonal correlation function have the symmetry of Poincare patch in AdS 2 , where the metric is In this coordinate, the special point is sent to infinity τ P = ±∞ and z = ∞. In summary, the Euclidean gravity configuration is the Euclidean black hole with a special point P with boundary conditions on the bulk fields on this point, see Fig. 2. In nearly AdS 2 setup, we interpret this as the special point at large z.
Next, we consider the Lorentzian continuation. The AdS 2 metric in Poincare coordinate is where we defined x ± = z ± t P . Because of the Poincare time translation symmetry of the SYK correlation function, we are interested in the Lorentzian geometry with this symmetry. Especially, the boundary condition at special point should be invariant under the Poincare time translation. This is interpreted as the end of the world line at large z with the same boundary condition with that on the special point P . We can think of this end of the world (EOW) brane as a shock wave that is created by the projection measurement on the left of the thermofield double state and falling to the bulk of AdS 2 spacetime [54], see (2.24). This corresponds to the Rindler Patch. The coordinate transformation t P = f (t) is extended to the bulk by x ± = f (y ± ) where x ± = z ± t P and y ± = X ± t R with the radial direction X in Rindler patch.
In summary, the Lorentzian configuration consists from the AdS 2 geometry with the end of spacetime at large z with the boundary conditions for bulk fields. The cutoff boundary is located on the constant X. The Lorentzian configuration are drawn in Fig. 2.
We can also evolve the SYK model with the mass deformed Hamiltonian (2.1). In this case, the location of physical boundary is oscillating around the constant z and the coordinate covers whole the Poincare patch. Therefore, we can see behind the original horizon in the evolution with deformed Hamiltonian as depicted in Fig. 4. In gravity side, this interaction is interpreted as a change of boundary conditions on the bulk field on AdS boundary. These are interpreted as quantum teleportation [3,4,52], where we measure the left side of TFD state and then apply the measurement dependent time evolution.
The underlying physics of this teleportation protocol is that we try to put each black hole microstate on a ground state of the deformed Hamiltonian to prevent the black hole generation. This is the gravity interpretation of preventing thermalization in the SYK. This essentially depends on how the ground state is close to the ground state of the deformed Hamiltonian and its gap. This motivate us to study the property of the mass deformed Hamiltonian. From next section, we study this Hamiltonian in various methods.

Large N , finite q analysis
In this section, we analyze the Hamiltonian in the large N limit.
Our starting point of the analysis is the Schwinger-Dyson equation for this model with collective degrees of freedom G, Σ [6,11]. According to [5], we also introduce these collective variables for off diagonal component. We study them in Euclidean time. The effective action in the large N limit is Adding mass term ETW brane ETW brane P P Figure 4: The gravity interpretation of evolution in different Hamiltonian. The left figure is the case where we evolve the state by the SYK Hamiltonian. The motion of the UV cutoff particle terminates at the finite Poincare time and correspondingly only the inside of the Rindler patch is visible from the boundary. The right figure is the case where we evolve the state by the mass deformed Hamiltonian. The motion of the UV cutoff particle extends to whole the Poincare time and whole the spacetime within the EOW brane is visible.
The derivation is shown in the appendix A. The Schwinger-Dyson equation arises as the equation of motion for this effective action. They become 3 and with Σ(τ, τ ) = J 2 G(τ, τ ) q−1 . It is also useful to rewrite the Schwinger-Dyson equation in the frequency space. In this representation, we can decouple G off from the diagonal part G, Σ. The equation becomes (3.9) We can determine G off (ω) after solving the equation for diagonal part G, Σ. In the finite temperature, the frequency ω is quantized to the Matsubara frequency ω n = 2π β (n + 1 2 ). Once we solve the Schwinger-Dyson equation, the energy can be calculated from the Green functions in the following way: This is derived from This is the exact relation between the energy and correlation function for the deformed SYK model even before the large N limit or the disorder average.
We can also rewrite the energy using the Schwinger-Dyson equation as Here we used Σ(τ 1 , τ 2 ) = −Σ(τ 2 , τ 1 ). Therefore, in the large N limit the energy becomes This expression is useful when we compute the free energy numerically.

Solving the model in the conformal limit
In this section, we study the ground state of the Hamiltonian H def . We can study the Schwinger-Dyson equation (3.8) numerically where the detail of the numerical calculation is shown in the appendix B. From the numerical analysis for various parameter regions of µ and q, we confirmed that the system has a mass gap above the ground state. We also find that the numerical solution agrees well with the correlation function that is obtained from the reparametrization f (τ ) = tanh(ατ ) of the SYK correlation function in the small µ limit. This is expected since the long time behavior of the SYK model is controlled by the reparametrization (or conformal symmetry) [5,11] and the mass term affects the long time behavior in small µ limit. Therefore, in this section we consider to solve the mass deformed theory using the approximate conformal symmetry of the SYK model.
The diagonal correlation function in the conformal limit is given by α is a function of µ with α(µ = 0) = 0, which we will determine later. The Fourier transformation of the conformal limit correlation function becomes We can easily confirm that in the limit ω α these reduce to the conformal limit of the SYK ground state correlation function This says that only the low frequency part ω α is affected by the mass term, as we expected. This also implies that the solution satisfies the Schwinger-Dyson equation (3.9) in the regime α ω J where we can ignore the mass term µ and the UV term ∂ τ . Now, we solve the Schwinger-Dyson equation at ω α. In this regime, the G c (ω) and Σ c (ω) are linear in ω. However, the slope is very large and we can ignore the first term ω in ω + Σ(ω) in (3.9). Therefore, we can approximate the Schwinger-Dyson equation for diagonal part as . (3.20) Because Σ(ω) is small in ω α, we can also ignore the first term Σ(ω). Then we solve the Schwinger-Dyson equation in the leading of ω expansion by inserting the expression for G c (ω) and Σ c (ω) (3.16) and (3.17). When we expand as G c (ω) = g c (α)ω + · · · and Σ c (ω) = σ c (α)ω + · · · , the Schwinger-Dyson equation gives This determines α as a function of µ as (3.23) The power of µ is given by 1 1−2∆ , which is always larger than 1. Therefore, in the low energy limit the physical mass gap is much smaller than the naive mass gap µ. This is in contrast with the two coupled SYK model [5] where the physical mass gap is much greater than the naive gap µ. We also compute the mass gap numerically and for small µ the numerics agrees with the conformal limit result (3.23).
Once we determine the conformal limit of the diagonal correlation functions, we can also determine the off diagonal correlation function. It is convenient to rewrite the Schwinger-Dyson equation as . (3.24) In the conformal limit, we can ignore the ω in the denominator and approximate G, Σ by the conformal limit G c (ω), Σ c (ω). Therefore, G off (ω) becomes  Figure 5: The plot of the mass gap E gap , which is defined as the exponential decay rate G(τ ) ∼ e −Egapτ of the correlation functions, for q = 4, J = 1 case. In the conformal limit, the mass gap is given by E gap = 2α∆. For small µ, the result in the conformal limit agrees with the numerics well.
The Euclidean time off diagonal correlator is obtained by the inverse Fourier transformation of G off (ω). This inverse Fourier transformation becomes We compare the conformal limit and the exact numerical solution for the Schwinger-Dyson equation in Fig. 6 and they show good agreements.
The τ = 0 value of the off diagonal correlator gives the expectation value of the spin operator S k = −2iψ 2k−1 ψ 2k . In the conformal limit, this becomes 4 Using G off (0), we can calculate the ground state energy: The first relation comes from the relation for the free energy 1 N ∂(βF ) ∂µ = iβ 2 G off (0) and specialize this relation to the ground state β → ∞. By integrating this differential equation, we obtain the ground state energy as where E 0 is the ground state energy of the SYK model. Using the relation H def = H SY K +µH M , we can also compute the expectation value of the SYK Hamiltonian under the ground state of the deformed Hamiltonian as The |G s (µ) has larger energy than the SYK ground state and the energy expectation value of |G s (µ) does not depend on s. Therefore we can prepare 2 N 2 (=dimension of the SYK Hilbert space) states from the mass deformation with the same energy expectation value.

variational approximation for the ground state
To study how the SYK "black hole microstate" is close to the ground state of the deformed Hamiltonian, we apply the variational method for the deformed Hamiltonian by the microstate |B s (β) . This is an SYK analog of variational approximation by smeared boundary states for mass deformations of (1 + 1)d CFT [26].
For variational approximation, we need to evaluate the mass deformed Hamiltonian in the microstate |B s (β) . Here we use the same collection of spins s = {s 1 , · · · , sN we can compute the expectation value of the mass deformed Hamiltonian H SY K + H M Bs as There correlation functions are evaluated in the state |B s (β) . Using the equation (2.17) and (2.18), we can represent this expectation value completely in terms of the SYK thermal correlation function: The first term is the thermal energy in the SYK model [11]: As usual, we minimize the energy evaluated on the trial wavefunction (3.34), to achieve the best approximation for ground state energy.

variational approximation in conformal limit
In the low energy limit, the partition function have the expansion [11] log is the specific heat of the SYK model and E 0 , S 0 are the ground state energy and the zero temperature entropy in the SYK model that is not calculated analytically. Therefore the energy expectation value becomes On the other hand, at low energy 2∆ . Therefore, the expectation value of the deformation term becomes Therefore, the total variational energy is Here we put e φ 0 = π 2 J 2 β 2 andη = µ(c ∆ ) 2 J α S . We should note that this potential is exactly the same with (2.41). The derivative becomes β∂ β = −2∂ φ 0 and the minimal value of the variational energy is the minimal value of the potential V . This potential has a unique minimal that is given by Therefore, the relation between β and µ becomes  Figure 7: Comparing the µ dependence T (µ) = β(µ) −1 from variational method for low energy approximation and exact numerical calculation for q = 4 and q = 6 case. The conformal limit is given by (3.41).

The variational energy becomes
Using the variational wave function, we can compute several physical observables. For example, we can compute the spin operator expectation value S k = −2i ψ 2k−1 ψ 2k , which is essentially the off diagonal correlation function at τ = 0. The half of the spin operator expectation value becomes Another observable we can compute is the energy of the SYK Hamiltonian H SY K that gives the energy of the ground state of the deformed Hamiltonian as an excited state of the SYK Hamiltonian. This becomes As a consistency check, we also solve the minimization condition for the trial energy (3.34) using the numerical solution for thermal SYK correlation functions. The comparison of numerics and the analytical results in conformal limit is shown in Fig. 7. Even Beyond the conformal limit, we can study both of the variational approximation and the ground state numerically. Especially, we can compare both results in the whole parameter region. In Fig. 8, we show the numerical results for the spin operator expectation value S k , ground state energy E 0 (µ) and energy in the SYK Hamiltonian H SY K for both of the exact ground state |G s (µ) and variational approximation |B s (β(µ)) . We found that these observables in |G s (µ) and |B s (β(µ)) are very close and |B s (β) is a good approximation for the ground state. We also checked that the true ground state energy never goes beyond that in the variational approximation, which is expected.

Comparison of variational approximation and Exact ground state
In the conformal limit, we have analytic expression both for the exact ground state and the variational approximation. By comparing the results, we can find that the variational approximation reproduce the correct scaling with respect to the mass parameter µ. On the other hand, the coefficients are different. This means that the variational approximation is not perfect even in the small µ limit. This is in contrast with the two coupled SYK model [5] where the observables in the exact ground state and the thermofield double state perfectly agree in the small mass parameter limit.
However, in the large q limit, the observables in |G s (µ) perfectly agrees with those in |B s (β(µ)) . Actually, we can study the large q limit analytically in the whole parameter regime and we can confirm that the variational approximation is perfect in any µ as we will see later.

Thermodynamics of the deformed SYK model
In this section we study the thermodynamic property of the deformed SYK model (2.1). In the complex SYK model with a similar deformation, an interesting phase structure was found [33,34] through the analysis of the large N free energy F N = − 1 N β log Z: the first order phase transition in µ-T plane 5 and the disappearance of the phase transition above some critical values of µ and the temperature T . The similar phase structure was also found in the two coupled real SYK model with equal random couplings [5]. It would be natural to expect a similar phase structure also in our setup.
The large N free energy can be evaluated by solving the Schwinger-Dyson equations (3.5), (3.8) and then evaluating the partition function on that solution. As we are interested in the phase structure at finite (µ, T ), we solve (3.5), (3.8) directly without any further approximation and numerically by discretizing τ direction. See appendix B for detail. The Schwinger-Dyson equations are discretized as (B.13) and the free energy is evaluated through (B.15). Here we have chosen the discretization parameter as τ = βm 2Λ (m = 1, 2, · · · , 2Λ) with Λ = 10 6 . For each µ, we have first solved the Schwinger-Dyson equation for T = 0.3 numerically by an iterative method [11, appendix G] with initial values for G and Σ = J 2 G q−1 chosen as G n = i ωn (ω n = 2π β (n + 1 2 )). Then we have decreased the temperature slowly by solving the equation for the temperature T − ∆T with the initial condition chosen as the solution obtained for the temperature T , with ∆T = 5 × 10 −5 . Once we reach a sufficiently small temperature, we solve the Schwinger-Dyson equation again by slowly increasing the temperature in the similar way. This recursive technique is similar to the technique employed in [5,33,34]. If we find two The results are summarized in Fig. 9. We find that the free energy for each µ interpolates two extreme behaviors: F = const. (i.e., gapped) for low temperature and F ≈ F SYK at high temperature, which is consistent with the structure of the deformed Hamiltonian (2.1). From the observations [5,33,34] we suspected that the system exhibits a first order phase transition in the intermediate temperature (for example, T ∼ 0.04 for µ = 0.2). However, we have not observed the aformentioned hysteretic behavior which would indicate the first order phase transition.
We further examine the presence of the second order phase transition by calculating the large N specific heat which would diverge at the second order phase transition point. See Fig. 10. Though the specific heat exhibits a peak at some temperature in the intermediate regime, we find that the peak is finite and smooth.
From these result we conclude that our model exhibits neither the first order phase transition nor the second order phase transition. 6 This result is rather surprising and we discuss possible  .1), here the horizontal axis the temperature T . Note that the universal increasing behavior at T ≈ 0 is a numerical artifact due to the fact that the numerical UV cutoff |ω n | < 2πΛ β is not large enough.
explanation in section 7.

Finite N analysis of the model
In this section, we study the mass deformed Hamiltonian (2.1) at finite N . We focus on the case with q = 4 and J = 1 of this model.
Since the canonical anti-commutation relation of ψ i , {ψ i , ψ j } = δ ij can be realized by the Gamma matrices Γ i as ψ i = 1 √ 2 Γ i , the Hamiltonian H def (2.1) for finite N is written as the following 2 N/2 × 2 N/2 matrix with J ijk random coupling chosen out of Gaussian distribution with the mean J ijk = 0 and the variance J 2   .1) and N = 30, with a single realization of J ijk . We observe that the shape of the eigenvalue density around the ground state exhibits a transition around µ ≈ 0.1 from a hard edge to a smooth decay, which is consistent with the behavior of E gap ; for µ 0.1 E gap ∼ µ 2 , which is significantly smaller than E gap ∼ µ.
Note that H def commutes with the following chirality (i.e. fermion number in ψ i ) matrix whose eigenvalues are ±1. Hence with an appropriate choice of basis, H def takes a block diagonal form 2 , regardless of the choice of J ijk . In Fig. 11 we display the eigenvalue density of H def for N = 30 and various values of µ.
p . Though these degeneracies are resolved by H SYK , the levels at different E p are not mixed for a sufficiently large µ, hence we obtain a blob structure.

Overlap
In section 5.2 we have realized that the spin ground state |B (↓,↓,··· ,↓) is a good variational ansatz to realize the true ground state energy of H def after the Euclidean evolution e − β 2 H SYK , with β being the variational parameter. In this section we would like to examine the agreement of these two states more directly, through the overlap of the states where |0 (+) is the ground state of H (+) and |B (↓,↓,··· ,↓) β is defined as with |B s defined in (2.9) and normalized as B s |B s = 1. Here β is chosen for each realization of J ijk such that the overlap (4.4) is maximized.
Note that |B (↓,↓,··· ,↓) β has a definite chirality Γ c |B (↓,↓,··· ,↓) β = +|B (↓,↓,··· ,↓) β for any values J ijk and N . This follows from the fact Γ the rising/lowering operator for S i , together with the following alternative expression of Γ c (4.2) The results are displayed in Fig. 12. For large µ, the Hamiltonian is dominated by H M whose ground state is |B (↓,↓,··· ,↓) , hence the overlap trivially approaches to 1. For small µ, the Hamiltonian is dominated by H SYK . Since the Euclidean evolution with β → ∞ is equivalent to the projection onto the ground state of H (+) SYK , the overlap should again approaches to 1. Note, however, that for N ≡ 4 mod 8 the spectrum of H (+) is two-hold degenerate. The degeneracy is resolved by a small perturbation by H M , and at the leading order in µ the ground state |0 (+) of H def is a certain linear combination of the two ground state of H SYK which is not necessarily the same linear combination obtained by the projection of |B (↓,↓,··· ,↓) . Hence we expect that the overlap is substantially smaller than 1. 8 The results in Fig. 12 are consistent with these expectations. On the other hand, for intermediate values of µ we have found that the overlap is not close to 1 any more even for N ≡ 4 mod 8, and the lowest value around µ = 0.01 significantly decreases as N increases. 7 If one is interested in the overlap between |B (↓,↓,··· ,↓) and the true ground state |0 , one has just to multiply the "probability of |0 to have Γ c = +1" to the results displayed in Fig. 12. Though we do not have an analytic expression, we observe for any N that this probability is almost 1 for µ ≥ 0.5 and not smaller than 0.5 also for the smaller values of µ. Especially the difference between |0 and |0 (+) does not matter when we consider | β B (↓,↓,··· ,↓) |0 (+) | 1 N (figure 13) in the large N limit. 8 Though we do not have a clear argument for this effect, we observe that the value of the overlap approaches some finite value as N increases from N = 12 to N = 28.  Figure 12: The maximized overlap averaged over the realizations of random coupling J ijk as Here the horizontal axis is µ.
Note, however, that as the dimension of the Hilbert space increases, the agreement of two vectors |φ ,|χ in the sense of | φ|χ | ≈ 1 becomes less likely to occur. For example the expectation value of the overlap of two randomly chosen unit vectors in d dimensional space can be evaluated as follows where in the second line we have realized the randomness of |e 1 ,|e 2 as |e 1 = U 1 |e , |e 2 = U 2 |e with random unitary transformations U 1 , U 2 and an arbitrary pair of fixed unit vectors |e ,|e .
In the third line, taking into account that the result is independent of the choice of |e ,|e , we have further replaced |e e| and |e e | with 1 agreement for large µ (µ > 0.2). On the other hand the two results are significantly different (by factor ∼ 100) for the smaller µ. However, it is not necessary to have an agreement in the first place since we have determined β(µ) through the two different quantities. Indeed, though the variational ansatz reproduced the ground stat energy of the deformed Hamiltonian well, there was a discrepancy in another observable | S k | (see Fig. 24; for a possible explanation for the discrepancy, see appendix D). This implies that |B ↓,↓,··· ,↓ β with β(µ) determined by minimizing the energy was actually not so a good approximation to the ground state itself.

Chaotic property
In [14] the authors conjectured that the Hawking-Page like transition of the model [5] is accompanied with the chaotic/integrable transition. Here we would like to test this proposal also for the current setup. In section 3.3 we have found that our model does not exhibits a phase transition in µ or in the temperature T . Hence, if the proposal is correct, our model should not exhibit a chaotic/integrable transition.
As a diagnostics of the quantum chaoticity, in this paper we adopt the level statistics which is relatively easy to study for finite N . It was conjectured that [42] if we quantize a classically chaotic system the fluctuation property of the resulting energy spectrum exhibits the same correlation among different levels as in the random matrix theory. Here the ensemble of the random matrix is determined by the time reversal symmetry of the Hamiltonian of the quantized system. Though a rigorous proof at fully quantum level is still lacking, this conjecture have been verified in various systems [43,44] and also proved at semi-classical level [55,56]. Hence one may use the presence of the RMT-like level correlation conversely as a reasonable definition of the quantum chaos.
Among various ways to characterize the level correlations, here we adopt the following quantity called the adjacent gap ratio [48,47,49,50]: where {E i } is the energy spectrum (E i ≤ E i+1 ) and (· · · ) in the right-hand side stands for the average over the spectrum. This quantity is evaluated for the random matrix theories with various type of the ensemble [47] as well as for the Poisson distribution which corresponds to the non-chaotic systems. By comparing the result obtained from the actual energy spectrum with these known values, one can diagnose whether the systems is chaotic or not.
As the Hamiltonian of our model is trivially separated (4.3) due to the conservation of chirality, the adjacent gap ratio should also be defined separately for the spectrum of each of where the spectrum {E i+1 . The average is taken over J ijk for each fixed i. Here we do not take the average over the spectrum; in this way we can diagnose the chaoticity of our model at each energy scale separately. The results are displayed in figures Fig. 15 and Fig. 16.
The time reversal symmetry of one dimensional fermion systems were studied in [57]. For N = 30, H (±) def has the same time reversal property for both µ = 0 and µ > 0 which corresponds to the Gaussian unitary ensemble (GUE) [57,58], hence we can safely compare our results with the adjacent gap ratio of GUE r GUE = 2 √ 3 π − 1 2 and that for the Poisson distribution r Poisson = 2 log 2 − 1. In contrast to the result obtained in [14], here we find that the adjacent gap ratio is close to r GOE over whole the spectrum, which implies that the system is chaotic for any values of µ and the energy scale (temperature); there are no chaotic/integrable transition. This is consistent with the proposal in [14].

Large N , large q analysis
In the large q limit, we can study the mass deformed SYK model analytically beyond the low energy approximation. In this section we study this limit to confirm the validity of the low energy approximation and the observation by the numerical analysis of finite q model in the region where we do not use the low energy approximation. In the large q limit, the G, Σ action reduces to the Liouville action: with the large q expansion and we also scale µ so thatμ = µq is kept finite in the large q limit. The derivation is shown in the appendix A. At small temperature and the late time of order τ ∼ q, this approximation  is not valid because of the exponential decay of the correlation functions. In this case, we also consider the solution in τ q regime and impose the matching condition between τ q and τ q solutions.

large q limit at zero temperature
At large q limit we can write the correlators as In the mass deformed theory, it is convenient to consider the to scale the mass term µ =μ/q and keepμ when q → ∞. The Schwinger-Dyson equation reduces to the following two equation: with the boundary conditions The general solutions of the equations (5.4) become with constants of the integration α,α, γ,γ. Each boundary condition (5.5) fixes the constants of integration in a following way This means This solution for τ 1 − τ 2 > 0 can be written as Using (3.10), we can compute the ground state energy of the deformed SYK model in the large q limit as follows.
In smallμ limit, we can approximate γ ∼μ 2J 1. In this limit, the ground state energy becomes The first term is the "ground state energy" of the SYK model at large q limit.
Given the ground state correlation function, we can compute the several physical observables again. The spin operator expectation value becomes The SYK energy evaluated on the ground state of the deformed Hamiltonian is (5.14) Inμ → 0 limit, using γ ∼μ 2J the first term becomes the SYK ground state energy. Therefore |G s (µ) serves an excited state of the SYK model with energy higher than the ground state bŷ µ 2q 2 . Inμ → ∞ limit, γ becomes ∞ and the |G s (µ) have the 0 energy in the SYK Hamiltonian, which is expected to the state |B s [27].

variational approximation in the large q limit
We can also study the variational approximation of the ground state of the deformed Hamiltonian by the SYK black hole microstate analytically even beyond the low energy approximation.
In large q limit, the trial energy (3.33) becomes Usingα = J sinγ, this can be rewritten as The R.H.S is monotonic onγ ∈ [0, π 2 ] that runs from 0 to ∞ and this have the unique solution. This is solved as Together with the relation with γ andμ (5.8), we can also write the matching condition as The inverse temperature is given by For smallγ,μ J ≈γ 2 and the temperature β(μ) is approximately On the other hand, for largeμ, we can approximateμ J ≈ 1 π 2 −γ and the temperature is approximated as Note that in the large β limit, using c ∆ ≈ 1 2 , α S ≈ 1 4q 2 and ∆ = 1 q the low energy approximation the low energy approximation for the relation 1 which completely agrees with the smallμ limit of the large q answer.
Using the matching condition (5.20), we find that the exact ground state energy (5.11) and the variational energy (5.16) actually exactly agree up to the order of q −2 . This means that in large q limit the black hole microstate |B s (β) is the same state with the ground state of the deformed Hamiltonian! Later we will confirm this fact by computing the overlap between |B s (β) and |G s (µ) at large q limit using the Liouville action.

large q limit at finite temperature
In this section we consider the large q limit at finite temperature. One motivation is to confirm the absence of the Hawking-Page type phase transition in the mass deformation in this paper at large q limit. In large q limit, Σ vary over a relatively short time, which is of order one. Moreover, (3.6) shows that Σ off is proportional to the delta function. On the other hand, G and G off varies with the time scale of order q. Using these separation of the time scales, we can approximate the convolution (3.4) as follows. Σ(τ ) is an odd function of τ , and we can approximate Σ(τ ) ∼ δ (τ ). Therefore, we can approximate the integral However Σ(τ, τ ) contains the factor 1/q and the equations already contain ∂ τ G and ∂ τ G off , we can ignore the term that contains Σ. Because we are considering the large τ regime, we can also ignore the term δ(τ − τ ) in the right hand side of the Schwinger-Dyson equation (3.6). Therefore, we obtain the equation The finite temperature solution is When we expand them in τ , we obtain In the following, we study the thermodynamical properties of the Hamiltonian H def in the large q limit. We study the inverse temperature regime of order q log q, q, √ q and 1. The derivations are skipped here and shown in appendix C.

Inverse temperature of order β = q log q
In this regime, it is convenient to use the parameter σ = qe −βµ , which is of order one quantity in this temperature regime. In this temperature regime, we can still use the large q expansion G(τ ) = 1 2 (1 + 1 q g(τ ) · · · ) and G off (τ ) = i 2 (1 + 1 q g off (τ ) · · · ) at early time. The solution for τ q becomes and for τ q The thermal energy, thermal free energy and the thermal entropy is We can also rewrite the free energy as whereμ = 2J sinh γ is a function of only µ. In this expression, it is clear that the free energy is the monotonic, smooth function of β in this temperature regime. This means that there are no phase transition in the large q limit. We observe the absence of the phase transition numerically in large N finite q case in section 3.3, and the large q analysis here is consistent with with this observation. This is contrast with the two coupled SYK model [5] where that model have a phase transition in the same temperature regime. The main difference from that model is that here the temperature is the monotonic function of σ. Because of this, we always have one solution for a given temperature and we do not have phase transition.

Inverse temperature of order β = q
In this order, the off diagonal correlator |G off (τ )| is smaller than 1/2 everywhere and we cannot use the same large q expansion for the off diagonal correlator as we did in the last subsection. We can still assume the large q expansion G(τ ) = 1 2 (1 + 1 q g(τ ) + · · · ) for the diagonal correlation function. The correlation function for τ q becomes e g(τ ) = α 2 J 2 sinh 2 (α|τ | + γ) , and for τ q, (5.36) The free energy becomes The first term, which is of order one, is the same with the free energy of the free fermionic oscillator. In the small β limit, this becomes 1 2 log 2 which is the leading of the thermal entropy of the SYK model at large q limit. Therefore, in this regime the entropy is increasing from low entropy regime of order β = q log q. In the large β limit, this can be expanded as which reproduce the free energy (5.33) in the order of β = q log q. In the high temperature limit, we can expand γ and F as As a check, we compare the large q results (5.33) and (5.37) for free energy with the free energy calculated from the numerical solution for the Schwinger-Dyson equation in Fig.18, which shows good agreement.
At this order, we obtain the same results with the [5]. Actually, we found that both models have the same Schwinger-Dyson equation at this order. In higher than this temperature we still have the same equation of motion and we only reproduce the former results in [5], but to make this paper to be self contained, we still continue the finite temperature analysis.

Inverse temperature of order β = √ q
In this regime, we can approximate the off diagonal correlation function G off as which is of order 1 √ q . The diagonal correlation function G is approximated as G(τ ) = 1 2 (1 + g(τ ) + · · · ) everywhere in τ ∈ [0, β] and the equation of motion for g becomes The same equation has also appeared in a different mass deformation of the SYK model [13].
The last term is of order 1/q, which seems to be ignorable. But at the time of order √ q, the other terms are also of the same order. This can be seen clearly after rescaling as Then, the equation of motion becomes The detailed analysis are in appendix C.

The partition function becomes
where h(k) is a function that we have not determined.  Figure 19: Plot of Lyapunov exponent as a function of √ qβµ = √ 2k with exact, small k and large k expansion.
In this regime, the chaos exponents increase from 0 to the maximal value 2π β . When k is large, the chaos exponents λ becomes 45) and for small k the chaos exponents becomes For finite k, we can numerically study the chaos exponent. The plot is shown in Fig.19 and the details are shown in the appendix C.

Temperature of order β = 1
In this limit we can ignore the mass term and we obtain the same physics with the large q SYK model [11]. The free energy becomes where F SY K is the free energy of the SYK model. The chaos exponents are maximal when 1 β √ q and then decrease to 2J in the high temperature regime β 1.

Computing the Overlap at large q limit
In the large q limit, we can compute the overlap using the Liouville on shell action. The strategy is to construct an analog of "Janus" solution [59] in the large q limit, where a similar holographic computation of the overlap is done in [60]. The overlap is represented as with an initial condition |0 , which only changes the normalization constant that should cancel between the numerator and the denominator. We can treat the Euclidean path integral for the overlap as a Euclidean time dependent coupling where τ runs in the range τ ∈ [− β 2 , ∞] and the time dependent coupling µθ(τ ) as depicted in Fig.20. After the disorder average, we can again obtain the effective action for G(τ 1 , τ 2 ), Σ(τ 1 , τ 2 ) variables with τ 1 , τ 2 ∈ [− β 2 , ∞] and the time dependent mass term −iµ τ ). Because this mass term explicitly depend on the Euclidean time τ , we do not have time translation symmetry and the solution depend on two times τ 1 , τ 2 . At τ = − β 2 , the state |B s impose the boundary condition ψ 2k−1 |B s = is k ψ 2k |B s , which leads to the boundary condition We also require that the solution approaches to the ground state solution of the deformed Hamiltonian at τ 1 , τ 2 → ∞. In the large q limit, the effective action reduces to the Liouville action. The difference with the ground state or thermal case is that we do not have time translation symmetry and the Green's functions depend on the two time variables as The field g satisfies the Liouville equation and g off satisfies the free field equation These equation of motion should be satisfied except for the line τ 1 = τ 2 where we impose The two time solutions are locally given by The matching condition from the variational method is e −γ = cosγ. (5.54) This also gives the relation Actually we can find the two time solution for maximal overlap. The solution is I : More compactly, we can write the solution as with the region dependent functions Now we can compute the overlap in the order of 1 q 2 using the Liouville on shell action. We denote S by the on shell action of the Liouville fields for overlap solution, and use S Bs for the on shell action for the state |B s (β) and S Gs for the ground state |G s (µ) . To compute the overlap, it is convenient to rewrite the Liouville action (5.1) using dimensionless coupling βJ and βµ.  with θ i = 2πτ i β for i = 1, 2. We take the derivative of the action S over J withμ fixed and the matching condition β = β(J ,μ). Then, we obtain Here we again used the fact that we can ignore the contribution from the variation of the field g, g off because of the equation of motion. In the third line, we use θ i = 2πτ i β and the equation of motion for g, g off . Now, using the the property of the two time solution we can integrate over τ 2 in the first term and we obtain τ ). (5.63) This means that the derivative of the on shell action only depends on the correlation function on τ 1 = τ 2 line, and especially that does not depend on the region III. Since this two time solution is equal to that of the boundary state |B s (β) in region I and identical to that of the ground state of the deformed Hamiltonian in region II, we obtain Since we can explicitly check that the overlap becomes 1 at J = 0, by integrating the above equation we obtain | B s (β)|G s (µ) | = 1 for general J andμ. Since the Liouville action capture up to 1 q 2 terms in the 1 q expansion, this overlap computation shows that the overlap behaves as e − N q 3 in large q expansion. In fact, we observed from the variational approximation that there is a finite difference between |G s (µ) and |B s (β(µ)) even in small µ regime.

Gravity interpretation
In this section, we consider the gravity interpretation of the mass deformed SYK model. Though we do not know the exact dual gravity of the SYK model, we can consider the similar gravity setup as we did for the microstate state |B s (β) [27]. Here we take the same approach with [27] where we consider the gravity configuration with the same symmetry with our SYK setup. First we consider the ground state |G s (µ) and its time evolution under the SYK Hamiltonian, and then consider the gravity interpretation.

Time evolution under the SYK Hamiltonian
In this section, we consider the time evolution of the ground state |G s (µ) under the SYK Hamiltonian H SY K . We can formulate this time evolution as time dependent mass term H def (u) = H SY K + θ(−u)H M where u is the Lorentzian time. This type of time evolution is called as quantum quench. A different type of quantum quench and black hole formation was studied in [5,15,61,62]. The quantum quench with time dependent mass terms are also studied in quantum field theories [63,64].
We saw that the ground state |G s (µ) has bigger energy expectation value than the ground state and is an excited state of the SYK model. Because of the similarity with the state |B s (β(µ)) , we also expect the similar thermalization for the state |G s (µ) . We solve this time evolution in the low energy limit where the SYK dynamics is governed by the Schwarzian action. For u < 0 with the Lorentzian time u, the reparametrization is given by f (u) = tan(α(µ)u), which is the Lorentzian version of the reparametrization to obtain the ground state correlation function. Then, we couple the reparametrization mode f (u) = tan(α(µ)t(u)) where t(u) is the reparametrization. For u > 0, because of the energy conservation, we impose where E 0 is the ground state energy and − N α S J {f (u), u} gives the energy increase from the ground state [65]. We have already evaluated the right hand side G s (µ)|H SY K |G s (µ) in (3.30) and the above equation is solved as The second equation determines the inverse temperature β in terms of µ 9 . We can also rewrite f (u) = A tanh( π β u+B)+C with three parameters A, B and C. These parameters are fixed by imposing the continuity for f (u) at u = 0 up to the second derivative, which becomes f (0) = 0, f (0) = α(µ) and f (0) = 0. This condition fix the reparametrization to be . Using the reparametrization (6.3), we can study the time evolution G > (u 1 , u 2 ) = G s (µ)|ψ i (u 1 )ψ j (u 2 )|G s (µ) using the reparametrization where ψ i (u) = e iH SY K u ψ i e −iH SY K u . The diagonal correlation function becomes This is exactly the thermal correlation function in Lorentzian time. The time evolution of the spin expectation value can be studied from the off diagonal correlation function as S k (u) = −2is k (t (u)) 2∆ G off (t(u), t(u)), which becomes The spin operator expectation value decays exponentially at late time. Therefore, the system loses the initial simple correlation pattern under the SYK time evolution and thermalizes. The term is close to one because π βJ is very small when µ J . Therefore, the time evolution is very close to that in |B s (β) , which is given in (2.23) 10 .

Gravity interpretation
As it is done in [27], we can consider the similar gravity configuration of our analysis. The ground state |G s (µ) is invariant under the evolution e −iH def t because it is the ground state of the deformed Hamiltonian H def . Because f (τ ) = tanh(ατ ) is the transformation from Poincare coordinate to the global coordinate [5], we expect the time translation symmetry in gravity side where the metric in this coordinate is given by Hence the system is gapped, we also expect the confined geometry where the emergent direction is capped off at some scale. Here we simply use the end of the world (EOW) brane picture on which the geometry terminates [66,67,68,69]. Because of the time translation symmetry the position of EOW branes should be static under the time translation along global time. We imagine that we have N bulk fields and at EOW branes we impose the boundary condition ψ 2k−1 = is k ψ 2k for the bulk fields as we did in the case of |B s (β) states.
When we evolve the ground state |G s (µ) by the SYK Hamiltonian, the system thermalizes. The evolution under the SYK Hamiltonian is given by the reparametrization (6.3). In gravity picture, this reparametrization gives the transformation from the global coordinate to the Rindler coordinate, which only covers a portion of global AdS 2 and has a horizon. Therefore we obtain the single sided black hole geometry with EOW brane from the ground state of the mass deformed Hamiltonian.
We can also interpret the similarity between |G s (µ) and |B s (β) in gravity. The symmetry of |G s (µ) is that in global time whereas the symmetry of |B s (β) is that in Poincare time and EOW branes are static under each symmetry. We can still match the Rindler patch in both geometries. Then, the EOW branes are falling from Rindler observer in a similar way. In this 10 In the two coupled SYK model, similar spin operator is constructed from left and right fermion as Under the decoupled Hamiltonian evolution, this behaves as S i (u) = S i (0) (cosh 2π β u) −2∆ . Though this shows the same exponential decay, the early time behavior is different from (2.23) and (6.5). Matching Rindler Patch ETW brane Figure 22: A cartoon of the gravity configuration. The left is the bulk interpretation of the |B s (β) and the right is that of the |G s (µ) . In the middle picture, we compare two geometries matching the Rindler patch of both geometries. From the Rindler observer, the EOW brane is falling. The Rindler observer feels the similar falling pattern for the EOW brane.
sense, two geometries are similar. Especially, we expect that the state |G s (µ) contains region behind the horizon.
It is also interesting to consider the protocol to escape the black hole interior [21,27] of single sided black holes with the black hole microstate |G s (µ) instead of |B s (β) . When we evolve the system by the SYK Hamiltonian, these correspond to single sided black holes. The escaping protocol [21] corresponds to evolving the ground state by the deformed Hamiltonian H def . We can apply the escaping protocol for finite time T and then turn off the mass term. This corresponds to insert the time evolution by H def before applying the SYK evolution as e −iH SY K t e −iH def T |G s (µ) . Therefore we just delay the black hole formation by inserting global AdS 2 region. When we apply the escaping protocol eternally, we shift the horizon infinitely and finally we obtain the geometry without horizon. This corresponds to the evolution e −iH def t |G s (µ) and as we observed this corresponds to the global AdS 2 patch. Therefore, in this case after eternally escaping the interiors we obtain the global AdS 2 with the EOW brane. The matching of spins s in the state |G s (µ) and those in the escaping protocol e −iH def T is important because the mismatch of the spins gives excited states of the H def . As we saw in the finite temperature analysis of the H def , high energy behavior is similar to that of the SYK model and chaotic. Therefore, when we have mismatch for order N spins, we expect that this  The SYK evolution, which is interpreted as the evolution without any double trace deformation, makes the black hole with EOW branes. We also evolve in backward by the SYK Hamiltonian. Middle: We apply the escaping interior protocol for finite amount of time T and then evolve by the SYK Hamiltonian. This is equivalent to shifting the horizon by insert the global AdS 2 patch. Right: We apply the escaping interior protocol for eternally. As a consequence, the horizons are shifted infinitely away from the original horizon. Finally we recover the global AdS 2 with the EOW brane.
mismatch leads to the black hole formation and failure of the escaping protocol. Therefore the state dependent deformation is important 11 to avoid the black hole generation. In this way, we can clearly understand the escaping protocol starting from the special microstates |G s (µ) .

Similarities and differences compared with Maldacena-Qi model
Because the model is similar to that of the eternal traversable model [5], it is good to compare with that. The Hamiltonian of the eternal traversable model is given by It is also important to choose the correct pair of fermions to make a spin operator.
with J L i 1 i 2 ···iq = J R i 1 i 2 ···iq . Here we introduce two copies of Majorana fermions ψ L i and ψ R i which satisfy the canonical commutation relation.
The similar thing is that both systems are gapped systems. This is natural because in both models we explicitly introduce the mass term in the Lagrangian. Both systems can be analyzed using conformal symmetry and the ground state has the same time translation symmetry that corresponds to the global time in AdS 2 . In the large q limit, the finite temperature behavior beyond the order of β ∼ q is the same with that of Maldacena-Qi two coupled model because we obtain the same equations.
When we consider the gravity interpretation, it is more surprising. In the traversable wormhole case the two side are connected in the deep interior. On the other hand, in our case the geometry is lost at the mass gap scale, which should happen in duals of confining phase [66,67,69,68]. This suggests that we may be able to understand the spacetime connectivity in a similar way to understand the confined geometry.
There are differences even in qualitative levels. The first big difference is the absence of the Hawking-Page like transition. There are many examples of mass deformation of the SYK model, tensor models or matrix models that show the Hawking-Page like transition [5,33,34] in the large N limit and it is surprising that we have not Hawking-Page like transition even at small mass range. We expect that this is reminiscent of the higher spin like nature of the SYK model, which suppress the order of transition.
Another difference is the size of the mass gap in the theory at low energy. In the two coupled SYK model, the physical mass gap is much larger than the parameter µ in the Lagrangian in small µ limit. Therefore the chaos helps to open a gap [51]. On the other hand, in our case the mass gap is much smaller than the naive gap µ. In our model the chaos suppress the mass gap, which seems to be more natural. We expect this is related to the absence of the Hawking-Page like transition. We will revisit this problem in the future [70].

Comparison with the Complex SYK model
It is also good to compare with the complex SYK model [71,72,73] because this model also takes the similar form of Hamiltonian. In the complex SYK model, the Hamiltonian is written in terms of the Dirac fermions c i , i = 1, · · · , N as Here A{· · · } is the antisymmetrization and the couplings J j 1 ···j q/2 ;k 1 ···k q/2 are independent complex variables with zero mean and the variance |J j 1 ···j q/2 ;k 1 ···k q/2 | 2 = J 2 (q/2)!((q/2)−1)!
The last term comes from the chemical potential µ for the generator of the global U(1) symme-tryQ = i c † i c i − N/2. When we rewrite the Dirac fermion by two Majorana fermions as c i = 1 √ 2 (ψ 2i−1 − iψ 2i ), the chemical potential term takes the same form with the mass term in the mass deformed SYK (2.1) with s k = 1 for all k [74].
The main difference is the existence of the U(1) symmetry. The complex SYK model have a soft mode that is associated to the U(1) symmetry whereas the SYK model do not have such a mode. The mass deformed SYK model has always a mass gap at zero temperature but the complex SYK model has a gapless excitation 12 . The chaos exponents are also studied in the complex model in the large q limit [74] and the µ dependence of the chaos exponent is different from the mass deformed SYK model in Fig.19.
One similarity is the specific charge Q = Q /N in the complex SYK model and the spin operator expectation value. In the complex SYK model, a natural correlation function is The specific charge is encoded in the correlation function as lim τ →0 + G(τ, 0) = − 1 2 + Q. By decomposing the Dirac fermion τ 2 ). Therefore, we can think of the specific charge Q as a counterpart of the spin operator expectation value S k = −iG off (0) in the mass deformed SYK model. A quantitative difference is that the specific charge in the complex SYK is not fixed in the IR [73], whereas the spin operator expectation value in the mass deformed SYK is determined by the IR conformal field theory data as (3.27) in small µ limit.

Possible microstates from the mass deformation
We show that we can prepare the 2 N 2 states of the form |G s (µ) from the mass deformation H def . In this paper we focus on the spin operator S k = −2iψ 2k−1 ψ 2k that is consist from an even index fermion and the odd index fermion. The way to construct the spin operator is not restricted to this form. For example, we can shuffle the index of even fermion as 2k → 2σ(k) where σ ∈ Sk 2 is the element of the permutation group Sk 2 , and then construct the spin operator S k = −2iψ 2k−1 ψ 2σ(k) . The mass deformation with S k gives a different set of states where the states have a spin operator expectation value in different directions. We can also construct with a pair of even index fermions. In this way, we can prepare many set of states as ground states of the mass deformed SYK in this paper.

Future problems
There are several future problems.
In this work we study the chaos exponent only at large q limit. It is interesting to do this at finite q numerically. We study the quantum quench problem in the small µ limit. At infinite µ, the ground state reduces to the infinite temperature boundary state |B s and in this regime real time evolutions are studied in [27,75] at finite N . It is also interesting future problem to study the real time evolution in finite µ both in large N and finite N .
In this paper we mainly study the SYK model side. Recently Jackiw Teitelboim (JT) gravity with EOW brane is studied [76]. It is a good problem to analyze the Jackiw Teitelboim gravity + matter theory with EOW brane and introduce the double trace deformation. When the brane is tensionless, JT + matter with EOW brane system just reduces to the orbifold of the traversable wormholes [5]. The analysis with the non zero tension EOW brane may lead to the bulk understanding of (the absence of) the Hawking Page like transition.
We did not find any energy/µ-dependence of the adjacent gap ratio for our model (1.1); there are no chaotic/integrable transition. This result is in contrast to the observation in [14] for the two coupled SYK model [5]. Indeed in the two coupled SYK model (7.1) the level correlation is qualitatively different in the two extreme regime µ → 0 and µ → ∞. In the limit µ → 0 the energy spectrum becomes a direct product of the energy spectrum of two SYK models {E m +E n } m,n≥0 . When the spectrum enjoys such direct product structure and there are no hierarchy between the level spacings of the two system (which is true in the current case), the two spectrums are completely mixed up. Hence there are no level repulsion between the adjacent levels even if each system has the RMT-like level correlations. In the limit of µ → ∞ the Hilbert space effectively splits into the eigenspaces of S. Within each eigenspace the direct product structure of the Hamiltonian is lost, and the levels have the RMT-like correlation. Hence one can expect the transition as µ increases. In our model (1.1), on the other hand, the picture at µ → ∞ is same as the two coupled SYK model while in the limit µ → 0 the system reduces to a single SYK model which is again chaotic.
To gain more insight on the mechanism of the Hawking-Page like transition and the chaotic/integrable transition (or their absence) and on how these two phenomena can be correlated, it would be very useful to repeat the same analysis for a generalization of the two coupled SYK model [5] such that the left coupling J L ijk and the right coupling J R ijk are chosen inde-pendently to each other. From the viewpoint of our model, this model is obtained by stating from the Hamiltonian (1.1) and then omitting all terms in H SYK which mix ψ 2i−1 's and ψ 2i 's. This model share the same features of both of the two coupled SYK model and our model. By rewriting the partition function in the large N limit by using the bi-local fields, one finds that the large N partition function is completely identical to the partition function of our model. On the other hand, the Hamiltonian of this model has the structure of direct product in the limit µ → 0 similar to the two-coupled SYK model, which strongly suppress the RMT-like level correlation in the small µ regime. It is worthwhile to test whether this model actually exhibits a chaotic/integrable transition at some finite µ or not. One can further consider an interpolation of the two coupled SYK model and this model by tuning the independentness of J L ijk and J R ijk continuously, where we observed that the Hawking-Page like transition disappears at some intermediate point before the two couplings become completely independent with each other. It would be interesting to clarify how the chaotic property as well as the other thermodynamic quantities behaves around this point. We would like to report these results in [70].
Note that it is subtle whether we should really classify a model which is almost the tensor product of two chaotic system as "integrable" although the nearest-neighbor level repulsions are highly suppressed. To clarify this point, it is worth to study other diagnoses of the quantum chaos such as the spectral rigidity or the spectral form factor (i.e. the long range correlation of the level fluctuations) and the OTOCs. Especially, while in the analysis of the level statistics one always has to take into account the finite N artifact, the OTOCs allow a direct large N evaluation [11] which would be more appropriate for the purpose of comparing the chaotic property with the large N Hawking-Page like transition.

A A derivation of the Large N equations
In this appendix, we give a derivation of the large N effective action and the Schwinger-Dyson equation of mass deformed SYK model. The deformed Hamiltonian is with mean J i 1 ···iq = 0 and variance J 2 By shifting the sign of ψ i and J i 1 ··· ,iq , we can set s k = 1 for any k = 1, · · · N/2 in the following derivation. The partition function becomes The integral over J i 1 ···iq is In the second line, the phase factor (−1) q appears from ((i) q 2 ) 2 . In the third line, we extend the sum from i 1 <···<iq to 1≤i 1 ,··· ,iq≤N . Because ψ i (τ ) is a Grassmann number, ψ i (τ ) 2 = 0 and the sum 1≤i 1 ,··· ,iq≤N survives when all of i 1 , · · · , i q are different. There are q! same contributions, we divide by q! and then the sum reduces to the sum in the second line. In the fourth line, we reorder the fermions and we get the sign (−1) q l=1 (q−l) , which becomes (−1) 2 becomes 1 because q is an even number. The partition now becomes Next, we further rewrite the partition function in terms of the correlation function G(τ, τ ) and the self energy Σ(τ, τ ). First we insert the delta functional τ >τ to (A.4): In the 2nd line, we replace the factor N i=1 ψ i (τ )ψ i (τ ) q by N q G(τ, τ ) q because we have the delta functional that relates them. Next, we represent the delta functional as the following integral 13 : Using this expression for the delta functional, we obtain 13 Strictly speaking, we need to take the correct contour to make the integral convergent.
Until now we do exactly the same transformation with that of the ordinary SYK model. From now, we further introduce the additional delta functional With this delta functional, we can replace the fermions in the mass term by G off (τ, τ ): (A.10) Next, we represent the delta functional as Using this expression for the delta functional, we get The fermion path integral gives the following functional determinant: Then, we get the effective action in terms of the G, Σ variables: (A.14)

A.1 large q expansion and Liouville action
In this section we derive the Liouville action at large q limit. The original Euclidean action is We define where G 0 (τ 1 , τ 2 ) = 1 2 sgn(τ 1 − τ 2 ) and G 0off (τ 1 , τ 2 ) = i 2 sgn(μ) is the two point function of free fermion with a Hamiltonian H = iµ k ψ 2k−1 ψ 2k with µ → 0 limit. In large q limit, the free fermion They satisfy Then, we can write them as We can expand the Pfaffian as Then, the action becomes To integrate out Σ field, it is helpful to introduce Then, this satisfies Then, the effective action becomes By integrating out Φ, we obtain the effective action as There is a nontrivial Jacobian when we change the integration variable from Σ to Φ, but this is g ab independent. We also omit the other terms that are g ab independent. Because G 0ab are constants except for τ 1 = τ 2 and g(τ, τ ) = 0, the effective action now becomes

B Numerical Solution to the Schwinger-Dyson equations
After the Wick rotation and the compactification of τ direction τ ∼ τ + β with ψ i (τ + β) = −ψ i (τ ), we can rewrite the partition function of the mass deformed SYK model as as explained in the appendix A. In the limit of N → ∞, we can evaluate the integrations over the bi-local fields by the saddle point approximation with G, Σ, G off , Σ off satisfying the saddle point equations (namely, the equation of motion for S eff ) In this section we explain how to solve the saddle point equations and evaluate S eff over the solution numerically. First we assume that the solution depends on τ, τ only through the difference τ − τ and satisfies the anti-periodicity G(τ − τ + β) = −G(τ − τ ) reflecting the anti-periodicity of ψ i , so that we can expand G, Σ, G off , Σ off in discrete fourier series G(τ ) = 1 β ∞ n=−∞ e −iωnτ G(ω n ) with ω n = 2π β (n + 1 2 ). Now, for the numerical computation, let us further discretize τ coordinate as which is equivalent to introducing an UV cutoff to ω n as −Λ ≤ n ≤ Λ − 1, so that each of G, Σ, G off , Σ off is a finite (2Λ) dimensional vector both in ω-space and τ -space and related as  and the saddle point equations are given by the ordinal derivatives of these terms by either the τ -components or the ω-components of G, Σ, G off , Σ off . It is convenient to perform G-derivative by τ -components and the derivative in Σ, G off , Σ off and we obtain (B.9) 14 Here we have renormalized the functional determinant S (1) eff such that the partition function correctly reproduces the partition function of N free fermions in the limit of J, µ → 0, as explained also in [11].
From the fourth equation in (B.9) we find that G off,n satisfies the following symmetry property: (B.10) For G n and Σ n , we find that we can consistently impose the following symmetry property though they may not be satisfied for general solutions. In the same way we can also impose the following reality relations: If we impose these symmetry properties, the Schwinger-Dyson equations (B.9) finally simplifies into the following pair of equations 14) The equation for G and Σ (B.13) can be solved numerically by using the same iteration technique as exploited in the undeformed SYK model [11, appendix G]. Once we obtain a set of solution (G, Σ) the large N partition function, or the large N free energy F = − 1 β log Z, can be evaluated from (B.4) with (B.3) as where we have also used the symmetry property of G (B.11).
Note that in the above formulation G off is merely an auxiliary field which does not contribute to the partition function, but just play a role to fix Σ off as Σ off,n = −iµ. Nevertheless G off itself is a physically meaningful observable G off (τ ) = ψ 2i−1 (τ )ψ 2i (0) and useful for the consistency check of the different approaches of the computations.

C Detail of the large q finite temperature analysis
In this appendix we give a detail of the large q analysis. As it was done in [5], we divide the range of the inverse temperature into four region that consists from the inverse temperature of order q log q, q , √ q and 1.
C.1 Inverse temperature of order β = q log q In this order, we fix σ = qe −βµ . (C.1) The general solution at early time is given by with the boundary conditions and at τ → ∞ we impose that (C.2) matches with the early time expansion (5.28) The condition at τ = 0 (C.3) gives The condition at τ → ∞ gives which lead to α =α,γ = γ + σ. (C.8) The parameter A is also determined as We ignored subleading terms in large q expansion 15 .
Here we ignore the term that contains Σ because Σ is of order q which can be ignored at the leading of the 1/q expansion. We also approximate G(τ ) ∼ 1 2 , which is the leading of the 1/q expansion that we mentioned above. Then we can solve this equation with the condition The equation for G(τ ) becomes Using the expansion G(τ ) = 1 2 (1 + g(τ ) q ) and G off (τ ) = i 2 µ( β 2 − τ ), we obtain the Liouville like equation k is now finite in our parameter regime β ∼ √ q, µ ∼ 1/q. The boundary condition forĝ is eĝ (± 1 2 ) = (βJ ) 2 , (C. 38) which is ∞ in the leading of 1/q expansion. Therefore, we should seek the solution which diverge at x = ± where we definedk = ke −ĝ 0 .

C.2.2 large k limit
In the large k limit, we can solve the equation (C.41) aŝ

Energy in SYK Hamiltonian
The plot of the half of the absolute value of the spin operator expectation value | S k |, which is equal to the τ = 0 off diagonal correlation function −iG off (0), as a function of µ. Right: The plot of the energy in the SYK Hamiltonian H SY K as a function of µ.
We have derived analytic formula for the conformal limit of spin operator expectation value (3.27), the ground state energy (3.29) and the energy in the SYK Hamiltonian (3.30). These expression contains Γ(1 − 4∆) = Γ(1 − 4/q), which diverges at q = 4. This does not mean that these observables are divergent but we should take into account the UV effects in the actual deformed SYK model. Numerically we observe that the exact value is finite but the scaling