Subsystem R\'enyi Entropy of Thermal Ensembles for SYK-like models

The Sachdev-Ye-Kitaev model is an $N$-modes fermionic model with infinite range random interactions. In this work, we study the thermal R\'enyi entropy for a subsystem of the SYK model using the path-integral formalism in the large-$N$ limit. The results are consistent with exact diagonalization [1] and can be well approximated by thermal entropy with an effective temperature [2] when subsystem size $M\leq N/2$. We also consider generalizations of the SYK model with quadratic random hopping term or $U(1)$ charge conservation.


Introduction
In these years, there is an increasing interest for the entanglement/thermal entropy of manybody systems in both condensed matter and high-energy physics. On the condensed matter side, entanglement/thermal entropy can be used to identify the critical theory [3][4][5][6] as well as topological ordered phases [7,8]. It can be used to distinguish the many-body localization (MBL) phase [9][10][11] and the thermal phase [12,13]. In MBL phase, the entanglement entropy of the eigenstate has area law scaling [14,15], while in the thermal phase, the reduced density matrix of the typical state takes a thermal form and therefore the entanglement entropy obeys volume law scaling [16]. Recent developments for the study of the information scrambling also focus on the entanglement entropy dynamics [17][18][19][20][21][22][23][24]. As an interesting example, the projective measurement can induce a transition from the entanglement entropy perspective in systems with random unitary evolution [25][26][27]. On the high energy side, the discovery of Ryu-Takayanagi formula [28][29][30] directly relates the entanglement/thermal entropy of a subsystem to holographic bulk geometry, which has been latterly generalized to the time dependent case [31] or with quantum corrections [32,33]. Recent refinement of the Ryu-Takayanagi formula gives an plausible solution to the information paradox [34][35][36][37][38][39].
However, for general interacting many-body system without holographic description or conformal symmetry, there is no efficient method to study the behavior of entropy, especially for strongly correlated systems. Fortunately, a solvable model describing N Majorana modes with infinite range random interaction has been proposed [19,40], which is related to the early work of Sachdev and Ye [41]. This model is now known as the Sachdev-Ye-Kitaev (SYK) model. The studies show that the SYK model is a non-Fermi liquid without quasiparticles and has emergent conformal symmetry at low energies [19,40]. The low energy Schwarzian action [19,40,42] also matches the result of NAdS 2 dilaton gravity [42,43]. Latterly, different generalizations have been proposed, including a complex fermion version [44,45] and the SYK model with extra random hopping terms [46][47][48], which have important applications in different fields. The hallmark of these models is that the computation of correlation functions is reduced to integral equations in the large-N limit. This also leads to the ability to compute quantum dynamics precisely [49][50][51][52][53].
In this work, we study the SYK model and its generalizations from the subsystem entanglement/thermal entropy perspective. We use the path integral method and derive the saddle point solution of the second Rényi entropy by taking the large N limit. We solve the saddle point equation numerically and obtain the scaling of the second Rényi entropy in various cases. In particular, we compare the results at extremely low temperature with the entanglement entropy result of SYK model studied by using exact diagonalization for finite size systems [1] and using the subsystem thermalization arguments [2] 1 . We find the results are consistent with the subsystem size is smaller than half of the system size.
The rest of the paper is organized as follows: In Section 2, we briefly review the SYK model and the subsystem thermalization arguments. In Section 3, we recall the path-integral formalism [1] for the Rényi entropy, where the action is written in terms of bi-local fields. In the large-N limit, we derive the saddle point solution of the Rényi entropy and analytically compute it in various limits. In Section 4, we solve the saddle point solution and present the numerical results for the second Rényi entropy of SYK models with different subsystem size and temperature. We find the argument in Section 2 works well despite small deviation from the analytical approximation [2]. We further study the generalizations of SYK model by introducing random hopping quadratic term or considering the complex SYK model with charge conservation law. We conclude our results in the Section 5. 1 There are also studies of entropy dynamics for coupled SYK models prepared in an thermofield double state [39,54,55], which is related to the black hole evaporation problem.

