Notes on the complex Sachdev-Ye-Kitaev model

We describe numerous properties of the Sachdev-Ye-Kitaev model for complex fermions with $N\gg 1$ flavors and a global U(1) charge. We provide a general definition of the charge in the $(G,\Sigma)$ formalism, and compute its universal relation to the infrared asymmetry of the Green function. The same relation is obtained by a renormalization theory. The conserved charge contributes a compact scalar field to the effective action, from which we derive the many-body density of states and extract the charge compressibility. We compute the latter via three distinct numerical methods and obtain consistent results. Finally, we present a two dimensional bulk picture with free Dirac fermions for the zero temperature entropy.


Introduction
The Sachdev-Ye-Kitaev (SYK) models [1,2] of fermions with random interactions have been the focus of much recent attention in both the quantum gravity and the condensed matter literature. The majority of this work has focused on the model with Majorana fermions, which has no globally conserved charge, other than the Hamiltonian itself. In this paper, we direct our attention to the model with N 1 complex fermions [3], a.k.a. the complex SYK model: ..<j q/2 , k 1 <...<k q/2 J j 1 ...j q/2 ,k 1 ...k q/2 A ψ † j 1 . . . ψ † j q/2 ψ k 1 . . . ψ k q/2 (1.1) where A{· · · } denotes the antisymmetrized product of operators. The couplings J j 1 ...j q/2 ,k 1 ...k q/2 are independent random complex variables with zero mean and the following variance: J j 1 ...j q/2 ,k 1 ...k q/2 2 = J 2 (q/2)! (q/2 − 1)! N q−1 . (1.2) One advantage of the antisymmetrized Hamiltonian is that it makes the particle-hole symmetry manifest. For example at q = 4, the Hamiltonian has the following form where C j 1 j 2 ,k 1 k 2 collects various terms arising from anti-commuting fermion operators; more explicitly, (1.4) Without C j 1 j 2 ,k 1 k 2 term, the Hamiltonian is not invariant under the particle-hole symmetry ψ † j ↔ ψ j . Using the same notation, we define the globally conserved U(1) charge Q by It is related to the ultraviolet (UV) asymmetry of the Green function (1. 6) where β = T −1 is the inverse temperature, ∆ = 1/q is the scaling dimension of the fermion operator, E ∈ (−∞, +∞) and θ ∈ (−∆π, ∆π) are the spectral asymmetry parameters related by the following formula e 2πE = sin(π∆ + θ) sin(π∆ − θ) . (1.9) Note the value of (E, θ) can not be fixed by the IR equations; so there is a one-parameter family of solutions in the scaling limit [1]. Ultimately, the actual value of (E, θ) is set by the value of specific charge Q. Although charge is a UV property of the system, the relationship between (E, θ) and Q is universal and independent of UV details [4,3]: (1.10) We will provide new derivations of Eq. (1.10) here (see Eqs. (2.34), (3.35) and Section 5). This universal relation is analogous to the Luttinger relation of Fermi liquid theory, which relates the size of the Fermi surface (an IR quantity) to the total charge (a UV quantity). The form of Eq. (1.7) also applies to fermionic fields with unit U(1) charge in asymptotically AdS 2 black holes, as was computed by Faulkner et al. [5]; the parameter E is then a dimensionless measure of the electric field near the AdS 2 horizon [3] (see Appendix G and Eq. (G.9)). For both SYK models and black holes, fields with U(1) charge p have the asymmetry factor e ±πpE .
Another key feature of the SYK models is the presence of a non-zero entropy in the zero temperature limit [4]: The function S(Q) is universal, in that it is determined only by the structure of the low energy conformal theory, and is independent of the UV perturbations to the Hamiltonian which are irrelevant to the low energy. Such a zero temperature entropy is not associated with an exponentially large ground state degeneracy. Instead, it signals an exponentially small manybody energy level spacing down to the ground state; see Section 2.5. For each given N , the entropy does go to zero at exponentially low temperatures. We will present a new derivation of S(Q) in Section 5 using a two dimensional bulk picture involving massive Dirac fermions on the hyperbolic plane. At finite but sufficiently low temperatures, the dynamics of the Majorana SYK model is governed by a collective mode with the Schwarzian action [2,6,7]. An analogous effective theory of the complex SYK model also includes a U(1) mode [8] where ϕ(τ ) is a monotonic time reparameterization obeying ϕ(τ + β) = ϕ(τ ) + 2π, and λ(τ ) is a phase field obeying λ(τ + β) = λ(τ ) + 2πn with integer winding number n conjugate to the total charge Q. Notation Sch(f (x), x) stands for the Schwarzian derivative (1. 13) In this effective theory, the U(1) and SL(2, R) freedom are actually decoupled, which can be demonstrated by the variable change λ(τ ) = λ(τ ) + iE 2π β τ − ϕ(τ ) . The action (1.12) is characterized by two parameters, K and γ, and these can be specified by their connection to thermodynamics. They depend upon the specific charge Q (or the chemical potential µ), but this dependence has been left implicit. The leading low temperature correction to the entropy in Eq. (1.11) at fixed N and Q is 14) and so γ is the T -linear coefficient of the specific heat at fixed charge, as in Fermi liquid theory. The parameter K is the zero temperature compressibility (1.15) Unlike S and E, the parameters K and γ are not universal, and depend upon details of the microscopic Hamiltonian and not just the low energy conformal field theory. The zero temperature entropy in Eq. (1.11) and the pair of soft modes as in Eq. (1.12) are also pertinent to higher dimensional charged black holes with AdS 2 horizons, and this is discussed elsewhere [9, 3, 10-12, 8, 13-18]. Key aspects of such black holes are summarized in Appendix G. We also note that supersymmetric higher dimensional black holes with AdS 2 horizons obtained from string theory have integer values for e N S [19,20], as does the SYK model with N = 2 supersymmetry [21] (which we do not consider here).
An important property of both complex SYK and charged black holes with AdS 2 horizons is the following relationship between the entropy S(Q) in Eq. (1.11) and the parameter E: (1. 16) This relationship first appeared in the study of SYK-like models by Georges et al. [4], building upon large N studies of the multichannel Kondo problem [22]. Independently, this relationship appeared as a general property of black holes with AdS 2 horizons in the work of Sen [23,24], where E is identified with the electric field on the horizon [5], as noted above. It was only later that the identity of this relationship between the SYK and black hole models was recognized [3].
We will obtain a deeper understanding of Eq. (1.16) in the present paper, based on the global U(1) symmetry associated with the conserved charge and the locality of effective action. Let us summarize our notation for thermodynamic quantities. These quantities are of the order of N : the total charge Q (which is integer for N even, and half-integer for N odd), action I, entropy S, and the associated free energy and grand potential. N -independent quantities include: the temperature T = β −1 , chemical potential µ, spectral asymmetry parameter E, specific charge Q, zero-temperature entropy S, charge compressibility K, and the T -linear coefficient in the specific heat γ. Except the first two, they are defined in the large N limit.

Outline of the paper
We begin Section 2 by setting up the formalism for the complex SYK model as a path integral over the two-time Green function and self energy. We introduce a definition of the conserved charge Q suitable for this formalism and then derive the known universal relation between Q and E (Eq. (2.34)). In Section 2.3, we find a general form of a local effective action I eff and derive Eq. (1.12). Section 2.4 is concerned with thermodynamic quantities and a discussion of what parameters come from the UV. In Section 2.5, we evaluate the path integral over λ and ϕ with action I eff exactly, which yields new results for the many-body density of states.
Section 3 sets up a renormalization theory of the complex SYK model. This will enable us to obtain another derivation of the relationship between the specific charge Q and the spectral asymmetry E.
In Section 4, we turn to the calculation of the parameters of the effective action, in particular charge compressibility K. We present three numerical computations that yield values of K in good agreement with each other. These computations and our analysis show that all energy scales contribute to the charge compressibility. A low energy analysis based on linear coupling, mentioned in Section 2.4, or conformal perturbation theory (see Appendix C) does not yield the correct value of K, even though such methods work [6,7] for the Schwarzian mode.
Section 5 presents a two dimensional bulk derivation of the zero temperature entropy of the complex SYK model. We show that the E-dependent value of the zero temperature entropy per fermion S can be obtained from a Euclidean path integral over massive Dirac fermions on hyperbolic plane H 2 . We show that the appropriate quantity is the ratio of fermionic determinants with different boundary conditions on the boundary of H 2 . Another bulk interpretation of the entropy appears in Appendix G, where we recall the connection to higher-dimensional black holes. In d + 2 dimensions (d 2), the AdS 2 arises as a factor in the near-horizon geometry of a near-extremal charged black hole. In this picture, S is related to the horizon area in the d extra dimensions, and, as we noted above, this S also obeys the differential relation (1.16).

Low temperature properties
In this section, we analyze the complex SYK model based on the (G, Σ) action. We provide a general definition of charge in this framework and prove its universal relation to the IR asymmetry of the Green function. Furthermore, we find the general form of effective action and evaluate the path integral over the low energy fluctuations, which yield new results for the many-body density of states.

Preliminaries
We start with a review of the basics. For convenience, we measure time in units of J −1 , which is equivalent to setting J to 1. For the Hamiltonian (1.1), we may consider either the partition function for a fixed charge or the grand partition function. The latter can be obtained from the (G, Σ) action: where σ(τ 1 , τ 2 ) = δ (τ 1 − τ 2 ) − µδ(τ 1 − τ 2 ) . (2.1) The Schwinger-Dyson equations are as follows: − (Σ + σ) Σ G = 1, Σ(τ 1 , τ 2 ) = G(τ 1 , τ 2 ) q 2 (−G(τ 2 , τ 1 )) q 2 −1 . (2. 2) The general idea of solving these equations in the IR limit is to ignore σ, which is localized at short times. However, care should be taken because the Fourier transform of σ contains the non-negligible, ω-independent term −µ. Fortunately, this term is absent from Σ, so we will use G and Σ as independent variables. Thus, σ moves to the second equation in (2.2), where it can be safely ignored as the equation is solved in the time representation.
Since the IR equations do not depend on µ, we get a family of solutions parametrized by a formally independent variable E. At zero temperature, · sin(2π∆) 2 cos π(∆ + iE) cos π(∆ − iE) .
We can also introduce a parameter θ to characterize the spectral asymmetry in the frequency domain: The spectral asymmetry parameters E and θ are related by the equations e −2iθ = cos(π(∆ + iE)) cos(π(∆ − iE)) , e 2πE = sin(π∆ + θ) sin(π∆ − θ) . (2.5) Using these relations, we can also express the prefactor b as a function of θ: b = 1 − 2∆ 2π · 2 sin(π∆ + θ) sin(π∆ − θ) sin(2π∆) . (2.6) The zero-temperature solutions can be extended to finite temperature: The subscript c here means "conformal". In the frequency domain with Matsubara frequencies ω n = 2π β n + 1 2 1, the Green function and self energy have the following form, (2.8) Given these exact solutions to the IR equations, it remains to be checked whether they extrapolate at higher energies to a solution to the full UV equations which depend upon µ. This has been examined numerically [1,25,4,26,27], and a consistent extrapolation exists for |µ| 0.24 [27,28]. The IR parameter E can be determined as a smooth, odd function of the UV parameter µ over this regime.
In addition to the emergent reparameterization symmetry that is present in the low energy limit of the Majorana SYK model, the complex SYK model has an extra emergent symmetry related to phase fluctuation: , where ϕ(τ ) is a monotonic time reparameterization with winding number 1 and λ(τ ) is a phase fluctuation with possibly arbitrary integer winding number. The symmetries are not exact in the presence of σ term in the (G, Σ) action (2.1). To make this point more transparent, it is useful to rewrite the action in terms of (G, Σ): (2.10) Now the first line of the r.h.s. of Eq. (2.10) is invariant under the symmetry transformation (2.9), while the second line changes. This point will be further discussed in Section 2.3.

Charge
For an explicit UV source field σ (cf. Eq. (2.1)) that arises from a microscopic Hamiltonian, the charge is conventionally defined by the UV asymmetry of the Green function as Eq. (1.6). In this section we will derive a formula for charge in (G, Σ) action framework for general source field σ (without assuming time translation symmetry) using ideas similar to those in Appendix C of Ref. [29].

"Flow" of Green function G
Let us consider the action (2.10) with σ(τ 1 , τ 2 ) depending on both times, not just on τ = τ 1 −τ 2 . If (G, Σ) is stationary (i.e. satisfies the Schwinger-Dyson equations) and β = ∞, then This can be established by considering an infinitesimal variation δλ(τ ) and the corresponding variations Only the σG term in (2.10) has a non-trivial variation, which is proportional to the l.h.s. of (2.11) if δλ(τ ) ∝ δ(τ − τ 0 ). On the other hand, the variation of the action must be zero since (G, Σ) is stationary. Following the ideas in Appendix C of Ref. [29], we may call the "current" flowing from τ 1 to τ 2 . To make a closer analogy to the aforementioned reference, let us substitute the expression σ(τ 1 , q 2 −1 (obtained from the Schwinger-Dyson equations) into the current formula: (2.14) Treating G and Σ as matrices indexed by (τ 1 , τ 2 ), we have Σ = −G −1 . If G were a unitary quasidiagonal matrix, the results in Appendix C of Ref. [29] would apply, and certain quantities would be quantized. However, here Green function G has non-trivial IR asymptotics violating the conditions of being quasidiagonal. Nevertheless, we will use similar ideas and definitions as the aforementioned reference. Note that Eq. (2.11) can be interpreted as the conservation of the current at each point τ 0 : as illustrated in Fig. 1 (a). It follows that the total current through a cross section τ 0 , (see Fig. 1 (b)) is independent of τ 0 . As explained below, this quantity is a natural generalization of the specific charge Q/N to general sources. We may call Q the "flow" of the matrix G as it depends solely on G through Eq. (2.14) with Σ = −G −1 . We also remark that the definition of the flow does not rely on the time translation symmetry. That is, the source σ(τ 1 , τ 2 ), and the Green function G(τ 1 , τ 2 ), may depend on both τ 1 and τ 2 rather than just τ = τ 1 − τ 2 . We now explain the interpretation of the flow Q as charge. Plugging the definition (2.13) of the current j into Eq. (2.16), we get This formula reduces to a simpler form when the source has the time translation symmetry, i.e. for σ(τ 1 , τ 2 ) = σ(τ ), where τ = τ 1 − τ 2 : The last expression in turn reduces to the conventional definition of the charge when σ(τ ) = δ (τ ) − µδ(τ ). In this case, for the Green function G(τ ) that is discontinuous at τ = 0, we use the average 1 2 (G(0 + ) + G(0 − )) to define its value at τ = 0. Thus, in agreement with Eq. (1.6). For extremely local UV sources such as δ (τ ) and δ(τ ), the charge is a local quantity. However, if we consider a general source σ(τ ), the r.h.s. of Eq. (2.18) includes contributions from all scales; see Fig. 1 (b) for a cartoon.

Invariance of the charge
We will show that the charge Q depends only on the UV and IR asymptotics of G(τ 1 , τ 2 ) and Σ(τ 1 , τ 2 ) (where Σ = −G −1 ) as well as some topological data. The UV asymptotics is determined by the δ (τ 1 − τ 2 ) term in Σ. To formulate the invariance, let (G 1 , Σ 1 ) and (G 2 , Σ 2 ) have the same asymptotics and in addition, let the following "relative winding number" be zero: If ν(G 1 , G 2 ) = 0, then (G 1 , Σ 1 ) can be continuously deformed into (G 2 , Σ 2 ). Here it is important to consider the winding number in frequency domain rather than time domain, because the Schwinger-Dyson equation Σ(ω) = −G(ω) −1 guarantees that a smooth path in (G, Σ) space will disallow both singularities and zeros of G(ω). This will not work for G(τ ), since the other equation Σ(τ ) = G(τ ) q 2 (−G(−τ )) q 2 −1 does not constrain zeros of G(τ ). To prove that the charge is invariant under such deformation, it is sufficient to consider infinitesimal, asymptotically trivial deformations. Let us use the formula where f is an arbitrary function such that This formula coincides with Eq. (2.17) for the step function f (τ ) = θ(τ − τ 0 ). The integral does not depend on the details of f because of the conservation law, namely Eq. (2.11). 1 More intuitively, f may be interpreted as a linear combination of step functions θ(τ − τ 0 ) with some weights for each cross section position τ 0 . In other words, we can blur the cross section, and this will not affect the flow. We proceed by anti-symmetrizing the integrand, we also get this equation: Note that the two terms cannot be integrated separately because the corresponding integrals are not absolutely convergent (since G(τ 1 , there is a logarithmic divergence in IR). Now, consider an infinitesimal variation δG and δ Σ = Σ (δG) Σ such that δG(τ 1 , τ 2 ) and δ Σ(τ 1 , τ 2 ) decay sufficiently fast as τ 1 − τ 2 → ∞: (2.25) Substituting δ Σ = Σ (δG) Σ into the first line of the last expression and using Σ(τ 4 , τ 3 ) = − dτ 2 dτ 1 Σ(τ 4 , τ 2 )G(τ 2 , τ 1 ) Σ(τ 1 , τ 3 ) in the second line, we get Now, we can regroup and integrate the terms containing f (τ 2 ) − f (τ 3 ) and f (τ 4 ) − f (τ 1 ) separately: Each integral contains a delta function that annihilates f (τ 2 ) − f (τ 3 ) in the first case and f (τ 4 ) − f (τ 1 ) in the second case. Therefore, we conclude that δQ = 0 for asymptotically and topologically trivial deformations.

Calculation of the charge
In fact, we can calculate the charge in the complex SYK model using the IR asymptotics. (The result for q = 4 was originally derived in Ref. [4] using a different method, see Appendix A for a detailed discussion.) We will use the antisymmetrized version of Eq. (2.18) to express σ in terms of Σ as we have done before: (2.28) The two terms in the last expression almost cancel each other at τ 1, but individually, the corresponding integrals are logarithmically divergent. To proceed, let us replace G(τ ) with where η is a small positive number. This change has little effect on the integrand in (2.28), but the two terms can now be separated. The corresponding integrals are equal to each other due to the symmetry τ → −τ . Thus, It is important that the symmetric-in-time regularization (2.29) is not symmetric in frequency, which has nontrivial consequences. The Fourier transform of G η is where ∆ = ∆ + η. Expanding the ω-independent coefficients in this expression to the first order in η, we explicitly see the asymmetry: where ψ(x) = Γ (x)/Γ(x) is the digamma function. Now, we are in a position to evaluate the integral over dω in (2.30), which should be understood as a principal value integral. For frequencies above some threshold ω 0 (such that ω 0 1 but η ln 1 ω 0 1), the difference between G and G η may be neglected. On the other hand, if |ω| < ω 0 , then G(iω) ∼ |ω| −1+2∆ , and hence, The last integral is simply 1/(2η), so the result is independent of η. Including the factor 1/(2πi) from (2.30) we obtain: which reproduces the result in Ref. [8].

Analogy with field-theoretic anomalies
In some sense, the calculation of the charge performed here (also see Appendix A for parallel discussions based on symmetric-in-frequency regularizations) has a flavor of perturbative anomalies in quantum field theory. Both describe a mismatch between the IR and UV. In the case of anomaly, there is no consistent UV cutoff respecting the symmetry; in our case, the UV behavior is well-defined but quantifiably different from the IR. By analogy with the Fermi liquid theory, one might expect the charge to be given by the first term in Eq. (2.34). However, that is not correct due to the non-trivial effect of regularization, which produces the second term, In Appendix A, we will further comment on this term and relate it to the Luttinger-Ward term in the standard analysis [30].

Covariant formalism and the effective action
At the low temperatures, β 1, the action (2.10) is almost invariant under the emergent symmetry (2.9). In other words, we can generate "quasi-solutions" of the Schwinger-Dyson equations by applying (ϕ, λ) transformations to the actual solution (G * , Σ * ): (2.36) Such quasi-solutions cost a small increase in action, 2 which we call the "effective action". More exactly, I eff [ϕ, λ] = I[G, Σ] − const, where the constant depends on convention: we may set it to I[G * , Σ * ] or use a different base value. The goal of this section is to determine the general form of I eff [ϕ, λ] to leading orders.

Covariant formalism
The form of the approximate solutions coincides with the transformation laws for functions changing from one frame to another. The latter is described by a diffeomorphism of the time circle together with a gauge transformation. For instance, a field ψ(x) with scaling dimension ∆ and charge p defined in frame "x" can be transformed to frame (y, φ) by the following formula: It is also straightforward to define the transformation laws for G and Σ, i.e. functions of two variables. Taking this viewpoint, the "quasi-solution" (2.36) may be interpreted as follows. In order to generate a quasi-solution (G, Σ) in the physical frame x = τ , we start with the frame (y, φ) = (ϕ, λ), where the Green function is given by G * . Then we pull back to the physical frame using the inverse transformation of (2.37), namely From this perspective, the effective action I eff [ϕ, λ] in fact measures the failure of (ϕ, λ) to be the physical frame. At first glance, introducing the notion of "frame" in this problem seems redundant because in the end we should write all results in the physical frame. However, the condition that the action is invariant under frame transformations (if we also transform σ) is helpful in determining the general form of the effective action I eff [ϕ, λ]. Now, let us consider expressing the (G, Σ) action in a general frame (ϕ, λ). (In this setting, G and Σ are arbitrary and not related to ϕ and λ in any particular way.) In the new frame, the fields are defined as follows: (2.39) Representing G, Σ, σ in terms of G ϕ,λ , Σ ϕ,λ , σ ϕ,λ and plugging into Eq. (2.10), we get (2.40) Naively, the ln det term transforms nontrivially under ϕ. However, the determinant needs some UV regularization anyway, and we will use a particular regularization to make it frameindependent: 3 The second factor is the free fermion partition function (formally equal to det(−σ)).
The expression (2.40) for the (G, Σ) action in frame (ϕ, λ) has the same general form as in the physical frame, but the UV source in different: (2.43) We will make a few comments on this transformation.
• First of all, the non-trivial change of the source lifts the degeneracy of quasi-solutions and induces the effective action I eff [ϕ, λ], which can be explained as follows. If σ ϕ,λ were given by the same expression as the standard source, σ std (ϕ 1 , ϕ 2 ) = δ (ϕ 1 − ϕ 2 ) − µδ(ϕ 1 − ϕ 2 ), then (G ϕ,λ , Σ ϕ,λ ) = (G * , Σ * ) would be a stationary point of the action (2.40) in frame (ϕ, λ), and therefore, its pullback (G, Σ) would satisfy the Schwinger-Dyson equations in the physical frame. Thus, the actual distinction between solutions and quasi-solutions can be attributed to the difference between σ ϕ,λ and σ std . In the first approximation, the effective action is obtained by integrating σ ϕ,λ − σ std against G * .

Diffeomorphism-invariant effective action
Now we are ready to determine the general form of the effective action. Let us consider the following quasi-solution: where ϕ maps the time circle to the standard circle of length 2π (i.e. ϕ(τ + β) = ϕ(τ ) + 2π) and G * is the IR solution of the Schwinger-Dyson equations with β formally set to 2π: To separate the U(1) and SL(2, R) degrees of freedom, let λ(τ ) = λ(τ ) + iE 2π β τ − ϕ(τ ) so that In general, the effective action contains local and non-local terms. 4 The local part describes interactions between the UV sources ε, µ and some IR data. We have discussed the origin of the fields ε, µ in last section; now let us search for the IR fields by the intermediate asymptotic where τ + = τ 1 +τ 2 2 , and the coefficients A, B are obtained by Taylor expanding the quasi-solution (2.49) w.r.t. small time separation τ 1 − τ 2 . These coefficient have the following explicit form: Thus, all relevant IR information is contained in the fields The local part of the action should have a covariant expression in an arbitrary frame x. We aim to find the effective action to β −1 order. In other words, the action should have the accuracy to recover the free energy or grand potential to T 2 order and capture, for example, the linear dependence of the specific heat. With the UV source (ε, µ) and IR data (A, O), the most general diffeomorphism-and gauge-invariant action is Let us discuss some of its features as well as defining the field ρ.
• In addition to the fluctuating fields ϕ and λ, the action depends on the global parameter E.
Its equilibrium value will be determined by finding an extremum (actually, a maximum) of I eff [E, ϕ, λ] in E with fixed external parameter β, µ.
• The function f is a general even function characterizing the charge response to the chemical potential. The gauge invariance (2.43) requires that any dependence on µ be through the combination µ − A. The G(E) term is of order 1 and related to the zero temperature entropy.
• We have expressed the action in a general frame x to emphasize its covariant properties.
In particular, we have used the notation O x and introduced a field ρ x (see [7], Eq. (167)): This can be further summarized in a matrix form. In fact, the pair (1, O) forms a representation of Diff(S 1 ), Similarly, the pair (ε, ρ) also forms a representation, Thus, the following combination transform as a 1-form which further implies the diffeomorphism invariance of the action (2.53).
The form of the ρ field may look obscure; however, it will naturally arise when we try to transform the Schwarzian action from the physical frame τ to a general frame x(τ ) using the chain rule and express everything in the new frame In other words, in the physical frame ε = 1, ρ = 0, and the last term in the effective action (2.53) is just the familiar Schwarzian action −α S Sch(x, τ )dτ . Now, let us restrict to the physical frame x = τ . Expanding f (µ − A) to the second order in A = i λ − 2π β E, we get where n is the winding number of the function λ, defined by the generalized periodicity condition λ(τ + β) = λ(τ ) + 2πn. The second line in the above equation is equivalent to the action (1.12), originally derived in Ref. [8], with To further simplify the effective action, let where λ has zero winding number. Then

Thermodynamics
We now use the effective action I eff [E, n, ϕ, λ] to compute the low temperature expansion of the grand potential Ω(β −1 , µ). If N is large, we may use the saddle point field configurations, ϕ(τ ) = 2π β τ , λ(τ ) = 0, n = 0, and find the extremum I * of I eff [E, ϕ, λ] with respect to E: The saddle point condition for E requires that This function has been calculated in Eq. (2.34), so one can compute G(E) as well.

Free energy and entropy
We can also find the free energy F (β −1 , Q) = Ω + µQ: where S = S(Q) is the "zero temperature entropy" per site and The formula (2.69) says that S(Q) and G(E) are related by the Legendre transformation, where Q and 2πE are the conjugate variables. It leads to the fundamental relation (1.16) between the entropy and the particle-hole asymmetry. Equation (2.68) is the usual Legendre duality between the free energy and the grand potential at zero temperature; it implies that µ 0 = dF 0 (Q)/dQ. At finite temperature, the chemical potential receives a definite correction: Similar relations hold for the low energy limit of charged black holes [5,17], as reviewed in Appendix G. The low T limit must be taken at fixed Q with µ obeying (G.6), to obtain a near-horizon metric that is conformally equivalent to AdS 2 .

Charge compressibility
As discussed in Refs. [6,7], the specific heat is determined by the prefactor α S of the Schwarzian action, which is related to the magnitude α G of the leading UV-sourced correction to the IR Green function. Specifically, where k c (2), ∆ and b are all IR data that can be obtained in the conformal limit. Now, for the complex SYK model we have one more thermodynamic coefficient to determine, namely the charge compressibility K. A natural question is whether the charge compressibility can be determined in a similar way by the same UV parameter α G . This possibility is based on the observation that the IR degrees of freedom A(τ ), B(τ ) in Eqs. (2.50), (2.51) satisfy the relation which might constrain the form of the effective action. Or is the charge compressibility independent of α S and requires a separate numerical study?
To answer this question, let us think about possible couplings between the IR degrees of freedom and some UV data. The idea of renormalization theory, as used in Ref. [7], is not to solve the actual problem in the UV (which is hard) but to replace it with a more tractable model with sufficiently many parameters that would reproduce the leading IR behavior and any possible corrections to it. The simplest term to include in the (G, Σ) action is the linear coupling dτ 1 dτ 2 σ(τ 1 , τ 2 )G(τ 2 , τ 1 ) of the UV source to the Green function, where the latter is represented by the asymptotic expression (2.50) at intermediate time intervals τ 1 − τ 2 with coefficients A τ 1 +τ 2 2 and B τ 1 +τ 2 2 . By smearing the actual, very singular source, nonlinear effects can be reduced. In this approximation, the effective action is a sum of terms proportional to A(τ ) dτ and B(τ ) dτ . Any contribution of the form A(τ ) 2 dτ is due to the second term via Eq. (2.72). But on the other hand, see Eq. (2.53). Therefore, the linear model predicts the following value of the charge compressibility K = −f (µ): However, a nonlinear coupling of the form 5 s(τ 1 , τ 2 , τ 3 , τ 4 )G(τ 1 , τ 2 )G(τ 3 , τ 4 ) d 4 τ can also generate a term proportional to A(τ ) 2 dτ . Let us denote this additional contribution by K non−linear so that Its actual value is not accessible without numerics. In Section 4 we will present numerical computations for the total K at half filling, namely µ = 0 and Q = 0. We would like to make a final remark on the ratio K linear /γ in (2.74). It agrees with Eq. (C.45) obtained from a different analysis following the perturbation theory done in Ref. [6]. This analysis relies on the UV parameter α G , see Appendix C for details. Similarly to K linear , it does not include the additional non-universal UV contributions to the compressibility.

Partition function at low temperatures and the density of states
We first overview some relevant results for the Majorana SYK model. An interesting time scale in the problem is given by the coefficient of the Schwarzian action, α S N . If the inverse temperature β is of this order of magnitude or greater, quantum fluctuations are strong. This regime was originally studied in Ref. [31] (see also [32]). The density of states (DOS) and the partition function for the pure Schwarzian model are as follows [33][34][35][36], where the energy E and the temperature β −1 are measured in units of (α S N ) −1 : These functions are defined up to an overall factor that depends on the normalization of the integration measure. The DOS and the partition function for the Majorana SYK model contain some additional factors. Up to a common overall constant, where, E 0 is the ground state energy. The factor N −3/2 in the partition function has been introduced to obtain the correct asymptotic behavior for β 1 fixed and N going to infinity: Note the absence of a ln N term. Indeed, the logarithm of the partition function admits a 1/N expansion, where different terms correspond to different classes of Feynman diagrams. In particular, the N 0 term is given by the sum of ladders closed into a loop, yielding the expression − 1 2 Tr ln(1 − K G ). Here K G is the exact ladder kernel; it has ∼ β eigenvalues that are not too small, whereas the − 3 2 ln β contribution is due to the eigenvalues close to 1. There is one more thing to take into account-variations between different samples: In Eq. (2.77), E 0 should be understood as the ground state energy for a particular realization of disorder. This may or may not be important, depending on the parameter q. Indeed, the sample-to-sample fluctuations are dominated by a particular Feynman diagram that contributes to ln Z but not to ln Z [7]. Assuming that N β 1, its value can be estimated as follows: Therefore, the fluctuations of the free energy are of the order of δF ∼ N 1−q/2 with no singular temperature dependence. We conclude that for β ∼ N , the sample-to-sample fluctuations are significant if q = 4 but not at larger values of q.
For the complex SYK model, the density of states is a function of two conserved quantities, charge and energy: where H and Q are defined in Eqs. (1.1) and (1.5), respectively, and Π Q is the projector onto the subspace with a given value of Q. For simplicity, we assume that N is even so that Q takes on integer values. The partition function for a fixed Q and the grand partition function are as follows: In analytical calculations, we will be interested in the case where E is close to E 0 (Q), the lowest eigenvalue of H with charge Q. We will show that up to a constant factor, or equivalently, up to a constant term. However, E 0 (Q) is difficult to compute with sufficient (say, 1/N ) precision; for q = 4, it depends on the realization of disorder. A simple though not very accurate approximation is as follows: We now derive Eq. (2.84) from the effective action. Note that the integration measure is defined up to some N -dependent factor. 6 We will use the factor N −3/2 that comes with Z Sch as previously explained. The additional normalization factors in our calculation are reasonably well motivated but not trustworthy, so the overall power of N in front of Z Sch will be checked independently.
Using the effective action (2.63), the grand partition function is expressed as follows: Let us focus on Z U(1) , which involves the variables n, E, and λ. The notation Dλ/U(1) indicates that we consider λ(τ ) up to an additive constant. The corresponding integral, may be interpreted as the partition function (per unit length) of a free particle with mass N K.
The integral over E is evaluated using the method of steepest descent. Since the integration path is parallel to the imaginary axis, and the symbol DE is understood as dE/i up to some real factor of the order of 1. Thus, For each value of the winding number n, the effective action attains its extremum at the value of E determined by the equation . Replacing the right-hand side with 2πf (µ) introduces an O(β −1 ) error in E and an O(β −2 ) error in I eff ; the latter is within the precision of the effective action model. 7 The value of ∂ 2 I eff /∂E 2 at the extremum is also almost independent of n. Applying the method of steepest descent and choosing the order 1 factor for future convenience, we get where, as in Section 2.4, The sum over n in (2.90) is evaluated using the Poisson summation formula: An important feature of the second line in (2.92) is that the argument of F is the integer charge Q being summed over, and not the mean charge QN ; consequently the entropic prefactor of each term in the sum is e N S(Q/N ) and not e N S(Q) . (Since dS/dQ = 2πE, the ratio of such factors in the Q + 1 and Q sectors is e 2πE .) So, when we multiply the second line in (2.92) by N −3/2 Z Sch (α S N β −1 ), we obtain an expression for Z(β −1 , µ) that is equivalent to (2.84). Finally, this yields (2.83), where the density of states at charge Q has a prefactor, e N S(Q/N ) , with entropy evaluated at the same charge Q.
On the other hand, if N is very large, the sum over n in (2.90) is reduced to the n = 0 term. Multiplying it by the same factor, we get The absence of a ln N term is consistent with 1/N expansion.

Renormalization theory
In this section, we describe the physics at intermediate time scales, 1 τ β, generalizing the ideas in Ref. [7] section 3. More exactly, we will study the renormalization of both symmetric and anti-symmetric perturbations to the Green function and the self energy.

General idea
The (G, Σ) action (2.10) is suited for the perturbative study near conformal point (G c , Σ c ), which is an exact saddle point for σ = 0. We will work at zero temperature, i.e. G c = G β=∞ , Σ c = Σ β=∞ , see (2.3). The actual UV source σ, consisting of a delta function and its derivative, is strong in the UV (i.e. for τ := τ 1 − τ 2 ∼ 1), and therefore, is hard to study without numerics. However, it is possible to introduce a weaker perturbation in a slightly extended UV region such that its effect at τ 1 reproduces the actual correction (δG, δ Σ) to the conformal solution. This method has been applied to the Majorana SYK model in Ref. [7] section 3, yielding a derivation of the Schwarzian action as well as the relation between its coefficient and the UV-sourced correction to the Green function.
One useful property of the Majorana SYK model is anti-symmetry under time reflection. Namely, the perturbation source δ (τ ) is an anti-symmetric function of time, and the ladder kernel that propagates the perturbation preserves this symmetry. As a consequence, the responses δG(τ ), δ Σ(τ ) are also anti-symmetric in time. However, that is not the case for the complex SYK model, see Fig. 2 for an illustration. The actual UV source σ(τ 1 , τ 2 ) = δ (τ 12 ) − µδ(τ 12 ) has both anti-symmetric and symmetric parts. More importantly, the ladder kernel (which will be studied later) mixes anti-symmetric and symmetric functions. The mixing effect will be characterized by a 2 × 2 matrix that generalizes the number k c (h) of the Majorana SYK model [6,7].
In general, renormalization theory determines how UV sources manifest themselves at intermediate scales, and thus, affect the IR physics. For instance, the interaction between the UV source and the IR deformation of the conformal solution due to reparameterization of time ϕ(τ ) contributes to the local part of the effective action for the ϕ field: it generates the Schwarzian term, which further determines such properties as specific heat. For the complex SYK model, the new ingredient is the perturbation due to chemical potential (or charge Q), sourcing the asymmetry of the Green function characterized by E or θ. The nontrivial relation (2.34) between θ and Q can be reproduced using renormalization, which further supports the statement that the charge is determined by the intermediate asymptotics of G.
To apply the renormalization theory for β = ∞, we write the charge as 18)), with σ being small at each individual time scale but possibly spanning multiple scales. This way, σ is regarded as a combination of infinitesimal perturbation sources. Focusing on a particular scale ξ = ln |τ |, we may characterize the cumulative effect of all sources at smaller scales by some value of E. The additional source at scale ξ contributes both to Q (via integral (3.1)) and to E (via renormalization). Thus, one can calculate dE/dQ, as elaborated in the following sections. The change in the asymmetry parameter E is propagated by the RG flow and further augmented by any sources present at larger scales.
3.2 Linear response to the perturbation σ