Review of the SYK model and its entanglement
The SYK q model describes q-body random interacting Majorana fermions, the Hamiltonian is written as: Here q is an even integer. i = 1, 2...N labels different Majorana modes with commutation relation {χ i , χ j } = δ i j . The J i 1 i 2 ...i q with different indices are independent Gaussian variables with zero means. Their variance is given as: In the large-N limit, the two-point function in thermal equilibrium G β (τ) = 〈T χ i (τ)χ i (0)〉 is computed by taking melon diagrams [19,40]. The Swinger-Dyson equation is written as Here ω n = 2π(n + 1/2)/β is the Matsubara frequency. In the low energy limit, by neglecting the iω n term, the Green's function (at J −1 τ β) is found to be proportional to sgn(τ)/|τ| 2∆ with ∆ = 1/q [40]. Here ∆ is identified with the scaling dimension of fermions. As a result, the spectral function for small frequency diverges as ω 2∆−1 for q ≥ 4. This is understood as a non-Fermi liquid behavior. For q = 2, the model is a quadratic random hopping model, and the single particle spectral shows a semi-circle law with finite weight as ω → 0. The free energy, or the on-shell action, is determined by the solution of (3) [40]: The thermodynamical entropy S and the energy E are then determined by taking derivatives. It is known that for q ≥ 4, there is extensive entropy at zero temperature limit [19,40]: which reflects the existence of exponential number of low-energy states in many-body spectral.
We have S(1/2) = 0, S(1/4) ≈ 0.2324N and S(0) = log(2)N /2. Now let's review the argument [2] for computing the subsystem entropy of the SYK model. We consider the full system is prepared in a state described by a density matrix ρ. Here ρ can be either a pure state or an ensemble. To define the entropy of a subsystem, we divide the system into A and B, where A contains M Majorana fermions and B contains N − M Majorana fermions. We take the large-N limit with fixed ratio λ = M /N . The n-th Rényi entropy of the subsystem A is given by The Von Neumann entropy is given by taking n → 1.
The energy E of the total system is given by E = tr(ρH). For the subsystem A, its Hamiltonian H A is defined as restricting H into the subsystem A: This is again an SYK model, but with effectiveJ = λ q−1 2 J because of the N dependence in (2). If we assume the ρ is symmetric under the permutation of indices. The energy of the subsystem is given by E A = M q N q E. The energy density in A is then E A /M = λ q−1 E/N . Using E 0 to denote the ground state energy of corresponding system, we have: Assuming the system thermalizes [56][57][58], the proposal [2] is to approximate ρ A by a thermal ensemble with with energy E A . For pure states, this approximation can only work if λ < 1/2. Then we have s is the entropy density at relative energy density ε A in the thermodynamic limit. As an approximation which is good away from ε A = 1 [63]: In the limit of λ → 0, the energy density becomes zero, which corresponds to a maximal mixed state with S (n) A = M log(2)/2. This works for any λ < 1 if we take q → ∞.

The formulation
In this section, we consider the path-integral representation for computing the subsystem Rényi entropy [1], restricted to a thermal ensemble ρ = e −β H /Z of the total system. Here We firstly consider the reduced density matrix ρ A . A pictorial representation of the density matrix ρ is given by: Here the free end of the solid line denotes the possible quantum states. Imposing the boundary condition on the end gives the matrix element of density matrix. The dashed line represents the interaction between subsystems. To find the reduced density matrix for system A, we trace out the subsystem B: Here the black dot represents the anti-periodic boundary condition for fermionic fields [64].
Then tr A ρ n A can be computed by sewing n reduced density matrix with an overall anti-periodic boundary condition. To be concrete, from now on we focus on the n = 2 case, which gives a lower bound of the Von Neumann entropy. We have: We then have S The generalization to the n-th Rényi entropy is straightforward. More explicitly, we have Here the boundary condition is given by Note that these boundary conditions break the time translation symmetry of the system explicitly.
To proceed, we should average over the disorder configuration of J i 1 i 2 ...i q as S However, this computation needs to introduce additional disorder replicas. In the large-N limit, both numerical analysis and analytical arguments suggest the validity of disorder replica diagonal ansatz for the SYK model to the leading order of 1/N expansion [42,45,65,66]. We assume that similar arguments apply here and make the approximation: It is then straightforward to integrate out J i 1 i 2 ...i q and introduce the G − Σ action according to the standard procedure [40]. The main difference is that we should introduce two sets of bilocal fields, in order to take different boundary conditions into account. Explicitly, we define Additional Lagrangian multiplier is introduced into the path integral as: After further integrating out the Majorana field, one obtains the G − Σ action for computing the subsystem entropy: with an action Here the label under the det indicates the boundary condition (14). The action is proportional to N . Consequently, we could use a saddle point approximation to the leading order of 1/N . The saddle point equation is Diagrammatically, this again corresponds to melon diagrams, while internal lines can be either in A or B, with the corresponding weight. Note that although the self-energy is the same for both A and B, the inverse for computing G A and G B is under different boundary condition, which leads to different results. After solving these equations self-consistently, we can then compute the action (19) and obtain the entropy S A . Finally, although we use the original SYK model as an example, it is straightforward to generalize the derivation to any SYK-like models with an effective G − Σ action 3 . We will present two examples in the next section.