Quadratic expansion of the (G, Σ) action
In this section, it will be convenient to treat bilocal functions f (τ 1 , τ 2 ) as operators (i.e. matrices indexed by τ 1 and τ 2 ) for which one can consider the product and the trace. A similar interpretation is also applicable to functions of four times. For example, With this in mind, we can express the (G, Σ) action (2.10) as follows: Here the power G q/2 is taken element-wise. Next, we expand the action I = I[G, Σ] to the second order in the variations around the conformal point, δG = G − G c , δ Σ = Σ − Σ c , ignoring the constant term I[G c , Σ c ]: The matrices W Σ and W G can also be expressed as follows: where the functional dependences of G on Σ and of Σ on G are given by the Schwinger-Dyson equations: (The equation Σ = Σ + σ is not used.) These are the explicit formulas and diagrammatic representations of those matrices: (3.9) (An arrow pointing from τ to τ denotes G c (τ, τ ).)

Ladder kernels
To calculate the effects of the perturbation source σ, we may first eliminate δ Σ from the quadratic action by evaluating it at the saddle point, δ Σ = W −1 Σ δG: We further take the saddle point with respect to δG and find its equilibrium value, which may be interpreted as the sum of ladder diagrams. The corresponding δ Σ is expressed in a similar way: (3.12) The ladder kernels K G , K Σ are products of W Σ and W G in different orders, and thus, have the same spectrum (excluding 0). Let us give their diagrammatic representations: (3.13)

Calculation of K G (h) and K Σ (h)
Due to SL(2, R) symmetry, K G and K Σ preserve power-law functions such as σ(τ ) ∼ τ 2∆−1−h . More exactly, we will consider perturbation sources of the form which generate the responses The goal is to find the 2 × 2 matrices It should be combined with the Fourier transform The relevant values of α are 2∆ − 1 + h for δG and 1 − 2∆ + h for δ Σ. Thus, The matrix W G (h) is obtained from (3.9); it is in fact independent of h: Note that this equation is related to the fact that K G (τ 1 , τ 2 ; τ 3 , τ 4 ) = K Σ (τ 4 , τ 3 ; τ 2 , τ 1 ). Therefore, K G (h), share the same eigenvalues.

Resonant response
Resonances occur at special values of h such that det(1 − K Σ (h)) = 0. In particular, h = −1, 0, 1, 2 are solutions of this equation, also see Appendix B for discussions on the other solutions. Among them, the h = 2 and h = 1 perturbation sources determine the coefficient α S in the effective action and the parameters E, Q, respectively. The dual values, 1 − h = −1, 0, correspond to IR degrees of freedom, namely, the fluctuating fields ϕ(τ ) and λ(τ ). Let h = h I be some resonance. The power-law source σ(τ ) ∼ τ 2∆−1−h I results in a divergent response, and therefore, has to be regulated. For this purpose, we multiply the source by a window function u(ln |τ |): Assuming that u has finite support, σ I vanishes in the IR so that any response at sufficiently large scales is due to RG flow. On the other hand, the window should be sufficiently wide and u(ξ) vary slowly with ξ, such that σ I (τ ) can be decomposed into power-law sources with h close to h I . Following the argument in Ref. [7] section 3.1 we conclude that at sufficiently large τ , A similar formula can be obtained for δG. Note that this result is independent of the details of window function u.

The h = 1 resonance
As already mentioned, this resonance is related to the parameter E and Q. So, let us find the residue of (K Σ (h) − 1) −1 at h = 1. First, we compute W Σ (1) and W Σ (1): Thus, (3.25) The matrix K Σ (1) has eigenvalues −(q − 1) and 1 as expected. The left and right eigenvectors associated with the eigenvalue 1 are: (3.26) By abuse of notation, we introduce the number not actually defining k(h). Thus,

Calculation of dQ/dE
As part of the renormalization scheme for Q and E, we calculate the variations of these parameters due to a perturbation source at a particular scale. More specifically, we consider the source (3.22) with h I = 1: The source also determines δ Σ through equations (3.23) and (3.28): This result may be interpreted as a change of the asymmetry parameter E. Indeed, it follows from Eq.

Computation of the compressibility
This section will present three different numerical approaches to computing the charge compressibility K of the complex SYK model. We will limit these computations to the particle-hole symmetric case, where Q = 0 and µ = 0. These computations will involve determination of the response of the particle-hole symmetric solution to small non-zero Q or µ: 1. In Section 4.1, we will compute the compressibility by an exact diagonalization, evaluating the ground state energy E 0 as a function of small Q.
The value of E 0 is determined by the UV structure of the model, and we therefore expect K to also be sensitive to the UV structure. This is as in Fermi liquid theory, where the compressibility involves a new Landau parameter, F s 0 , and is not determined by just the quasiparticle effective mass m * . In contrast, in both Fermi liquid theory and the SYK models, the T -linear coefficient of the specific heat is determined by low energy physics: in Fermi liquid theory by m * , and in the SYK model by the leading low energy deviation of the conformal solution [2,6].
2. In Section 4.2 we will numerically compute K by an alternative method: full numerical solution of the Schwinger-Dyson equations of the SYK model.
3. Finally, numerical approach in Section 4.3 employs diagonalization of the two-particle kernel.
The values of K obtained in these subsections are in excellent agreement. Throughout the whole section we will recover J.

Exact diagonalization
In this subsection we perform exact diagonalization of the complex SYK Hamiltonian for q = 4. Note that there is little dependence of K on N . We also show the N dependence of E 0 in Fig. 3(b).

Schwinger-Dyson equation
Here we briefly review the numerical solution of the Schwinger-Dyson equations for the complex SYK model. This has already been discussed in many papers, see particularly Refs. [6,8]. Our main purpose is to show that this method gives compressibility K very close to the result obtained in Section 4.1 from exact diagonalization.
We solve Schwinger-Dyson equation numerically using the well-known method of weighted iterations . For non-zero chemical potential it is convenient to start iterations with the conformal answer, regulated at the boundaries τ = 0 + and τ = β − , and with the θ parameter corresponding to specific charge Q close to expected numerical value. This prevents iterations from falling into exponentially decaying solution. We find Q numerically using the formula For large βJ we can use equation (2.34) to find parameter θ in conformal solution (2.7). We plot an example of exact numerical G(τ ) and its conformal fit G c (τ ) in Fig. 4. The grand potential can be computed from the expression [8] − from which one can obtain the entropy as  where Q is computed numerically using (4.2). Finally compressibility in units J can be obtained numerically by using the formula Numerically we fix J = 1 and compute the ratio Q/µ for small T and small µ. We first approximate the result to the zero temperature to obtain K as a function of small µ, as shown in Fig. 5, left. Then we approximate such K(µ) to µ = 0 ( Fig. 5.b, right). We did computations for q = 4 and used 10 8 grid points for the two-point function. The value of K we found is

Kernel diagonalization
This type of numerics was first done in Ref. [6] for the antisymmetric kernel 8 . In Appendix C.1 we discuss analytical approach for kernel diagonalization. (also see Ref. [8] Appendix F). The fluctuation analysis here is complementary to that in Section 3 in the sense that here we expand the fluctuations around the exact saddle while in the Section 3 we expand around the conformal saddle.
We remind that we are working on the saddle with Q = 0, where the general expressions for the kernel (3.13) have additional symmetry, i.e. they commute with the operator that switches two times, and thus we may analyze the kernel on the subspace of antisymmetric and symmetric functions separately. For this purpose, let us consider the symmetrized antisymmetric and symmetric kernels 9 where we fix β = 2π so all angles take values in the interval [0, 2π]. Since these kernels are invariant under the translation of all four times, i.e. they commute with the operator D = i(∂ θ 1 + ∂ θ 2 ), one can look for the eigenfunctions of the kernels Ψ A/S h,n (θ 1 , θ 2 ), which are simultaneously eigenfunctions of the operator D: (4.8) Let us also define variants of the kernels with the parameter n accordingly, such that φ A/S h,n (θ) are the eigenfunctions with eigenvalue k A/S (h, n), more explicitly, To numerically diagonalize kernels K A/S n (θ, θ ) in the space of antisymmetric/symmetric functions (on the discretized coordinates θ, θ ), it is more convenient to impose the symmetry explicitly, namely, we use (K A n (θ, θ ) − K A n (θ, −θ ))/2 and (K S n (θ, θ ) + K S n (θ, −θ ))/2 in the actual calculation.
We expect to find the highest eigenvalue of the kernels K A n (θ, θ ) and K S n (θ, θ ) for large βJ , (b) Plot for nα S K (n) for q = 4. One can see that within computational accuracy α S K (n) almost does not depend on n, confirming the expectation (4.11). We use 10 8 grid points for numerical computation of G(θ) and 10 5 grid points for the kernel discretization in θ and θ directions, so the kernel becomes a 10 5 × 10 5 matrix.
These eigenvalues correspond to h = 2 and h = 1 modes. The Schwarzian coupling α S and compressibility K is related to α A K and α S K through the formulas 10 where α 0 = 2πq cot(π/q)/((q − 1)(q − 2)). We compute numerically k S for q = 4 and different values of βJ and n. The plot of k S for q = 4 and n = 1 is represented in Fig. 6(a). By fitting the data points by polynomial in 1/βJ we obtain From this fit we find that α S K = 9.2 and from (4.12) we obtain for q = 4 K ≈ 1.04/J . (4.14) This agrees very well with K obtained from the Schwinger-Dyson equation (4.6) and exact diagonalization. We also plotted nα S K (n) in Fig. 6(b), where α S K (n) obtained from fitting k S for different n and α S K (1) = 9.2. One can see that within computational accuracy α S K does not depend on n in agreement with expectation (4.11).
Following the discussions in Ref. [6], one might expect that the numerical result of α S K can be related to the deviation of the exact Green function from the conformal one similar to the case for α A K (cf. Ref. [6] Eq. (3.88)). We present a calculation following this procedure in Appendix C. The result does not agree with the numerical value of K but agrees with the K linear (see (2.74)) On the other hand, the numerical result for α S from the anti-symmetric kernel agrees perfectly with the theoretical computation [6]. The reason for the disagreement for K presumably related to the fact that K is a UV sensitive quantity, and the naive perturbation theory for the symmetric kernel in 1/βJ series does not work well, e.g. the integrals obtained from higher corrections to the Green function have uncompensated power-law divergences which then contribute to the first order 1/βJ term, changing the final result. One sign of such a breakdown of perturbation theory is visible in our numerical results for eigenvectors of the symmetric kernel. They agree with the conformal kernel eigenfunctions everywhere except UV region, whereas for the antisymmetric eigenfunctions the agreement is perfect everywhere; see .
The divergence of the eigenfunctions of the symmetric kernel in UV region is captured in the large q limit (see Appendix D).