Analysis in different limits
For λ → 1, all Majorana modes are in B. Consequently, the path-integral reduces to a thermal ensemble with inverse temperature 2β and the final result S (2) A reproduces the thermal Rényi entropy of the full system. A more non-trivial limit is λ → 0. For λ = 0, the contribution from I (2) and I β cancels, leading to S A = 0 as expected. To get the leading order contribution for small λ, we need to expand the action (19) around λ = 0 and then computing the action using the λ = 0 solution of Green's functions (20). The action to the leading order of λ is written as: Here we have used the Schwinger-Dyson equation (20) at λ = 0. To proceed, we need to know the specific form of the Green's functions. For G B , at λ = 0, the solution is block diagonal: while for G A , the boundary condition gives off-diagonal components: Here the additional term imposes the boundary condition (14). We then need to compute (21) using these solutions of Green's functions. The important observation is that: 1. The Green's function G B is of the same form as two copies of G β , which is the Green's function on a thermal ensemble.
2. The Green's function G A is of the same form as that of the Kourkoulou-Maldacena states [67]. Consider pairing the Majorana fields as c j = , the Kourkoulou-Maldacena states is given by applying imaginary evolution to the eigenstate of n j = c † j c j : As a result, e −δS (2) corresponds to the difference between 〈K M |K M 〉 = 〈s|e −β H |s〉 and tre −β H with Majorana fermion number 2M . However, we know that for different s j , 〈K M |K M 〉 should be the same, and proportional to corresponding tre −β H with a possible factor that does not depend on the interaction strength [67]. We could just determine the factor by considering J → 0 4 . This leads to S (n) A = M log(2)/2 for small λ and any β J, consistent with analysis in the previous section.

Numerical results for the SYK-like models
Equations (19) and (20) give a complete answer to the subsystem entropy in the large-N limit. However, for general λ it is hard to study them analytically due to the lack of time translational invariance. Instead here we perform a numerical study. In this section, we first describe our numerical procedure, and then presents numerical results for different SYK-like models.
Numerically, we need to first solve (20) iteratively. We discretize the imaginary time τ ∈ [0, 2β) to L points, then both G s and Σ s (for s = A, B) become L × L matrices. We typically take L ∼ 200-300. After the discretization, the Schwinger-Dyson equation is written as: Here, the boundary condition (14) is taken into account by the bare Green's function G 0 s . Explicitly, we have [55]: and zero otherwise. Here sgn(x) is the sign function with sgn(0) = 0. Importantly, to obtain accurate numerical result with small number of points, we discretize G 0 s and take the inverse to get (G 0 s ) −1 , instead of using the discretized version of ∂ τ , which would lead to an oscillating Green's function.
After solving (20), we plug the results into (19) to compute the action. To ensure the convergence, we add to the action. An extrapolation to L → ∞ is performed finally by evaluating the action for different L.

The original SYK model
A /N for β J = 50 and q = 4 for different λ (blue circle). The green square is the exact diagonalization result [1]. The yellow triangle is thermal entropy with the effective temperature obtained by numerically matching the energy density (8). We also plot the analytical approximation (9) using a red line and λ log 2/2 using a purple line.
We firstly present our numerical results for the original SYK model. In Figure 1, we show a typical result for β J = 50 with q = 4 (blue circle). We take large β J so that it corresponds to extremely low temperature. The green squares are results of the ground state entanglement entropy from exact diagonalization with N = 44 Majorana fermions [1]. We also compare S (2) A with the thermal Rényi entropy for the canonical ensemble (yellow circle), where the effective temperature of the corresponding ensemble is determined by matching the energy density (8) numerically. The analytical approximation (9) is shown in red line.
We find that for λ 0.4, the numerical result for S A matches the exact diagonalization result well. When λ = 1/2, we find S A (1/2) = 0.334N /2. When the full system is prepared in a thermal ensemble, we find the subsystem entropy matches the thermal entropy well, except very close to λ = 1. The analytical formula (9) also gives a very good approximation, especially for small λ 0.4.
We then study the dependence of subsystem entropy on q and temperature T . We firstly consider the q dependence. Here we use the standard convention [40] and we fix the βJ = 50/ 2, where J is related to J by J 2 = 2 q−1 J 2 /q. This corresponds to β J = 50 for q = 4. According to (5), the thermal entropy even at zero temperature would approach N log 2/2 in the large-q limit.
Consequently, the full curve for S (2) A will approach S (2) A = M log 2/2 as q → ∞, as shown in Figure 2 (a). Similarly, as we increase the temperature, S A will approach the same upper bound (See Figure 2 (b)). For small q or low temperature, S A = M log 2/2 is non-monotonic. In particular, when q = 2, the zero temperature entropy density is zero and S (2) A in the limit β → ∞ is fully symmetric around λ = 1/2. The numerical result is shown in Figure 2 (a) and is consistent with an analytical result derived from the random matrix theory [1] (More details of this method is given in the appendix). Note (a) that for all cases, the slope for S (2) A at small λ collapses to λ log 2/2, as proved in the previous section.