Bulk picture and zero-temperature entropy
In this section, we find the zero-temperature entropy S of the complex SYK model by considering a massive Dirac fermion in AdS 2 . The actual calculation is done in the Euclidean case, that is, on the hyperbolic plane. The asymmetry of the Green function (2.7) may be interpreted as a phase factor with an imaginary phase, 2πiE, suggestive of an imaginary U(1) field acting on the Dirac fermion. (It corresponds to a real electric field in AdS 2 .) The partition function in the presence of such a field yields the dependence of S on E, and hence, on Q via Eq. (3.34). We will find that the S so obtained is exactly equal to that obtained from direct computations for the complex SYK model [4,3,8].
Our computation of S should be contrasted with that for higher-dimensional charged black holes [37,5,3,[13][14][15][16][17][18]38], summarized in Appendix G. In the latter case, the value of S in Eq. (G.4) is determined by the horizon area and has no direct connection to the parameters of the SYK model. The present section interprets S as the contribution of fermionic fields; such matter fields [5] only make a subdominant contribution to thermodynamics in the conventional higher dimensional AdS/CFT correspondence.

General idea
For illustrative purposes, we will use the Majorana SYK model, Among many methods of calculating its zero-temperature entropy S = S Majorana , the formula can be derived by evaluating 1 2 ln det(− Σ) with proper regularization [39,2] (see also appendix E). Indeed, S is defined as the zeroth order term in the 1/β expansion ln Z N = − E 0 N β + S + O(β −1 ), where ln Z may be approximated by minus the (G, Σ) action at the saddle point. As explained in appendix E, the double integral part of the action has β and O(β −1 ) terms but no constant term.
For the complex SYK model, Z should be understood as the grand partition function, and S should be replaced by its Legendre transform, G(E) = S(Q) − 2πEQ. We will derive a formula similar to (5.2) by considering ln Z in the β → ∞ limit and extracting the constant term: 11 For Σ asymmetric in time and frequency, the direct calculation of det(− Σ) is fraught with regularization difficulties. This is where the bulk picture offers a crucial advantage, replacing the tricky UV regularization with a simple subtraction of a boundary contribution.
In an abstract sense, the bulk is an artificial system that mimics most important properties of the real one. It may also be regarded as a heat bath for a small subset of sites [7]. The following argument seems to apply to all large N systems, but we will focus on the Majorana SYK model for simplicity. Consider adding an extra site to N existing ones and modifying the couplings J j 1 ,...,jq accordingly, multiplying them by N N +1 q−1 2 ≈ 1 − q−1 2N . In the thermodynamic limit, the logarithm of the partition function is proportional to N , and its change by the stated procedure is just ln Z N . Calling the original N sites a "bath", we get: where "full" refers to the bath and the extra site together, but with the couplings unchanged. In the β → ∞ limit, ∂ ln Z N ∂ ln J = − E 0 N β + O(β −1 ); hence, the last term in the above equation may be neglected.
To calculate ln Z full − ln Z bath , we may write the Hamiltonian as H full = H bath + i χ ξ, where χ represents the extra site and ξ is a certain operator acting on the other sites. When N is large, ξ is Gaussian, meaning that the bath is completely characterized by the correlation function T ξ(τ 1 ) ξ(τ 2 ) = −Σ(τ 1 , τ 2 ) while higher correlators are obtained by Wick's theorem. This suggests the replacement of the real system by a collection of Grassmann variables Ψ j with a quadratic action I = − 1 2 j,k B jk Ψ j Ψ k , where the indices take values on the time circle (for the extra site) and some abstract locations (for the bath). The full matrix B has this structure: with σ and B bath being square and skew-symmetric, and Y rectangular. Using this artificial model, we get where we have used the identity det ( A B C D ) = det D · det(A − BD −1 C). While the previous description leaves many possibilities for choosing B bath , the nicest one is a Majorana fermion with mass M = 1 2 −∆ on the hyperbolic plane. All its properties follow from those of the Dirac fermion, studied in the next subsection and appendix F. In this preliminary discussion, we use the Poincare half-plane model with the metric ds 2 = (dτ 2 + dy 2 )/y 2 (y > 0). A Majorana spinor ψ has two components, ψ ↓ and ψ ↑ . Solutions of the equation of motion have this asymptotic form: The boundary condition ψ − (τ ) = 0 is chosen, which prescribes a sufficiently fast decay near the boundary. We will refer to it as the "Dirichlet b.c." and to the condition ψ + (τ ) = 0 as the "Neumann b.c.". Assuming that only the first term in (5.7) is present, we can promote the asymptotic coefficient ψ + (τ ) to a field and identify it with the field ξ(τ ) characterizing the bath. This is reasonable because the correlator ("+" for Dirichlet, "−" for Neumann) (5.8) matches ξ(τ 1 )ξ(τ 2 ) = −Σ(τ 1 , τ 2 ) if the "+" sign is chosen. The part of the action involving the boundary fermion χ(τ ) is Since we are interested in low temperature properties, or large time scales, the χ∂ τ χ term may be neglected. Thus, χ becomes a Lagrange multiplier field, forcing ψ + to vanish. This indicates a change from the Dirichlet to Neumann boundary condition. The corresponding asymptotic coefficient ψ − (τ ) may be identified with χ(τ ), whose correlator is −G(τ 1 , τ 2 ). To summarize, the zero-temperature entropy of the Majorana SYK model is where [· · · ] reg denotes the constant term in the 1/β expansion. The partition functions Z bath and Z full correspond to a Majorana fermion on the hyperbolic plane with the Dirichlet and Neumann boundary conditions, respectively. For the complex SYK model, one should consider G(E) instead of S and use a Dirac fermion. The calculation will follow. We note that this procedure is similar to that used to compute the influence of double trace operators on the free energy in the AdS/CFT correspondence [40][41][42][43][44][45].

Dirac fermion on the hyperbolic plane
Now we describe a realization of the auxiliary "bath" system for the complex SYK model. The abstract action I bath = −Ψ † B bath Ψ is chosen in the form Specific to two dimensions, the spin connection factors into a scalar and a constant matrix: (Further details, such as the expressions for the Dirac matrices γ 1 , γ 2 and the spin matrices Σ ab = 1 4 γ a , γ b , can be found in appendix F.) The Majorana case differs in that ψ ↓ , ψ ↑ are real, the U(1) gauge field A is absent, and the action has an overall factor 1 2 . We use the Poincare disk model of the hyperbolic plane H 2 : 14) The U(1) gauge field A is imaginary (but becomes real upon the analytic continuation from the hyperbolic plane to the global anti-de Sitter space sharing a diameter of the Poincare disk). More specifically, Thus, the model is characterized by the Dirac mass M and the field strength E. We also need to specify a boundary condition. To this end, we note that a general solution of the Dirac equation (γ c ∇ c + M )ψ = 0 has this asymptotic form near the boundary: where The dependence on the polar angle ϕ in Eq. (5.17) is a consequence of gauge choice: we use the local frame (vielbein) shown in Fig. 8, whose orientation relative to the tangent vector ∂ ϕ depends on ϕ. For the bath model, we postulate the Dirichlet boundary condition, ψ − (ϕ) = 0. But when the bulk fermion is coupled to a boundary fermion, the correct condition is Neumann, ψ + (ϕ) = 0. The Euclidean propagator for each boundary condition, with the matrix structure is calculated in appendix F, see Eq. (F.47). In particular, when both x 1 = (r 1 , ϕ 1 ) and x 0 = (r 0 , ϕ 0 ) approach the boundary, the propagator becomes where for 0 < ϕ 1 − ϕ 0 < 2π, we have

Subtraction of infinities and the "spooky propagator"
We are now in a position to evaluate the thermodynamic quantity Each of the two terms in the square brackets suffers from a UV divergence and the divergence due to infinite volume. The former is canceled due to the subtraction of the terms and the latter due to the regularization [· · · ] reg , which amounts to the subtraction of a boundary contribution. The two terms exactly cancel each other if M = |E|. For M > |E|, it is convenient to take the derivative with respect to M using the relation (5.22) between M and ∆: In the last expression, −C + may be regarded as a propagator of a ghost particle. For this reason, we call the difference C sp = C − −C + the "spooky propagator". The function C sp (x 1 , x 0 ) has no singularity at x 1 = x 0 and may be interpreted as the bulk fermion propagating from point x 0 to the boundary, where it mixes with the boundary fermion, and then moving to point x 1 . 12 This is an explicit formula: Let us complete the calculation of G(∆, E) using Eq. (5.24). We have The area of the hyperbolic plane is obviously infinite, but it can be made finite by regularization. Indeed, consider the disk D r of radius r centered at the origin. It has the following area and boundary length: Hence, [Tr C sp ] reg = −2π Tr C sp (0, 0). Plugging this in (5.24), we get: . In conclusion, we rewrite Eq. (5.30) as follows, 32) and note that it is consistent with the combination of (2.66) and (2.34): Indeed, both equations give the same result for the mixed derivative if we use the fact that ∂(2θ)/∂∆ = −iπ tan π(∆ + iE) − tan π(∆ − iE) .
where Π ν λ is the projector onto the eigenspace of the SL(2, R) Casimir operator with the eigenvalue λ(1 − λ). The operators Π ν λ are defined by integral kernels that depend on pairs of points x 1 , x 0 ∈ H 2 ; the normalization is such that Π ν λ (x, x) = 1. We will make the connection to the Plancherel factor explicit by deriving (5.34) from (5.35), bypassing the full calculation of the Dirac propagator. As explained in appendix F, the components of a Dirac spinor have different effective spins ν, equal to ↓ = − 1 2 − iE and ↑ = 1 2 − iE. The Dirac operator is represented by the matrix where ∇ + and ∇ − are certain differential operators changing the value of ν by 1 and −1, respectively. (Here the subscripts "±" have nothing to do with boundary conditions.) The Casimir operator is expressed in terms of ∇ ± by Eq. (F.25), so both 4∇ − ∇ + for ν = ↓ and 4∇ + ∇ − for ν = ↑ are equal to 1 4 + E 2 − Q. Using this and the formula ∆ = 1 2 − √ M 2 − E 2 , we obtain the following expression for the propagator: Let us first calculate the matrix element involving ν = 1 2 − iE spinors with Dirichlet boundary condition (indicated by the subscript "+"), The general idea is to use the Casimir eigendecomposition (5.35); the role of boundary condition will become clear later.
For the task at hand, it is convenient to transform Eq. (5.35) to a different form, which generalizes to complex values of ν: where the contour Γ is illustrated in Fig. 9(a). It is obtained by a deformation of the vertical line Re(λ − ν) = 1 2 and consists of the line from 1 2 − i∞ to 1 2 + i∞ and circles surrounding the poles of tan π λ − 1 2 − ν in the strip 1 2 < Re λ < 1 2 + Re ν or 1 2 + Re ν < Re λ < 1 2 (depending on the sign of Re ν). The rewriting is based on this representation of the Plancherel factor, 40) and the symmetry Π ν λ = Π ν 1−λ , which allows one to extend the integral in (5.35) from a half-line to a full line. More explicitly, +∞ 0 ds s sinh(2πs) cosh(2πs) + cos(2πν) The discrete series contribution (i.e. the second term in (5.35)) can be treated as residues of the same integrand, which leads to the expression (5.39). Note that when λ and ν are arbitrary complex numbers, Π ν λ is no longer an orthogonal projector. Formally, it is just a function of x 1 , x 0 ∈ H 2 , and 1 ν should likewise be interpreted as a (generalized) function, namely, g(x 0 ) −1/2 δ(x 1 − x 0 ), where g is the determinant of the metric tensor. Given this caveat, we will proceed with caution. It is true that However, the following corollary holds only for the Dirichlet boundary condition and is qualified by a restriction on λ: Indeed, we should require that the left-hand side of the above equation be well-defined, meaning the absolute convergence of the corresponding integral: To check this condition, let us use polar coordinates, x = (r, ϕ). As r tends to 1, the propagator C ↑↓ + (x 1 , x) scales as (1 − r) 1−∆ , whereas Π ν λ (x, x 0 ) has terms proportional to (1 − r) λ and (1 − r) 1−λ . Since g(x) ∼ (1 − r) −2 , the convergence condition is exactly as indicated in Eq. (5.43).
We now apply the decomposition of identity (5.39) together with Eq. (5.43): Note that the contour Γ passes between the poles of the integrand at ∆ − = ∆ and ∆ + = 1 − ∆. The propagator C ↑↓ − with Neumann boundary condition cannot be obtained in the same way, but we can use analytic continuation in M . Suppose that M > |E| initially. As M changes to −M avoiding the branch cut between E and −E, the numbers ∆ + and ∆ − are swapped, and the propagator C ↑↓ + turns into −C ↑↓ − for the original value of M . On the right-hand side of (5.45), the analytic continuation should involve a deformation of the integration contour such that it avoids the moving poles, see Fig. 10. Thus, The "spooky propagator" C ↑↓ sp = C ↑↓ − −C ↑↓ + is given by the integral of the same function along the difference contour Γ sp ∼ Γ−Γ, which consists of circles wrapping the points λ = ∆ − (clockwise) and λ = ∆ + (counterclockwise) as shown in Fig. 10(c). Hence, the spooky propagator is determined by the residues of the integrand at ∆ − and ∆ + : The calculation of the other diagonal element of the propagator matrix, −C ↓↑ (in all three variants) is completely analogous; we just need to use ν = − 1 2 − iE. Restricting to coincident points and using the normalization condition Π ν λ (x, x) = 1, we obtain the final result: which is equivalent to Eq. (5.34).

Discussion
One of the main new physical consequences of our computations on the complex SYK model is the many-body density of states in Eq. (2.83). For each total charge Q, the energy dependence of the density of states is the same as in the Schwarzian theory, with a ground state energy E 0 (Q), and a zero temperature entropy S(Q/N ) determined by the value of Q. Although this result is natural from the physical point of view, we derived it from the effective action (1.12), which describes an ensemble with fluctuating Q. The presence of the particle-hole asymmetry parameter E in the action was essential for the consistency of that calculation. The other parameters in the effective action in Eq. (1.12) are the charge compressibility K, and γ, the coefficient of the T -linear specific heat at fixed Q. While the value of γ was determined by a low-energy analysis using conformal perturbation theory [6,7], we have shown here that a similar procedure does not apply for K. This is highlighted by the UV divergence in the eigenmodes of the symmetric sector of the two-particle kernel shown in Fig. 7. It is necessary to account for high energy contributions to obtain the correct value of K, and we presented three such computations in Sections 4.1, 4.2, and 4.3; the numerical values so obtained were consistent with each other. These distinct behaviors of γ and K are analogous to those in the Fermi liquid theory: the quasiparticle effective mass m * determines the specific heat, but an additional Landau parameter, F s 0 , is needed for the compressibility. We presented a new computation of the zero temperature entropy S of the complex SYK model in Section 5. The entropy was shown to be equal to the difference in the logarithm of the partition function of a massive Dirac fermion on H 2 between Neumann and Dirichlet boundary conditions, in a manner similar to the influence of double-trace operators in the usual AdS/CFT correspondence [40][41][42][43][44][45]. This bulk approach correctly reproduced the Q and E dependence of S in the SYK model.
The above computation of the entropy should be contrasted with that in higher dimensional black holes whose near-horizon geometry has an AdS 2 factor (reviewed in Appendix G), where the entropy is given by the horizon area in the higher-dimensional space, and arises from degrees of freedom unrelated to the fermions. While all entropies mentioned here obey Eq. (1.16), the functional form of S(Q) is different for the higher-dimensional black holes [3]. Probe Dirac fermions can be added to such higher-dimensional black holes [5,48], and their Green function agrees with those of the SYK model [9,3]; however such fermions only contribute O(1) entropy in the distinct large-N limit of the usual AdS/CFT correspondence.

A Luttinger-Ward analysis and the anomalous contribution to charge
In this section, we will discuss frequency domain derivations of the charge formula (2.34) for general q following the strategy in Ref. [4] (GPS) appendix A. Here we aim to provide an alternative route to the discussions in Section 2.2 that may be more transparent to the readers familiar with Luttinger's theorem and Luttinger-Ward functional. We also draw attention to the comparison with perturbative anomalies in quantum field theory.

A.1 IR divergence and anomaly
Instead of Feynman propagator used in Ref. [4], we will work with the imaginary time Green function for convenience and express the charge as the following integral We proceed by the standard Luttinger-Ward procedure (see Ref. [30]), i.e. inserting the identity 1 = ∂ z (G(z) −1 + Σ(z)), which leads to the expression This manipulation is similar to the manipulations done in Eq. (2.28). However, instead of further anti-symmetrizing the integrand, we split the two terms in braces with an explicit cut-off: where the principal value is implemented by a symmetric cut-off in frequency-domain: We emphasis that the regularization is crucial in the discussion here: the integral I 1 and I 2 are logarithmically divergent, so their value depends on the regularization-scheme. More explicitly, the logarithmic divergence arises from the IR asymptotics of the Green function (non-Fermi liquid behavior) G(z) ∼ z 2∆−1 . On the contrary, the analogous integrals in the standard Luttinger-Ward analysis for the Fermi liquid are well defined (i.e. with no divergence) and one can further prove the second integral actually vanishes due to the existence of the Luttinger-Ward functional [30]. The situation here is very similar to the perturbative anomalies in quantum field theory (e.g. see the discussions in [49] chapter 19, a similar comment was also made in Ref. [50]). A particularly simple example is the two dimensional massless QED, where the Feynman diagrams formally satisfy the Ward identity both for vector and axial current. However, the regularizations will make it impossible to have both gauge invariance and the axial current conservation.
In this section, we choose to use a regulator (A.4) (following GPS) that let I 1 term inherit the physical meaning as in the Fermi liquid, while let I 2 term carries the anomalous contribution. As a comparison, the time domain symmetric regulator used in Section 2.2 will set I 2 = 0 but shift the value of I 1 integral. In any case, the sum I 1 + I 2 is regularization-scheme independent and determines the physical charge.

A.2 Calculation of I 1 integral
Let us first evaluate the I 1 integral explicitly In the second step we bend the contour close to the real axis so that we can proceed using the analytic properties of Green function, thus (arg G(∞ + iη) − arg G(iη)) . (A.6) A.3 Anomalous Luttinger-Ward term I 2 at q = 4 Now, we calculate I 2 integral in the present regularization-scheme. Before moving to the evaluation of I 2 for general q, let us first review/simplify and remark on the detailed calculations performed in Ref. [4] appendix A for q = 4 model. In the reference aforementioned, I 2 is expressed using spectral function: where the domain {ω} is defined as {ω} := {ω 1 , ω 2 > 0, ω 3 < 0} ∪ {ω 1 , ω 2 < 0, ω 3 > 0} 13 . The spectral function ρ(ω) = − 1 π Im G R (ω) has the following IR asymptotics The UV behavior of ρ(ω) has to be determined by numerics. However, I 2 integral only depends on the IR asymptotics and therefore universal. A simple argument is as follows. Without the cutoff, the z integration in Eq. (A.8) will run into a logarithmic divergence at small frequency, which can be seen by power counting of the IR asymptotics of ρ(ω). Now assume we consider a variation of the spectral function δρ(ω) that does not change the IR asymptotics of ρ(ω), e.g. δρ(ω) ∼ ω −1/2+s with s > 0 at ω → 0. Then the corresponding variation of the integral δI 2 is free of IR divergence as the ω integrations contribute a term asymptotic to z −1+s at small z and the z integration is IR finite now. Therefore, for the variation δI 2 , there is no obstruction to take the η → 0 limit first, namely replacing the principal value by an integration along imaginary axis. Then we can integrate z first by deforming the contour to right half plane and picking up the residues: where x = ω 0 and y = ω 1 + ω 2 − ω 3 in I 2 . Next, we finish the d 4 ω integration and have That is to say, for a variation δρ(ω) that does not change the IR asymptotics of spectral function, the integral I 2 is unchanged. This conclusion also allows us to substitute the exact ρ(ω) by its IR form while calculating I 2 . Thus, where s ± are defined in Eq. (A.9) and characterize the spectral asymmetry. Finally, we evaluate the explicit integrals, first for d 4 ω and then dz: We may call I 2 the "anomalous" Luttinger-Ward term as it arises from a formally vanishing integral. As we mentioned before, its counterpart in Fermi liquid is well-defined and indeed vanishes. The anomaly discussed here also shares some similarity with the perturbative anomaly in quantum field theory.