< l a t e x i t s h a 1 _ b a s e 6 4 = " 7 B Y / x + V U n s 3 F Y o n I A 1 k 2 e 3 n g Y r c = " > A A A C x n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k W o m 5 K o o M u i m y 4 r 2 g f U I s l 0 W k P T J E w m S i m C P + B W P 0 3 8 A / 0 L 7 4 x T U I v o h C R n z r 3 n z N x 7 / S Q M U u k 4 r z l r b n 5 h c S m / X F h Z X V v f K G 5 u N d M 4 E 4 w 3 W B z G o u 1 7 K Q + D i D d k I E P e T g T 3 R n 7 I W / 7 w T M V b t 1 y k Q R x d y n H C u y N v E A X 9 g H m S q I u y t 3 9 d L D k V R y 9 7 F r g G l G B W P S 6 + 4 A o 9 x G D I M A J H B E k 4 h I e U n g 5 c O E i I 6 2 J C n C A U 6 D j H P Q q k z S i L U 4 Z H 7 J C + A 9 p 1 D B v R X n m m W s 3 o l J B e Q U o b e 6 S J K U 8 Q V q f Z O p 5 p Z 8 X + 5 j 3 R n u p u Y / r 7 x m t E r M Q N s X / p p p n / 1 a l a J P o 4 0 T U E V F O i G V U d M y 6 Z 7 o q 6 u f 2 l K k k O C X E K 9 y g u C D O t n P b Z 1 p p U 1 6 5 6 6 + n 4 m 8 5 U r N o z k 5 v h X d 2 S B u z + H O c s a B 5 U 3 M P K w f l R q X p q R p 3 H D n Z R p n k e o 4 o a 6 m i Q 9 w C P e M K z V b M i K 7 P u P l O t n N F s 4 9 u y H j 4 A Q m y P z w = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " b / c M j / 3 M a k p f v Q l 5 k y s P O U F a e E 8 = " > A A A C x n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k W o m 5 K o o M u i m y 4 r 2 g f U I s l 0 W o e m S U g m S i m C P + B W P 0 3 8 A / 0 L 7 4 x T U I v o h C R n z r 3 n z N x 7 / T g Q q X S c 1 5 w 1 N 7 + w u J R f L q y s r q 1 v F D e 3 m m m U J Y w 3 W B R E S d v 3 U h 6 I k D e k k A F v x w n 3 R n 7 A W / 7 w T M V b t z x J R R R e y n H M u y N v E I q + Y J 4 k 6 q L s 7 1 8 X S 0 7 F 0 c u e B a 4 B J Z h V j 4 o v u E I P E R g y j M A R Q h I O 4 C G l p w M X D m L i u p g Q l x A S O s 5 x j w J p M 8 r i l O E R O 6 T v g H Y d w 4 a 0 V 5 6 p V j M 6 J a A 3 I a W N P d J E l J c Q V q f Z O p 5 p Z 8 X + 5 j 3 R n u p u Y / r 7 x m t E r M Q N s X / p p p n / 1 a l a J P o 4 0 T U I q i n W j K q O G Z d M d 0 X d 3 P 5 S l S S H m D i F e x R P C D O t n P b Z 1 p p U 1 6 5 6 6 + n 4 m 8 5 U r N o z k 5 v h X d 2 S B u z + H O c s a B 5 U 3 M P K w f l R q X p q R p 3 H D n Z R p n k e o 4 o a 6 m i Q 9 w C P e M K z V b N C K 7 P u P l O t n N F s 4 9 u y H j 4 A R M 2 P 0
Before we end this subsection, we briefly discuss the connection between the entropy computed from the path integral method in the limit β → ∞ and the entanglement entropy of the ground state. The path integral formulation introduced in the previous section is a standard method to compute the entropy for a subsystem. In systems with local interaction, at finite β, the entropy for a subsystem is dominated by the thermal entropy and therefore is extensive. As we increase β to ∞, the thermal entropy density decreases to zero and the entropy is mainly coming from the entanglement entropy between the subsystem and its complement. The path integral formulation provides one method to compute the ground state entanglement entropy in the limit β → ∞. However, the above physics needs to be modified in SYK 4 model since it has extensive zero temperature entropy. Therefore the zero temperature entropy of the subsystem and the entanglement between it and its complement can have same order and both of them contribute to the entropy computed from the path integral formulation. Since the zero temperature entropy is proportional to the subsystem size, the entropy is not symmetric at λ = 1/2. Only when λ < 1/2 is it roughly equal to entanglement entropy of the ground state as presented in Fig. 1. Notice that this symmetry is fully restored in the SYK 2 model since the zero temperature entropy is zero.