A.4 Dimensional Regularization
It may be useful to also present a "dimensional regularization" version of the calculation for the anomalous term I 2 . More explicitly, we use the following form for Green function where η is a small positive number that will be taken to zero in the end. Note that the regulator here is symmetric in ω. 14 In other words, the present regulator has the same symmetry as the hard cut-off used above and they should give the same value for I 2 .
The small shift in the scaling of Green function will induce a small shift in both the scaling and the prefactor in the spectral function: The shift in the prefactor follows from the analyticity of Green function G η (z) on the upper half plane. On the other hand, the shift in scaling saves the I 2 integral from IR logarithmic divergence and allows us to do the dz integration first. Therefore, We can proceed by inserting the explicit expressions for ρ and ρ η with an intermediate transition scale ω Λ above which the integrand identically cancels, (A.17) 14 In contrast to the one used in Section 2.2 which is symmetric in τ .
The actual value of ω Λ is not important and will not enter the final result. Now the dω integration are straightforward to perform. After taking the limit η → 0, we have 18) A.5 Generalization to q > 4 It is straightforward to generalize the explicit calculations to q > 4 using the same regularizationscheme. First of all, the argument that the I 2 term only relies on the IR asymptotics of the spectral function (A.9) still applies. Therefore, we end up with an explicit integral that generalizes Eq. (A.12): (A.19) It is useful to note that the d q ω integration is the "multivariate Beta function". Or more explicitly, we can use the following formula Thus, the integral simplifies and we insert the expressions for s ± to get (A.21) Putting together with I 1 , we reproduce the charge formula A.6 q → 2 limit and the semicircle We would like to check if the charge formula (A.22) is valid at q → 2 limit. Namely, whether is consistent with the result obtained by directly diagonalizing the quadratic random Hamiltonian. For q = 2, we can solve the Schwinger Dyson equation Figure 11: Spectral function ρ(ω) for the q = 2 model has the shape of semicircle as expected (after rescaling the x or y axis).
exactly with the UV term. We work in zero temperature and have recovered the J dependence in the equation. The solution is given by where the sign is determined by the requirement that the spectral function is non-negative and the Green function in the lower half plane is the complex conjugate of the upper half plane.
The spectral function is nonzero on an interval −2J − µ < ω < 2J − µ and has the semicircle form with area 1 as expected and shown in Fig. 11. The charge is determined by the following simple formula where the −1/2 is the shift required by the definition. Pictorially, Q is equal to the area of ABOC in Fig. 11. Next we determine phase factor θ by matching the IR asymptotics Fig. 11) .
Now we are ready to relate the charge Q to θ using the semicircle geometry shown in Fig. 11 Q = area of ABOC = 2 π

B Operator spectrum
The solutions of the equation det(1 − K G (h)) = 0 contain important information about the OPE of two fermions ψ † j (τ ) ψ j (0). For instance, at θ = 0, the matrix is symmetric, and therefore, the eigenvectors of Using the reflection symmetry h ↔ 1 − h (cf. Eq. (3.21)), we restrict our discussion to h 1/2 and label them in ascending order in this section, i.e. 1/2 h . . are solutions of k A/S (h) = 1 respectively. In particular, the anti-symmetric sector, corresponding to the solutions of k A (h) = 1, reproduces the scaling dimensions of the operators appearing in χ j (τ ) χ j (0) OPE of the Majorana SYK (which is determined by the equation k c (h) = 1 in the notation of Ref. [6,7]). The leading one h A 0 = 2 corresponds to the Schwarzian sector and responsible for the energy fluctuation. Analogously, the leading one in the symmetric sector h S 0 = 1 is related to the U(1) charge in the complex SYK model as we discussed in Section 3. For general θ, the matrix W Σ (h) has no symmetry and the symmetric/anti-symmetric sectors generally mix via 2 × 2 ladder kernel K G (or K Σ ). Let us denote the two eigenvalues by k A (h, θ) (anti-symmetric branch) and k S (h, θ) (symmetric branch) as a generalization of the notation k A/S (h). Their explicit formulas are as follows, To illustrate, we plot k A (h, θ) (blue lines) and k S (h, θ) (red lines) as functions of h for ∆ = 1/4 and θ = π/8, π/6, π/4 together with θ = 0 (black lines, as a reference) in Fig. 12. Here are some comments: Figure 12: Plots of k A (h) and k S (h) for ∆ = 1/4 and θ = π/8, π/6 and π/4. Black lines are for reference, represent the value of k A/S (h, θ) for θ = 0.

For integer value
An immediate corollary is that k A (2, θ) = 1 and k S (1, θ) = 1 for all θ, i.e. the scaling dimensions of the energy and charge operator are protected (as well as their dual field with h = −1 and 0 respectively).
2. General solutions of k A/S (h, θ) = 1 depend on θ. In Fig. 13 we plot the value of h A/S 1 (θ) as functions of θ for ∆ = 1/4. Fig 12, we notice that for ∆ = 1/4, there is a critical value θ c = π/6 which is determined by the following equation

From
Above the critical value, namely for θ > θ c , the solutions h S n 1 disappear. In other words, the only solutions for k S (h, θ θ c ) = 1 are h S 0 = 1 and its dual 0. To explain why Eq. (B.6) is relevant, let us analyze the pole structure of k S (h, θ). Naively the expression (B.5) has simple poles at h = 2∆ + m with m ∈ Z 0 due to the overall factor Γ(2∆ − h). However, the expression in big parentheses in (B.5), has zeros that cancel some of the poles. Indeed the poles at odd m, i.e. at h = 2∆ + 1, 2∆ + 3, . . . are canceled in k S (h, θ). Furthermore, when cos(2θ) < 1 − 2∆, the poles at even m i.e. at h = 2∆, 2∆ + 2, . . . are also canceled. At critical value cos(2θ) = 1 − 2∆, there is a discontinuity for k S (h, θ c ) at h = 2∆ + 2k for k 1. Explicit calculation yields (for ∆ = 1/4) . (B.8) Parallel discussions apply to the anti-symmetric branch k A (h, θ) where an additional set of solutions to the equation k A (h, θ) = 1 emerges when θ > θ c , as shown in Fig. (12.c). Technically, this is related to the additional set of poles at h = 2∆ + 2k.

C Low energy contribution to K
Our computations of fluctuations in this appendix, and in Section 4.3, follow the methods of Ref. [6]; these methods are related to those in Section 3. The key difference is that in Section 3 we expand the (G, Σ) action w.r.t. the conformal saddle, while in this appendix and Section 4.3 we expand (G, Σ) action w.r.t. the exact saddle. We will use subscript " c " and " exact " to emphasis the contrast. As we discussed in Section 4, we expect that all energy scales contribute to the numerical value of the compressibility K, and a low energy conformal perturbation approach (similar to that used successfully for the specific heat in Ref. [6]) does not yield the correct value of K, instead it reproduces K linear discussed in Section 2.4.2 (cf. Eq. (2.74)). The UV divergence of the symmetric kernel eigenmodes in Section 4.3 provides explicit evidence for this claim. Nevertheless, we will present the low energy analysis of the symmetric kernel here, as it could be useful for other investigations.

C.1 Effective action for fluctuations around the saddle point
In this section we consider the (G, Σ) action for the complex SYK model with zero chemical potential and derive effective action for quadratic fluctuations around the saddle point of the action. In this section we recover J.
The (G, Σ) action for the complex SYK with zero chemical potential is (cf. Eq. (2.1)) The crucial difference from the Majorana SYK model is that now the bilocal fields G(τ 1 , τ 2 ) and Σ(τ 1 , τ 2 ) are not necessarily antisymmetric under exchange of variables τ 1 ↔ τ 2 . The saddle point G exact , Σ exact of the action (C.1) is the exact solution of the Schwinger-Dyson equations (2.2). Now we consider small fluctuations around the exact saddle point G exact , Σ exact : G(τ 1 , τ 2 ) = G exact (τ 12 ) + δG(τ 1 , τ 2 ), Σ(τ 1 , τ 2 ) = Σ exact (τ 12 ) + δΣ(τ 1 , τ 2 ) , (C.2) and expand the (G, Σ) action up to quadratic terms. Next we integrate out Gaussian fluctuations of the δΣ field and obtain the Gaussian action for the fluctuations δG, which is convenient to parametrize as δG(τ 1 , τ 2 ) = |G exact (τ 12 )| 2−q 2 g(τ 1 , τ 2 ) where we also decomposed fluctuations g(τ 1 , τ 2 ) on symmetric and antisymmetric parts g(τ 1 , τ 2 ) = g S (τ 1 , τ 2 )+g A (τ 1 , τ 2 ), so g S/A (τ 2 , τ 1 ) = ±g S/A (τ 1 , τ 2 ) and introduced the antisymmetric and symmetric kernels K A/S exact whose explicit expressions have been shown in Eq. (4.7). We copy the formulas here with the emphasis of subscript " exact " (C.4) In the large βJ limit the exact Green function G exact in the kernels K A exact and K S exact can be approximated by the conformal solution G c , so one obtains conformal kernels K A c and K S c . The spectrum of the conformal kernels can be computed exactly and was given in (B.3) where h labels the SL(2, R) representations: for the antisymmetric case h = 2, 4, 6, . . . and for the symmetric case h = 1, 3, 5, . . . and there is also a principal series h = 1/2 + is, s ∈ R + for both cases. Note k A c (h = 2) = 1 and k S c (h = 1) = 1, i.e. if one uses conformal kernels in (C.3) then the effective action for these special modes is zero, which would also indicates instability [53]. Actually, this problem appeared because we replaced the exact kernels by conformal ones. The exact eigenvalues k A and k S which correspond to h = 2 and h = 1 modes differ from 1 by 1/(βJ) corrections. In order to find these corrections one uses 1/(βJ) correction to the conformal Green function [6,7] and α G is a UV dependent constant, which can be found numerically. Next using the corrected Green function in the kernels one finds corrections to their eigenvalues.
We find the 1/(βJ) correction to the conformal kernels using perturbation theory [6]. At the first order one computes diagonal matrix elements of the perturbation. In our case the perturbation to the conformal kernels consists of two parts and reads The part of V which involves f 0 (τ 12 ) is called the rung term and the part with f 0 (τ 13 ) is called the rail term. The corrections to the eigenvalues are simply where for |Ψ A and |Ψ S we take unperturbed conformal eigenfunctions, so K A c |Ψ A = k A c |Ψ A and K S c |Ψ S = k S c |Ψ S . For the antisymmetric case the correction for h = 2 mode was already found in [6,7] and reads (k A = 1 + δk A ) In the next subsection we are going to compute correction to h = 1 eigenvalue of the symmetric kernel. We expect to find a similar to (C.9) form, but with some other coefficient α S K . We will find that such obtained value α S K does not agree with the numerical computation in the section 4.3. The reason for this disagreement is hidden in the fact that when we use conformal Green function in the kernels instead of the exact one, we implicitly assume that the eigenvalues k A (h = 2) and k S (h = 1) are not affected by the UV domain (τ < 1/J or τ > β − 1/J), where the conformal Green function G c diverges, but the exact G exact goes to 1/2. And in fact this turned out to be correct for the antisymmetric case. But it is not correct for the symmetric case. On a technical level the exact h = 1 eigenfunctions Ψ exact 1,n (τ 1 , τ 2 ) of the symmetric kernel do not approach completely conformal eigenfunctions Ψ 1,n (τ 1 , τ 2 ) at the βJ → ∞ limit. There is always a discrepancy in the UV domain. The exact eigenfunctions grow as βJ at the coincident points τ 1 = τ 2 , whereas the conformal ones approach a constant. This effect is nicely captured in the large q limit, which we discuss in the appendix D.

C.2 Symmetric sector
The rung and rail integrals in the symmetric sector are (we omit factor − 2α G βJ for brevity) where the conformal kernel, h = 1 wave functions and the correction to the conformal propagator are , where we fix β = 2π so all angles take values in the interval [−π, π]. It is convenient to represent different integrals by Feynman diagrams. Let us introduce Feynman rules. We denote propagators as A very useful tool for computation in a conformal theory is the star-triangle identities [54] +π −π dθ 0 (sin 2 θ 01 2 ) α 1 (sin 2 θ 02 2 ) α 2 (sin 2 θ 03 where α 1 + α 2 + α 3 = 1 and . (C.14) Graphical representation of the star-triangle identities is 15) In addition to the star-triangle identities we also list a couple of useful integrals. One is the integral, which gives the delta-function π −π dθ 0 sgn(θ 10 )sgn(θ 02 ) (sin 2 θ 10 2 2α). The other useful integral, where α 1 , α 2 , α 3 are arbitrary real numbers is The integral (C.17) can be computed by setting one angle to zero and projecting to a line, where one can use Feynman parameters.
Since the correction to the conformal propagator f 0 (θ 12 ) is obtained from three point function of two fermions with the operator of dimension h = −1, it can be represented as the integral Since the integrals δk S rung and δk S rail are logarithmically divergent we introduce a soft cutoff η → 0 by multiplying f 0 (θ) by (sin 2 θ 2 ) η , so the new graphical representation for f η 0 (θ 12 ) is Finally we notice that f 1−∆ f ∆ = −4(q − 1)α 0 and the conformal kernel can be depicted as

C.3 Computation of the rung integral
In case of the rung integral we can simply integrate over θ 3 and θ 4 to obtain Next, it is convenient to use decomposition which can be derived from multiple-angle formula and properties of the Chebyshev polynomials. Therefore the integral (C.22) takes the form (in what follows we assume that n > 0, so |n| = n).
For k = 1 we have c 1 = n 2 and we find using (C.18) For k > 1 there is no divergence and we can set η = 0 and compute n k=2 where H n = n m=1 1 m is the Harmonic number. Therefore we finally find

C.4 Computation of the rail integral
In this case we represent h = 1 eigenmodes Ψ S 1,n (θ 1 , θ 2 ) in the form .
(C. 41) We notice that this can be written as where k A c (h) is defined in (C.5). So finally we find for the rail integral Combining the rung and rail integrals and restoring the factor − 2α G βJ we find for the corrected symmetric kernel eigenvalue Using equations (4.12), (C.9) and (C.44) we find where γ = 4π 2 α S . Unfortunately, the result (C.44) gives only a partial contribution for corrected symmetric kernel eigenvalue at order 1/βJ , and (C.45) agrees with the K linear value in (2.74). Presumably integrals obtained from higher corrections to the Green function have uncompensated power-law divergences which then contribute to the first order 1/βJ term, changing the final result.
We also notice that the computation we did can be drastically simplified if we take large n limit as in [6]. Doing this we indeed find the result (C.44). The caveat of such approach in our case is that a priori we can not guarantee that the final result contains only linear in n terms, because each integral is divergent, so in principle we could miss some other corrections in n.

D Large q for symmetric kernel
In the large q limit, we take q → ∞, keeping J = √ qJ/2 q− 1 2 fixed. Then in the first order the Green function reads G(θ) = sgn(θ) 2 1 + g(θ) q + . . . , e where we set β = 2π and v is a function of βJ , which is found from the equation πv/ cos πv 2 = βJ , so v → 1, when βJ → ∞. Now consider the symmetric kernel Using the large q Green function (D.1) we find in the leading q order In this expression the dependence on q remained only in 1/q prefactor, the rest depends only on βJ , which is our parameter. Therefore, evidently, the eigenvalues of the large q kernel In the leading order in βJ this gives for the kernel eigenvalues k S q=∞ = We compare the large q wave functions (D.12) and numerical results for h = 1 and n = 3 in Fig. 14.

E Zero temperature entropy and ln det term
For the Majorana model at large N , the free energy can be obtained by extremizing the (G, Σ) action: Inserting the solution of the Schwinger-Dyson equation G and Σ, we have . We would like to calculate F to the first order in β −1 and extract the linear coefficient S, namely the zero temperature entropy. Moreover, the linear term is only present in F 1 , not in F 2 . One way to see this is to calculate the F 2 integral in the conformal limit β → ∞, and try to extract the finite piece in βF 2 : The integral has a UV divergence that contributes to the ground state energy F 0 . We are interested in the finite piece, which can be obtained after a regularization. For instance, one can use a cut-off and evaluate the integral which has vanishing constant piece. The absence of even power in terms in the above integral is due to the θ → −θ symmetry of the integrand. Thus, we conclude that the zero temperature entropy is only contained in F 1 , i.e. the ln det term.

(E.5)
More explicitly, we consider the partial sum with a cutoff n Λ ∼ β → ∞ and single out the n Λ independent term as the zero temperature entropy The analytic formula for the finite piece S turns out to be easier to be found by evaluating its ∆ derivative S (∆), which amounts to summing digamma functions ψ(x) = Γ (x)/Γ(x) by the following formula Finally, integrating S (∆) to the desired position with boundary value, we have the entropy formula for the Majorana model .
This procedure has been described in [39,2]. The emphasis here is to support the claim S = 1 2 ln det(−∂ τ − Σ), which will be given a bulk interpretation.

F.1 Dirac operator and spinors in two dimensions
In a Lorentzian space, the Dirac Lagrangian is We will work in Euclidean signature. In the case of flat space, the Wick rotation takes x 0 to x 2 = ix 0 . The Dirac matrices γ c and spin matrices Σ ab = 1 4 γ a , γ b are (Λ 0 represents an infinitesimal counterclockwise rotation.) The Euclidean action is The Majorana case differs in that ψ ↓ , ψ ↑ are real and that the action has an overall factor 1 2 . Note that the Dirac spinor ψ splits into two irreducible representations of the universal cover of SO(2) (or the Lorentz group): ψ ↓ has spin − 1 2 , and ψ ↑ has spin 1 2 . A general ν-spinor ξ is one-dimensional and transforms as Λ 0 ξ = −iνξ. In the absence of electromagnetic field, the covariant derivative may be written as Here, (ω αbc ) is a spin connection defined relative to some local orthonormal frame (vielbein), whereas (ω α ) may be regarded as a vector potential such that To take advantage of the splitting of the tangent space under SO(2), let us replace the local orthonormal frame ( e 1 , e 2 ) with e + = 1 2 ( e 1 − i e 2 ) , e − = 1 2 ( e 1 + i e 2 ) .

F.2 Spinors on H 2 and SL(2, R) symmetry
The hyperbolic plane H 2 is described by the Poincare disk model with the metric (F.16) The local frame is chosen to be proportional to the coordinate frame: This choice is called the "disk gauge" in Ref. [47]. It is a special case of conformal gauge, which is defined for any metric in the conformal form, ds 2 = f (z, z) dz dz. In this gauge, (F.18) The operators ∇ + , ∇ − commute with isometries of the underlying manifold. Let us introduce modified polar coordinates (u, ϕ) such that A ν-spinor with angular momentum m is proportional to e i(m−ν)ϕ (in the disk gauge), hence, We now discuss SL(2, R) symmetry following Ref. [47]. The abstract symmetry generators L −1 , L 0 , L 1 , have two different realizations. The more natural one is by Killing vector fields, acting on spinors as follows: The second set of operators (commuting with the first) is denoted by L R n ; they change the spin by −n: and that the commutation relations between these operators are just a special case of (F.12). This is an explicit expression for the Casimir operator: The eigenvalues of Q may be parametrized as ∆(1−∆). The joint eigenspace of this and the angular momentum operator, m = −i∂ ϕ + ν, is spanned by the ν-spinors ξ ν ∆,m , ξ ν 1−∆,m defined by the formula We will also need the asymptotics at the origin, where one term usually dominates. In the special case m = ν, we have ln |z| for z → 0 . (F. 30) Here and below, Z stands for a point in H 2 with coordinate (z, z). This way, we distinguish general functions like ξ ν ∆,m from analytic functions of z. The operators L R ±1 act on the basis spinors as follows:

G Higher dimensional black holes in asymptotically AdS space
This appendix will begin by recalling the basic thermodynamic properties of charged spherical black holes in global AdS d+2 [37,5] in d + 2 spacetime dimensions (d 2). We denote the T = 0 radius of the black hole by R h . Then we will discuss the universal structure of the theory of such black holes at temperature T 1/R h [3,[13][14][15][16][17][18], where the effective action in Eq. (1.12) applies. See Ref. [17] for more details.
In the AdS/CFT correspondence, AdS d+2 spacetimes are dual to conformal field theories (CFTs) in d + 1 spacetime dimensions. With an Einstein-Hilbert gravitational action, the CFT is the large N maximally supersymmetric SU(N ) Yang-Mills theory for d = 3, and a suitable large N limit of the ABJM theory for d = 2. We place the CFT on a sphere, S d , and add a chemical potential conjugate to a global U(1) symmetry. This induces a total charge Q on S d . In the holographic description, the asymptotically AdS d+2 spacetime crosses over to a charged black hole with a with a near-horizon AdS 2 × S d geometry [37,5,56], which was identified with the physics of SYK models [9].
The Einstein-Maxwell theory of a metric g and a U(1) gauge flux F = dA has Euclidean action where κ 2 = 8πG N , R d+2 is the Ricci scalar, L is the radius of AdS d+2 , and g F is a gauge coupling constant. The properties of the black hole are fully specified by the temperature T and the chemical potential µ. The later is specified by a boundary condition on the time component of the U(1) gauge potential A lim r→∞ A τ (r, τ ) = iµ .

(G.2)
At T = 0, let the radius of the black hole horizon equal R h , the total charge equal Q, and the chemical potential equal µ 0 . These quantities are related to each other by where s d is the area of the d-dimensional surface of a unit sphere. We will treat R h as the independent variable below, the dependence on Q and µ 0 = µ(T = 0) follows from the above.
Moving to non-zero T 1/R h , we find the entropy S(Q, T ) has the form in Eq. (1.14) (we do not track factors of N in this Appendix) with .

(G.4)
Note that the entropy is simply given by the area of the horizon in the higher-dimensional geometry. The contribution of fermion determinants to the action is subdominant in the large N limit of the AdS/CFT correspondence, unlike the computation in Section 5.
The low T behavior of the chemical potential is given as follows, . (G.6) We can now verify that the Maxwell relation is obeyed as T → 0, which then implies the fundamental identity in Eq. (1.16). Finally, we also note the value of the compressibility (G.8) Now we turn to the universal structure for T 1/R h . The near-horizon metric takes the AdS 2 × S d form with metric and gauge field (at T = 0, see Ref. [17] for T > 0) (G. 10) and the dimensionless parameter E determining the strength of the near-horizon electric field is the same as that in Eq. (G.6). The Green function of a massive Dirac fermion at the AdS 2 boundary was computed in Ref. [5,48], and was found to have the same spectral asymmetry as that in Eq. (1.7), also determined by E. Several works [9,3,[10][11][12]8,[13][14][15][16][17][18] have discussed the nature of the theory of AdS 2 horizons, which is applicable to the higher dimensional black holes at T 1/R h . Under these conditions, all modes which are non-constant on S d can be ignored, and we can focus on the fluctuations of the AdS 2 sector. These fluctuations are also described by the 0 + 1 dimensional Schwarzian theory [10][11][12]8] in Eq. (1.12). A recent analysis [17] has shown that the parameters K, γ, and E appearing in this Schwarzian theory are precisely those specified by the thermodynamic results recalled in this appendix.