SYK model with random hopping
In this subsection, we consider a specific generalization of the SYK model, which is the hybridization of SYK 2 and SYK q with q ≥ 4 [46][47][48]: Near the SYK q fixed point, the SYK 2 term is relevant, which governs the low energy physics as a Fermi liquid. Consequently, there is no extensive zero temperature entropy. Below we focus on the q = 4 case. For small V /J 1, the crossover from Fermi liquid (described by SYK 2 ) to the non-Fermi liquid (described by SYK 4 ) happens at a temperature T ∼ V 2 /J, where the thermodynamical entropy increases rapidly from S/N ∼ 0 to S/N ∼ O(1).
The path integral formulation for computing the entropy with this additional random hopping term is given by modifying the action (19) and thus (20). The additional term in the (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 B Y / x + V U n s 3 F Y o n I A 1 k 2 e 3 n g Y r c = " > A A A C x n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k W o m 5 K o o M u i m y 4 r 2 g f U I s l 0 W k P T J E w m S i m C P + B W P 0 3 8 A / 0 L 7 4 x T U I v o h C R n z r 3 n z N x 7 / S Q M U u k 4 r z l r b n 5 h c S m / X F h Z X V v f K G 5 u N d M 4 E 4 w 3 W B z G o u 1 7 K Q + D i D d k I E P e T g T 3 R n 7 I W / 7 w T M V b t 1 y k Q R x d y n H C u y N v E A X 9 g H m S q I u y t 3 9 d L D k V R y 9 7 F r g G l G B W P S 6 + 4 A o 9 x G D I M A J H B E k 4 h I e U n g 5 c O E i I 6 2 J C n C A U 6 D j H P Q q k z S i L U 4 Z H 7 J C + A 9 p 1 D B v R X n m m W s 3 o l J B e Q U o b e 6 S J K U 8 Q V q f Z O p 5 p Z 8 X + 5 j 3 R n u p u Y / r 7 x m t E r M Q N s X / p p p n / 1 a l a J P o 4 0 T U E V F O i G V U d M y 6 Z 7 o q 6 u f 2 l K k k O C X E K 9 y g u C D O t n P b Z 1 p p U 1 6 5 6 6 + n 4 m 8 5 U r N o z k 5 v h X d 2 S B u z + H O c s a B 5 U 3 M P K w f l R q X p q R p 3 H D n Z R p n k e o 4 o a 6 m i Q 9 w C P e M K z V b M i K 7 P u P l O t n N F s 4 9 u y H j 4 A Q m y P z w = = < / l a t e x i t > (b) action is:

< l a t e x i t s h a 1 _ b a s e 6 4 = " b / c M j / 3 M a k p f v Q l 5 k y s P O U F a e E 8 = " > A A A C x n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k W o m 5 K o o M u i m y 4 r 2 g f U I s l 0 W o e m S U g m S i m C P + B W P 0 3 8 A / 0 L 7 4 x T U I v o h C R n z r 3 n z N x 7 / T g Q q X S c 1 5 w 1 N 7 + w u J R f L q y s r q 1 v F D e 3 m m m U J Y w 3 W B R E S d v 3 U h 6 I k D e k k A F v x w n 3 R n 7 A W / 7 w T M V b t z x J R R R e y n H M u y N v E I q + Y J 4 k 6 q L s 7 1 8 X S 0 7 F 0 c u e B a 4 B J Z h V j 4 o v u E I P E R g y j M A R Q h I O 4 C G l p w M X D m L i u p g Q l x A S O s 5 x j w J p M 8 r i l O E R O 6 T v g H Y d w 4 a 0 V 5 6 p V j M 6 J a A 3 I a W N P d J E l J c Q V q f Z O p 5 p Z 8 X + 5 j 3 R n u p u Y / r 7 x m t E r M Q N s X / p p p n / 1 a l a J P o 4 0 T U I q i n W j K q O G Z d M d 0 X d 3 P 5 S l S S H m D i F e x R P C D O t n P b Z 1 p p U 1 6 5 6 6 + n 4 m 8 5 U r N o z k 5 v h X d 2 S B u z + H O c s a B 5 U 3 M P K w f l R q X p q R p 3 H D n Z
and consequently To study this crossover from the entropy perspective, we could either consider tunning V /J with fixed temperature T /J (Shown in Figure 3 (a)) or tunning T /J with fixed temperature V /J (Shown in Figure 3 (b)). We find that the slopes for S (2) A at small λ collapse to λ log 2/2, since the argument in the last section also works here. In the low-temperature Fermi liquid phase dominated by SYK 2 , the entropy of the full system is almost zero and S (2) A is almost symmetric around λ = 1/2, while for the non-Femri liquid phase the entropy of the total system is always large.

Complex SYK model with chemical potential
Here we consider another generalization with charge conservation called the complex SYK model. Instead of Majorana fermions, we consider complex fermions with annihilation (creation) operator c i (c † i ) satisfying {c i , c † j } = δ i j with i = 1, 2, ..., N . The Hamiltonian reads The total charge Q = i c † i c i is conserved [Q, H c ] = 0. One could consider an additional chemical potential term by adding H c − µQ, which changes the filling Q = 〈Q〉 /N . The Q dependence of the zero temperature entropy has been worked out in [44,45]. When the chemical potential becomes sufficiently large in the low temperature limit, a first order transition is found numerically, where charge Q jumps from finite to zero [68]. Here we again focus on the q = 4 case.
For the complex fermion case, we define the subsystem A by containing M complex fermion modes. With the chemical potential µ, the G − Σ action is now written as The saddle point equation becomes For the particle-hole symmetric case with µ = 0, this is proportional to the Majorana case (19) by a factor of 2, since in this case G s (τ, τ ) = −G s (τ , τ). For a small subsystem A, we expect it is maximally entangled under the density constrain due to its strong correlation with the rest of the system. This means that for a single qubit, in particle number basis, the density matrix should be ρ A = Q|1〉〈1|+(1−Q)|0〉〈0|. This leads to a Rényi entropy − log(Q 2 + (1 − Q) 2 ), smaller than log 2 away from half filling Q = 1/2. The numerical result is shown in Figure 4 (a). We find that when µ/J becomes larger, the slope of S In Figure 4 (b), we compare the fitting of initial slope to the prediction from − log(Q 2 + (1 − Q) 2 ). The charge Q is determined numerically using a thermal ensemble. We find the formula gives a good approximation. The difference should be attributed to the fact that the relation between µ and Q may receive corrections from the replicated system.

Conclusion
In this work, we use the path integral approach to analytically derive the self-consistent equation governing the entropy of a subsystem for a thermal ensemble of the SYK-like models in the large-N limit. We numerically solve the equation and study the scaling of subsystem Rényi entropy. In particular, we focus on the second Rényi entropy and we find results of SYK model at extremely low temperature are consistent with exact diagonalization [1] and can be well approximated by thermal entropy with an effective temperature [2] when λ ≤ 1/2. We also study the scaling of Rényi entropy as we vary the parameters q and temperature.
We further study the generalizations of the SYK model by introducing the quadratic random hopping term or considering the complex fermion version with U(1) symmetry. For the random hopping case, we see the crossover from the low-energy Fermi liquid to the non-Fermi phase at intermediate temperature. For the complex SYK case, we observe the chemical potential dependence of entropy for λ → 0. We also find that a first order transition exists when tuning the subsystem size at large chemical potential. j = n − 1 − j pair up to give real values. We see that the general form is complicated, and the way to analytic continue to other α values is not obvious.
Another interesting case is the Rényi entropy at n → +∞. Using The value S