title: Non-unitary dynamics of Sachdev-Ye-Kitaev chain
authors: Chunxiao Liu Pengfei Zhang Xiao Chen
text:  
  
  
  
  
  Department of Physics, University of California Santa Barbara, Santa Barbara, California, 93106, USA
  
  
  PengfeiZhang.physics@gmail.com
  Institute for Quantum Information and Matter and Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, California 91125, USA
  
  
  chenaad@bc.edu
  Department of Physics, Boston College, Chestnut Hill, MA 02467, USA
  
  
  
  
  
  
      We construct a series of one-dimensional non-unitary dynamics consisting of both unitary and imaginary evolutions based on the  Sachdev-Ye-Kitaev model. Starting from a short-range entangled state, we analyze the entanglement dynamics using the path integral formalism in the large N limit. Among all the results that we obtain, two of them are particularly interesting: (1) By varying the strength of the imaginary evolution, the interacting model exhibits a first order phase transition from the highly entangled volume law phase to an area law phase; (2) The one-dimensional free fermion model displays an extensive critical regime with emergent two-dimensional conformal symmetry.  
  
  
  Non-unitary dynamics of Sachdev-Ye-Kitaev chain
      Xiao Chen
      December 30, 2023
  ===============================================
  
  
  
  
  
  
  
  § INTRODUCTION
  
  Recent years have witnessed tremendous breakthrough in many-body quantum dynamics. For a closed many-body quantum system decoupled from the environment, under the unitary dynamics, the interaction in the system can lead to chaos and thermalize all the small subsystems. The total wave function acts as its own heat bath and this phenomenon is referred as quantum thermalization <cit.>. 
  
  The irreversible thermalization process can be avoided if we allow non-unitary evolution, which naturally arises in open quantum systems.  Recently it is observed that a unitary dynamics subjected to repeated measurement can exhibit non-thermal phases if we follow the quantum trajectory of the many-body wave function. More strikingly, by varying the measurement rate, there is a continuous entanglement phase transition <cit.>. In the phase with slow measurement rate, the state remains highly entangled and the entanglement entropy obeys volume law scaling, while in the phase with fast measurement rate, the entanglement entropy obeys area law scaling. Notice that the explicit form of the measurement is not important and can be either a projective measurement or a more generalized weak measurement <cit.>. 
  
  Motivated by these findings, we consider the following question: For a non-unitary dynamics U∼exp(-i H t) governed by a non-Hermitian Hamiltonian
  
      H=H_1-ig H_2,
  
   can we realize an entanglement phase transition by varying g? In the above equation, both H_1 and H_2 are Hermitian Hamiltonians. More specifically, in this paper, we consider the following interaction: H_1 is a Hamiltonian describing interaction between different sites and H_2 is a Hamiltonian defined at each site.  H_2 for example can describe the coupling of the onsite degrees of freedom to an external field. For such non-unitary dynamics, in the limit g=0, we expect that the steady state will typically saturate to a highly entangled state with a volume law scaling, while in the limit g→∞, this becomes a purely imaginary evolution and the steady state is a trivial product state with zero entanglement entropy. 
  In a strongly interacting system, it is not obvious if there is a phase transition occurring at finite g. To address the above questions, we consider a one-dimensional (1D) non-unitary dynamics constructed from Sachdev-Ye-Kitaev (SYK) model <cit.> and explore the possible phase transition in it. 
  
  The SYK model is a fermionic system with random all-to-all interaction <cit.>. This model can be analytically solved in the large N limit. Due to this special property, this model has vast applications in different fields including high energy, condensed matter physics and quantum information theory.  Many variants of the SYK model have been constructed to study quantum chaos, quantum gravity and non-Fermi liquid analytically <cit.>. In particular, there are studies on the entanglement entropy of the SYK model <cit.>, where transitions to the replica wormhole <cit.> solution are found in the long time limit.
  
  In this paper, we will use the SYK model to construct a set of 1D chain models and explore the non-unitary dynamics harbored in them. We study the entanglement dynamics by using the path integral formalism, which can be obtained self-consistently by virtue of the large N nature of the SYK model. By varying g, we observe different entanglement scaling behaviors which correspond to different saddle point solutions.  We hope these results in the large N limit could shed light on the more generalized phase transition in interacting systems at finite N where analytical tools are lacking.    
  
  The rest of the paper is organized as follows. In Sec. <ref> we define the SYK chain models, derive the path integral formalism for the Rényi entropy of these models, and write down the corresponding saddle point equations which are amenable to numerical study. Then in Sec. <ref> we apply the formalism to two models: the interacting model with inter-cluster SYK_4 coupling and intra-cluster SYK_2 coupling, and the non-interacting model with inter-cluster SYK_2 coupling and intra-cluster SYK_2 coupling. We find that the steady state of the interacting model exhibits either a volume-law or an area-law phase as the coupling λ is varied, and the two phases are separated by a first order transition, while that of the non-interacting model exhibits critical behavior for any finite coupling λ. We finally summarize and discuss these results in Sec. <ref>. Some derivation details, as well as results for another model with both inter- and intra-cluster SYK_4 interactions, are given in the Appendices.
  
  
  
  
  
  § MODEL AND METHOD
  
  
  We consider non-unitary time evolution of SYK chain generated by the following non-Hermitian Hamiltonian:
  
      H ≡ J H^p_inter-cluster - i V H^q_intra-cluster≡ J  ∑_x H^p_x,x+1 - i V  ∑_x H^q_x,
  
  where H^p_inter-cluster is the inter-site interaction with coupling strength J and H^q_intra-cluster is the onsite interaction with coupling strength V. From now on, we drop their subscript for conciseness. At each site, there are N Majorana fermions. H^p_x,x+1 denotes p-body random interaction between neighboring sites x and x+1 and H^q_x is q-body random interaction at site x (See the cartoon in Fig. <ref>). Their explicit forms are
  
      H_x,x+1^p= i^p/2∑_i_1<i_2<⋯<i_p/2,j_1<j_2<⋯<j_p/2 J^x,x+1_i_1i_2⋯ i_p/2j_1j_2⋯ j_p/2χ^x_i_1χ^x_i_2⋯χ^x_i_p/2χ^x+1_j_1χ^x+1_j_2⋯χ^x+1_j_p/2, 
          H_x^q= i^q/2∑_i_1<i_2<⋯<i_q V^x_i_1i_2⋯ i_qχ^x_i_1χ^x_i_2⋯χ^x_i_q.
  
  where the subscript i_m,j_n=1,2,…,N is the Majorana flavor index at each site and the superscript x=1,2,...,L is the site of the SYK chain; we assume periodic boundary condition for the chain so that L+1≡ 1. J^x,x+1_i_1⋯ i_p/2j_1⋯ j_p/2 and V^x_i_1i_2⋯ i_q
  are time-independent Gaussian random variables with vanishing mean and the following variance
  
      |J^x,x+1_i_1⋯ i_p/2j_1⋯ j_p/2|^2 = (p/2)!(p/2-1)!/2N^p-1,        |V^x_i_1i_2⋯ i_q|^2 = (q-1)!/N^q-1.
  
  Note that both H^p_x,x+1 and H^q_x are Hermitian Hamiltonian. Therefore in the time evolution governed by exp(-i Ht), H^p_x,x+1  provides the unitary evolution while H^q_x is responsible for the imaginary evolution. We describe the path integral formalism for this non-unitary evolution in the next subsection. 
  
  
  
  
  
  
   §.§ Rényi entropy under non-unitary evolution 
  
  We are interested in the entanglement dynamics of certain pure state. We prepare an initial state |ψ(t=0)⟩ = ⊗_x=1^L |{0}⟩_x which is the state annihilated by all the complex fermions defined by c^x_j = χ^x_2j-1+i χ^x_2j as c^x_j|{0}⟩_x = 0, for all x ∈ 1,2,...,L and j= 1,2,...,N. Since our initial state is a product state between different sites, there is no spatial entanglement. 
  
  We then evolve |ψ(t=0)⟩ under the Hamiltonian (<ref>) to entangle different sites. At time t, the state reads
  
      |ψ(t)⟩ =  1/√(Z[{J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1i_2...i_q}](t))e^-i t H|{0}⟩
      =1/√(Z[{J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1i_2...i_q}](t)) e^-it J H^p - t V H^q |{0}⟩,
  
  where 
  
      Z[{J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1i_2...i_q}](t) = ⟨{0}| e^i t H^† e^-i t H|{0}⟩
  
  is introduced to ensure normalization of the state |ψ(t)⟩. From now on, we keep the dependence of random parameters implicit. It is useful to write √(Z(t))|ψ(t)⟩ and Z(t) as path-integrals:
  
      ∫_b.c. Dχ^x_i(τ) exp(-∫ dτ{1/2∑_x,iχ^x_i∂_τχ^x_i+∑_xVH^q_x[χ^x_i]+f(τ)∑_xJH^p_x,x+1[χ^x_i]})
  
  Here we introduce f(τ) to specify whether the real-time evolution is forward or backward: f(τ)=i for forward evolutions e^-iH^p t and f(τ)=-i for backward evolutions e^iH^p t [Note that the direction of real-time evolution is contour-dependent, i.e f(τ) is a contour-dependent quantity.]. We use b.c. to denote proper boundary conditions specified by the pictorial representations following Ref. <cit.>:
  
      √(Z(t))|ψ(t)⟩=
      
      {[thick,scale = 0.75,baseline=([yshift=-2.7pt]current bounding box.center)]
      
        (1.0,0.5) arc(90:-90:0.5 and 0.5);
      
       [mid arrow]  (1.,0.5)  –  (-1.0,0.5);
       [mid arrow]  (1.,-0.5)  –  (-1.0,-0.5);
      
       [dotted]  (-0.1,-0.5)  –  (-0.1,0.5);
       [dotted]  (-0.6,-0.5)  –  (-0.6,0.5);
       [dotted]  (0.4,-0.5)  –  (0.4,0.5);
       [dotted]  (0.9,-0.5)  –  (0.9,0.5);
      
         (1.5,0) circle (1.5pt) node[left];
      
        (2,0.) node|{0}⟩;
        (-0.9,0.8) nodet;
        (1.4,0.8) node0;
      
        (0.1,1) nodeχ^x_e;
        (0.1,-1) nodeχ^x_o;
      },      
      Z(t)=
      
      {[thick,scale = 0.75,baseline=([yshift=-3.1pt]current bounding box.center)]
      
        (1.0,0.5) arc(90:-90:0.5 and 0.5);
        (-1.0,0.5) arc(90:270:0.5 and 0.5);
      
       [mid arrow]  (1.,0.5)  –  (.0,0.5);
       [mid arrow]  (1.,-0.5)  –  (.0,-0.5);
      
        [mid arrow]  (-1.,0.5)  –  (.0,0.5);
       [mid arrow]  (-1.,-0.5)  –  (.0,-0.5);
      
       [dotted]  (-0.0,-0.5)  –  (-0.,0.5);
       [dotted]  (-0.5,-0.5)  –  (-0.5,0.5);
       [dotted]  (0.5,-0.5)  –  (0.5,0.5);
       [dotted]  (1,-0.5)  –  (1,0.5);
       [dotted]  (-1,-0.5)  –  (-1,0.5);
      
         (1.5,0) circle (1.5pt) node[left];
          (-1.5,0) circle (1.5pt) node[left];
      
        (2,0.) node|{0}⟩;
         (-2,0.) node⟨{0}|;
        (0,0.8) nodet;
        (1.4,0.8) node0;
        (-1.4,0.8) node0;
      
        (0.1,1.1) nodeχ^x_e;
        (0.1,-1.1) nodeχ^x_o;
      }
  
  where we have separated out the Majorana modes with even/odd indices as χ^x_e/o and x∈[1,L]: they represents Majorana fermions χ^x_2j/χ^x_2j-1. The solid lines represent the evolution and the black dots denote the initial state at t=0, at which the even and odd Majoranas are related due to c^x_j|{0}⟩_x=0. The dotted lines represent interaction between fermions, which contains coupling between different sites. For √(Z(t))|ψ(t)⟩, additional quantum state is attached to the free ends.
  
  We are interested in the Rényi entropy of |ψ(t)⟩. We divide the SYK chain into subsystems A and B, with A containing sites x=1,2,...,L_A. The reduced density matrix ρ_A = Tr_B ρ(t) is obtained by tracing over the degrees of freedom in B. A pictorial representation of ρ_A reads:
  
      Z(t)ρ_A(t)=
      {[thick,scale = 0.75,baseline=([yshift=-4pt]current bounding box.center)]
      
       [red] (1.0,0.5) arc(90:-90:0.5 and 0.5);
       [blue] (1.0,0.3) arc(90:-90:0.5 and 0.5);
      
       [mid arrow,red]  (1.,0.5)  –  (-1.0,0.5);
       [mid arrow,red]  (1.,-0.5)  –  (-1.0,-0.5);
       [mid arrow,blue]  (1.,0.3)  –  (-1.0,0.3);
       [mid arrow,blue]  (1.,-0.7)  –  (-1.0,-0.7);
      
      
       [dotted]  (-0.1,-0.7)  –  (-0.1,0.5);
       [dotted]  (-0.6,-0.7)  –  (-0.6,0.5);
       [dotted]  (0.4,-0.7)  –  (0.4,0.5);
       [dotted]  (0.9,-0.7)  –  (0.9,0.5);
      
         (1.5,0) circle (1.5pt) node[left];
         (1.5,-0.2) circle (1.5pt) node[left];
      
        (2,-0.1) node|{0}⟩;
        (-0.9,0.8) nodet;
        (1.4,0.8) node0;
      
        (0.1,1) nodeχ^x_e;
        (0.1,-1) nodeχ^x_o;
      
      
       [red] (4.0,0.5) arc(90:270:0.5 and 0.5);
       [blue] (4.0,0.3) arc(90:270:0.5 and 0.5);
      
       [mid arrow,red]  (4.,0.5)  –  (6.0,0.5);
       [mid arrow,red]  (4.,-0.5)  –  (6.0,-0.5);
       [mid arrow,blue]  (4.,0.3)  –  (6.0,0.3);
       [mid arrow,blue]  (4.,-0.7)  –  (6.0,-0.7);
      
      
       [dotted]  (5.1,-0.7)  –  (5.1,0.5);
       [dotted]  (5.6,-0.7)  –  (5.6,0.5);
       [dotted]  (4.6,-0.7)  –  (4.6,0.5);
       [dotted]  (4.1,-0.7)  –  (4.1,0.5);
      
         (3.5,0) circle (1.5pt) node[left];
         (3.5,-0.2) circle (1.5pt) node[left];
      
        (3,-0.1) node⟨{0}|;
        (5.9,0.8) nodet;
        (3.6,0.8) node0;
      
        (4.9,1) nodeχ^x_e;
        (4.9,-1) nodeχ^x_o;
       
       [blue,dashed] (-1,0.3) arc(180:0: 3.5 and 0.5);
       [blue,dashed] (-1,-0.7) arc(180:360: 3.5 and 0.5);
      }
  
  Here the red line represents the A subsystem and the blue line represents the B subsystem. The nth Rényi entropy is then defined as 
  
      S^(n) = - logTr_A ρ_A^n/n-1=-log(Z_n(t,L_A)/Z^n(t))/n-1,
   
  where Z_n(t,L_A) is given by sewing n copies of reduced density matrix ρ_A, as shown in Fig. <ref>, where we have labeled contour 𝒞 anti-clockwise by a single time variable s ∈ [0,4nt]; note also that the contour sews even and odd Majorana fields χ^x_2j/χ^x_2j-1 (See Appendix <ref> for details). The contour itself again defines the path-integral for computing the Rényi entropy, with boundary conditions indicated by dashed lines and black dots. From the contour, it is obvious that there is a symmetry by interchanging A and B: S^(n)(t,L_A)=S^(n)(t,L-L_A). In addition, we have Z_n(t,0)=Z_n(t,L)=Z^n(t), which is consistent with the fact that the full system is in a pure state.
  
  To proceed, we need to average over random variables J^x,x+1_i_1...i_p/2j_1...j_p/2 and V^x_i_1...i_q  after computing logTr_A ρ_A^n for each realization of J^x,x+1_i_1...i_p/2j_1...j_p/2 and V^x_i_1...i_q. Additional disorder replicas are required to accomplish this. In the SYK models that we consider, however, the computation can be simplified by working within a disorder replica diagonal ansatz: both numerical analysis and analytical arguments <cit.> suggest that at leading order of 1/N, the Rényi entropy can be approximated by 
  
      S^(n) = - logTr_A ρ_A^n/n-1≃  - logTr_A ρ_A^n/n-1 = -1/n-1(logZ_n(t,L_A) - logZ_n(t,0)),
  
  and the result becomes exact in the large N limit. 
  
  
  
  After disorder average, to compute Z_n(t) one could further introduce the bilocal fields G_x and Σ_x <cit.>. Leaving details for Appendix <ref>, we find
  
      Z_n(t)  = ∫_𝒞D G_xD Σ_x e^-NS_n[ G_x, Σ_x],
  
  with
  
      S_n[ G_x, Σ_x] =- N/4∑_x{log[ ∂_s,x -  Σ_x ] - ∫^4nt_0ds ∫^4nt_0 ds' [  Σ_x  G_x -  f(s)f(s')J^2/p( G_x  G_x+1)^p/2P -V^2/qG_x^qP]}.
  
  As aforementioned, the time parameters s,s' ∈ [0,4nt] are parameterized on the contour 𝒞. P(s,s') is defined as
  
      P(s,s') =  {[                                                                             1,                                  s,s'  both parameterizing even or odd fields,;                                                                             0, s  parameterizing odd (even) fields and  s'  parameterizing even (odd) fields. ].
  
  
  In the large N limit, since the action (<ref>) is proportional to N, we obtain the saddle point solution
  
  
      G_x(s,s')   =   (  ∂_s,x-  Σ_x^J -  Σ_x^V)^-1(s,s'),
      Σ^J_x(s,s')    =   J^2/2 f(s)f(s')  G_x^p/2-1(s,s')(  G_x+1^p/2(s,s')+ G_x-1^p/2(s,s'))P(s,s'),
      Σ^V_x(s,s')    =    V^2  G_x^q-1(s,s')P(s,s').
  
   
  Therefore the Rényi entropy is
  
      S^(n) = 1/n-1(S_n,saddle(L_A) - S_n,saddle(0)),
  
  where S_n,saddle is the on-shell action (<ref>) obtained from the solution of the Eqs. (<ref>). Explicitly, we have
  
      S_n,saddle[ G_x, Σ_x] =- N/4∑_x{log[ ∂_s,x -  Σ_x ] - ∫^4nt_0ds ∫^4nt_0 ds' [ p-1/pΣ_x^J  G_x+q-1/qΣ_x^V  G_x ]}.
  
  
  
  
  
   §.§ Brownian dynamics
  
  We can also generalize the above formalism to a non-unitary Brownian SYK model in which J^x,x+1 and/or V^x are/is independent Gaussian random variables in time, with vanishing mean and the following variance
  
      J^x,x+1_i_1...i_p/2j_1...j_p/2(τ) J^x,x+1_i_1...i_p/2j_1...j_p/2(τ') = δ(τ-τ')  (p/2)!(p/2-1)!/2JN^p-1,    V^x_i_1...i_q(τ)V^x_i_1...i_q(τ') = δ(τ-τ')  (q-1)!/VN^q-1,
  
  i.e. the correlation is non-vanishing only at equal-time. A variant of this model was studied in Ref. <cit.>.
  Previously, there has been extensive study on the Brownian unitary circuit models <cit.>. These models can provide analytical solutions for the quantum dynamics even with a small onsite Hilbert space and can sometimes show different dynamics from systems without randomness <cit.>. Here we plan to explore the influence of the temporal correlation on the non-unitary quantum dynamics.
  
  One can apply the entire procedure of the last subsection to computing the Rényi entropy of the Brownian SYK model with minimal change. Note that in the Brownian model, a time-dependent random variable admits only equal-time correlation t=t', which happens at the redefined times s' = 2jt + s and s'=2(j+1)t - s on the contour for Z_n(t). Therefore we need to introduce an additional projection to the Brownian term
  
      P_B(s,s') =  ∑_j=1^2nδ(s-s'- 2j t) + δ(s+s' - 2jt).
  
  As an example, when the inter-site interaction becomes Brownian, the saddle point equation (<ref>) becomes
  
  
      Σ^J_x(s,s')    =   J/2 f(s)f(s')  G_x^p/2-1(s,s')(  G_x+1^p/2(s,s')+ G_x-1^p/2(s,s'))P(s,s')P_B(s,s').
  
   
   Similar modification applies when V^x becomes Brownian. The on-shell action can still be computed as in (<ref>). 
  
  
  
   §.§ Numerical details
  
  
  The complex form of the saddle point equations (<ref>) hampers analytical treatment for generic parameters of J, V, t, and  L_A, and numerical methods are exploited instead. We will hereby focus on the second Rényi entropy with n=2. We discretize s ∈ [0,8t) to 4𝖫 points, and represent G_x and Σ_x as 4𝖫× 4𝖫 matrices. Numerically, for moderate J and V, we take 𝖫= O(1)× t and finite size scaling for 𝖫 is performed as the final procedure to obtain the 𝖫→∞ result. 
  
  One solves the numerical version of Eqs. (<ref>) iteratively: at the nth step, one plugs the Green's function obtained from the previous step, G^(n-1)_x, into the numerical form of Eqs. (<ref>) and (<ref>), then use Eq. (<ref>) to generate new Green's function G^(n-1)_x,new, which is to be used at the (n+1)th step according to G^(n)_x=(1-ρ) G^(n-1)_x+ρ G^(n-1)_x,new, where 0<ρ<1 is some weight that controls the rate of convergence. The saddle point solution is found when the Green's function series {G^(n)_x} converge within numerical precision. The numerical version of (<ref>) is
  
  
      (G_x)_ij   =   [(G^0_x)^-1- Σ_x^J - Σ_x^V]^-1_ij,
      (Σ^J_x)_ij   =   J^2/2 f_if_j (G_x)^p/2-1_ij( (G_x+1)_ij^p/2+(G_x-1)_ij^p/2)P_ij(Δ t)^2,
      (Σ^V_x)_ij   =    V^2 (G_x)_ij^q-1P_ij(Δ t)^2,
  
   
  where Δ t=2t/𝖫. P and f are the straight-forward discretization of their definitions (see Eq. (<ref>) and the paragraph below Eq. (<ref>)):
  
      P   =   [baseline=4.5ex,x=0.2cm,y=0.2cm] (0,0) – (8,0) – (8,8) – (0,8) – (0,0);
              [gray] (0,0) rectangle (1,1);
              [gray] (0,7) rectangle (1,8);
              [gray] (7,0) rectangle (8,1);
              [gray] (7,7) rectangle (8,8);
              [gray] (1,1) rectangle (3,3);
              [gray] (1,5) rectangle (3,7);
              [gray] (5,1) rectangle (7,3);
              [gray] (5,5) rectangle (7,7);
              [gray] (3,3) rectangle (5,5);
              [gray] (0,3) rectangle (1,5);
              [gray] (3,0) rectangle (5,1);
              [gray] (7,3) rectangle (8,5);
              [gray] (3,7) rectangle (5,8);  ,       f_j = {[                                             i, j ∈(𝖫/2, 𝖫)∪(3𝖫/2, 2𝖫) ∪(5𝖫/2, 3𝖫)∪(7𝖫/2, 4𝖫),;                                            -i,                                     otherwise. ].
   
  where we have pictorially presented the 4𝖫× 4𝖫 matrix of P, with the gray (white) representing unit (vanishing) entries . The partial operator ∂_τ in the continuum version (<ref>) is now discretized in Eqs. (<ref>) to the non-interacting Green's function (G^0_x)^-1
  
      x ∈ A:     (G^0_x)_ij≡(G^0_A)_ij= 1/2sgn(i-j)  for i,j ∈{1,2,...,2𝖫} or ∈{2𝖫+1,...,4𝖫},
      
      x ∈ B:     (G^0_x)_ij≡(G^0_B)_ij= (G^0_A)_i-𝖫,j-𝖫,
  
  where sgn(x) is the sign function with sgn(0)=0, the row and column indices are understood in the modulo sense i+4𝖫≡ i, and the unmentioned entries are identically zero. Note this substitution ∂_τ→ (G^0_x)^-1 is crucial to ensure numerically accurate result with small number of points. We introduce the following terminology: a solution is replica diagonal if for all x, the Green's function G^0_x is non-zero only for elements that are non-zero in G^0_x[The replica diagonal solution here refers to diagonality with respect to the Rényi index n. This is different from the diagonality in the disorder replica space mentioned in Subsection <ref>.], and a solution is replica quasi-diagonal if this property is satisfied for all x in the subsystem bulks (i.e. interior) A, B and is violated on the subsystem boundaries ∂ A and ∂ B.
  
  We consider two different types of initial conditions for numerical iteration (not to be confused with the Green's function solution at t=0):
  
      Type 1(A).    G_x^(0)=G^0_A   for x∈[1,𝖫],
         Type 1(B).   G_x^(0)=G^0_B   for x∈[1,𝖫], 
         Type 2(D).    G_x^(0)=G^0_A   for x∈ A   and   G_x^(0)=G^0_B   for x∈ B,
         Type 2(QA).    G_x^(0)=G^0_A   for x∈A,   G_x^(0)=G^0_B   for x∈ B and   
          G^(0)_x = 0.9G^0_A+0.1 G^0_B    for x∈∂ A,
         Type 2(QB).    G_x^(0)=G^0_A   for x∈  A,   G_x^(0)=G^0_B   for x∈B and   
          G^(0)_x = 0.9G^0_B+0.1 G^0_A    for x∈∂ B,
  
  By definition, the initial condition Type 2(D) is replica diagonal, and the initial condition Type 2(QA/QB) is replica quasi-diagonal. When q≥ 4, the replica diagonal/quasi-diagonal property is preserved under the iteration and the initial condition Type 2(D) (Type 2(QA/QB)) always leads to replica diagonal (quasi-diagonal) self-consistent solution. We have also checked that other initial conditions would not lead to new saddle points. 
  
  After solving Eqs. (<ref>), we compute S_2,saddle(L_A) as a function of the subsystem size L_A using the discretized version of Eq. (<ref>). The final result is obtained by using Eq. (<ref>)
  
      S^(2)/N = S_n=2,saddle(L_A) - S_n=2,saddle(0).
  
  
  
  
  
  
  
  
  
  § PHYSICAL PROPERTIES OF THE MODEL 
  
  
   In this section we present numerical results along with analytical understandings  for our model with different (p,q). We focus on models with (p,q)=(4,2) and (2,2) in the next two subsections. For the (4,2) model, we find that the steady state is a volume-law entanglement state for V≪ J and is an area-law entanglement state for J ≪ V. These two phases are separated by a first order transition. On the other hand, for the (2,2) model, we find the 
  steady state is always in a critical phase for finite V/J. We will briefly mention properties of the (p,q)=(4,4) model in Appendix <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Inter-cluster SYK_4  with intra-cluster SYK_2
  
  
  
  The (p,q)=(4,2) model consists of Hermitian, inter-cluster SYK_4 interaction with strength J and non-Hermitian, intra-cluster SYK_2 interaction with strength V.  We prepare a product state and let it evolve under this Hamiltonian described by Eq. (<ref>). We study the scaling of the second Rényi entropy S^(2) of the state as a function of J, V, subsystem size L_A, and evolution time t.
  
  
  
    §.§.§ V≪ J
  
  
  
  
  
  
  
  
  Before discussing the non-unitary dynamics, we first consider the unitary dynamics at V=0. Since the system has local interaction, we expect the entanglement entropy for a subsystem grows linearly in time and saturates to a constant proportional to L_A at late time. We compute S^(2) by numerically solving the saddle point solution and present the results in in Fig. <ref>. In summary, the entanglement dynamics satisfies
  
      S^(2)(t,L_A)={[               v_E t  if   t< s_0L_A/v_E;             s_0 L_A if    t≥ s_0L_A/v_E ].
  
  where the entanglement velocity v_E depends on the magnitude of J and the transition occurs at t=s_0L_A/v_E.
  Since the initial state we consider has energy expectation value ⟨ E ⟩=0, the final steady state is indeed a maximally entangled thermal state at infinite temperature with s_0=log 2/2. Notice that the entanglement dynamics in the large N model is distinct from that in a chaotic system with finite or small N, in which the energy conservation law leads to a diffusive dynamics for S^(2)(t), i.e., S^(2)(t)∼√(t) <cit.>. This slow dynamics disappears in the SYK chain model in the large N limit. 
  
  The above behavior still holds when V takes a small finite value, although with s_0< log2/2. For instance, we present the V=0.3 result in Fig. <ref>, which is analogous to the V=0 result. This indicates that when V≪ J, there exists a phase in which the steady state entanglement entropy has volume law scaling.  
  
  Here we explain this entanglement scaling from the saddle point solution perspective. The S^(2) in Fig. <ref> and Fig. <ref> is determined by the minimum value of the various saddle point solutions which are obtained by choosing the initial conditions Type 1(A/B) and Type 2(D/QA/QB) as mentioned in Sec. <ref>. As shown in Fig. <ref>, starting from an initial condition of Type 1, S^(2) quickly saturates to a volume law scaling, including a branch with S^(2)∝ L_A (Type 1(A)) and a branch with S^(2)∝ L_B (Type 1(B)). On the other hand, initial condition Type 2 always leads to a constant S^(2) that is insensitive to subsystem size. Furthermore, for small V≲ 0.3, we observe that the replica quasi-diagonal saddle obtained using Type 2(QA/QB) and the replica diagonal saddle obtained using Type 2(D) give slightly different results. Both solutions give S^(2) that grows linearly in time, however the replica quasi-diagonal saddle has a lower S^(2) value (see Fig. <ref>). Interestingly, as V increases, the replica quasi-diagonal saddle gradually merges with the replica diagonal saddle, and as V>0.3, both initial conditions Type 2(Q) and Type 2(D) give the same replica diagonal solution.  We will henceforth refer to the saddle point solutions by the initial conditions that lead to them.
  
  The above analysis suggests that early time S^(2) is determined by the Type 2 saddle which grows linearly in time, and that late time S^(2) is determined by the Type 1 saddle which eventually saturates to a constant that is proportional to L_A. The time dynamics is presented in Fig. <ref> and Fig. <ref> for V=0 and V=0.3 respectively and is consistent with the result in Eq. (<ref>). In addition, in Fig. <ref>, we verify that no finite size effect exists in both saddle solutions.  
  
  We can further understand the Type 1 and Type 2 saddles by analyzing the Green's functions. Here we take V=0.2 as an example. The S^(2) computed from different saddles are plotted in Fig. <ref>, and the subsystem bulk Green's functions for different saddles are shown in Figs. <ref>–<ref>. The Green's function corresponding to Type 2 saddles has a replica diagonal form (Fig. <ref>) inside the subsystem bulk. This means there exists self-consistent replica diagonal/quasi-diagonal solution of the Green's functions when computing Z_n(t). Compared to the Green's functions used for computing the normalization Z(t) (with proper time labeling), the Green's functions for Z_n(t) are modified only in the neighborhood of subsystem boundaries, since only near the boundaries is there a change of the Schwinger-Dyson equation. Consequently, terms in Z_n(t) and nZ(t) that involve only Green's functions in the subsystem bulk A and B cancel, and the Rényi entropy receives contribution only from the Green's functions on the subsystem boundaries ∂ A and ∂ B. This leads to a S^(2) which is independent of the subsystem size. On the contrary, if we take a replica non-diagonal saddle as in Fig. <ref> and Fig. <ref>, the Green's functions become highly replica non-diagonal in one of the subsystems (A or B, depending on the initial conditions being Type 1(A) or Type 1(B)). Consequently, there is non-vanishing contribution to the entropy even deeply inside the subsystem bulk. This may 
  give rise to a S^(2) proportional to the subsystem size. We mention in passing that the Green's functions also distinguish replica diagonal and replica quasi-diagonal saddles due to their difference at the subsystem boundaries. As an example, we show in Fig. <ref> the Green's functions of the replica quasi-diagonal saddle on subsystem boundaries for various V. The corresponding Green's functions in the subsystem bulk (not plotted) all have the same profile as in Fig. <ref>. The boundary Green's function profiles confirm that as V increases, the replica quasi-diagonal solution converges to the replica diagonal solution.
  
  Interestingly, the  simultaneous existence of such replica non-diagonal and replica diagonal/quasi-diagonal solutions is directly related to the recent resolution of the black hole information paradox. In the language of gravity theories, such highly replica non-diagonal solution is known as a “replica wormhole” <cit.>. In the story of the information paradox, without such replica wormhole, the entropy of the system would grow unboundedly, which ultimately violates unitarity. We find similar behavior exists in our model: without finding the Type 1 saddles, we would conclude that the entropy of our system grows linearly, finally exceeding the maximal entropy L_Alog 2/2. When the “replica wormhole” saddle point is taken into account, the entropy becomes well-behaved.
  
  
  
    §.§.§ V≫ J
  
  
  Now we turn to the other limit J/V→ 0. First, the J=0 case defines purely imaginary evolution governed by the intra-cluster SYK_2 interaction. The system evolves into its ground state after long time evolution. Since the SYK_2 model is a free fermion model and has vanishing zero temperature entropy <cit.>, the ground state is a simple product state, i.e.,
  
      |Ψ_ GS⟩=|ψ_ GS^1⟩|ψ_ GS^2⟩⋯ |ψ_ GS^L⟩
  
  where |ψ_ GS^x⟩ denotes the SYK_2 ground state at site x and there is no entanglement between different sites. When J≠ 0, we expect that the steady state entropy S^(2) still satisfies area law and confirm this result in Fig. <ref>. This scaling behavior can be understood as follows: As long as J≪ V, we can perturbatively study S^(2) by taking the Type 2(D) saddle solution. The leading behavior of S^(2) receives corrections only from the subsystem boundaries and can be written in closed form (see Appendix <ref> for derivation):
  
      S^(2) = δ S
      =  4J^2/p∫^t_0 ds ∫^2t_t ds' f(s)f(s')
      G_c^p
  
  where G_c is the conformally invariant Green's function
  
      G_c(s,s') =  1/2t V sinπ(s-s')/2t,
  
  for the ground state of SYK_2 valid for s-s' ≫ϵ≈ V^-1[This means the outer integral of Eq. (<ref>) has to be regulated as ∫^t-ϵ_ϵ ds.]. 
  Plugging G_c into Eq. (<ref>), we get
  
  
      S^(2)
      = a J^2/V^2 +  b log(  V t) /V^2 t^2
      + c J^2/V^4 t^2+  1/t^2O((ϵ/t)^2),   t≫ V^-1,
  
  where the prefactor of the logarithmic term b = J^2/3π^2 V^2 is cutoff-independent, while a, c are cutoff-dependent finite constant. The first term of  S^(2) is responsible for the area law scaling behavior of the steady state and the J^2/V^2 scaling form is confirmed numerically in Fig. <ref>.  The second term contains the leading time dependence at small J/V which shows a decaying behavior that is linear in log t/t^2, see the right inset of Fig. <ref>. We also verify the explicit form of the coefficient b, see the left inset of Fig. <ref>.
  
  
  
  
    §.§.§ The transition at the intermediate regime
  
  The above analysis suggests that there exists a volume law phase when V≪ J and an area law phase when V≫ J. We further explore the intermediate regime and find that there is a phase transition separating these two phases. In Fig. <ref>, we analyze the Type 2 saddle solution for S^(2) as a function of time and we find that it changes dramatically across V/J = 0.39. When V/J=0.38, it grows linearly in time while when V/J=0.39 or 0.4, it decays with time and saturates to a constant. This result indicates that the transition is first order and occurs close to 0.39. In the volume law phase, the steady state entanglement entropy is determined by the Type 1 saddle solution, which scales linearly with the subsystem size up to L/2. The slope decreases as we increase V/J. As we enter the area law phase, S^(2) is bounded by the Type 2 saddle solution which decays with time and saturates to a finite constant. 
  
  We further analyze the volume law phase close to  the transition point, V/J≲ 0.39, and find that the behavior of S^(2) is slightly different from that of the V≤ 0.3 region at early time. We take V/J=0.36 as an example, see Fig. <ref>. When the subsystem size is small, S^(2)∼α L_A, with the coefficient α=α(t) a decreasing function in time. This decaying  behavior originates from the Type 1 saddle solution and is  responsible for the spike observed in Fig. <ref>. At late time, this coefficient saturates to a finite constant and the steady state still exhibits volume law scaling. 
  
   We summarize the main result for the (4,2) model in the phase diagram shown in Fig. <ref>. Three distinct early time behaviors of S^(2)(L_A) (dashed and dotted lines)  as well as two distinct late time behaviors (red lines) are sketched. Early time entropy always contains a plateau, corresponding to a Type 2 (replica diagonal/quasi-diagonal) solution. The late time, steady state entropy defines either a volume law phase or an area law phase, corresponding to a Type 1 (replica non-diagonal) solution or a Type 2 solution, respectively. In the regime 0.3<V<0.39, we observe a non-monotonic entanglement dynamics at early time for small subsystem size. This behavior has also been observed in the non-unitary dynamics of many systems with small N [We numerically checked several one-dimensional non-unitary dynamics governed by non-Hermitian Hamiltonian and we found the same non-monotonic behavior. ].
  
  
  
  
  
  
  
  
  
  
    §.§.§ Brownian version
  
  
  By assuming equal-time (and otherwise vanishing) correlation for the disorder variables J and/or V (see the definition in Eq.(<ref>)), a Brownian (p,q)=(4,2) model can be written down and studied analogously. Following the general recipe of Subsection <ref>, we consider two Brownian versions: 1) the (4_B,2) model, in which only the inter-cluster SYK_4 interaction J becomes Brownian, and 2) the (4_B,2_B) model, in which both J and the intra-cluster SYK_2 interaction V are Brownian.
  
  
  We find that the (4_B,2) model is in many ways similar to the (4,2) model: the steady state exhibits volume-law entanglement for V≪ J and area-law entanglement  for V≫ J. The two phases are separated by a first order phase transition. We further compare the (4_B,2) model and the (4,2) model as follows:
  
  
    * In the V≪ J regime, both models satisfy the early and late time dynamics (<ref>). The J=0 limit has the same coefficient s_0 = log 2/2. The  (4_B,2) model has a much smaller coefficient v_E.
  
    * The initial conditions Type 2(D) and Type 2(QA/QB) converge to distinct saddle point solutions in the (4,2) model when V<0.3 J. They converge to the same solution in the (4_B,2) model,  which is replica diagonal for all V/J.
  
    * In the V≫ J regime, the steady state of the (4_B,2) model has the scaling form S^(2)∝ J/V. This is different from the S^(2)∝ J^2/V^2 scaling of the (4,2) model and is attributed to the presence of J^-1 in the Brownian correlation relation (<ref>).
  
    * The first order transition occurs at V/J∼ 0.07 for the (4_B,2) model, while it occurs at V/J ≃ 0.39 for the (4,2) model. 
  
   
  Unlike the (4_B,2) model, the (4_B,2_B) model does not exhibit phase transitions. The entanglement dynamics in the entire parameter range of V/J is described by Eq. (<ref>), and the steady state is a volume-law entanglement state [In this case, the replica diagonal saddle point can be solved analytically and its entropy is always linear in t.]. 
  
  
  
  
  
  
  
   §.§ Inter-cluster SYK_2 with intra-cluster SYK_2
  
  
  
  
  Another variant of the model is the (p,q)=(2,2) model. This free fermion model exhibits critical phases distinct from the previous interacting model. We again first consider the limit V=0. The entropy S^(2) as a function of subsystem size and time is shown in Figs. <ref> and <ref>. We observe that S^(2) at early time increases linearly with √(t), and at late time exponentially saturates to a steady state value (see the inset of Fig. <ref>). Numerical evidence suggests that the subsystem scaling of this steady state value agrees with the following expression (see Fig. <ref>)
  
      S^(2)=LS^(2)_SYK2(L_A/L),
  
  where S^(2)_SYK2(λ) is the eigenstate Rényi entanglement entropy for a subsystem with Nλ Majorana modes of a single-cluster SYK_2 model with N Majorana modes <cit.>:
  
      S^(2)_SYK2(λ) = 3/2λlog 2 +
          (1/2- λ) log(√(1 + 4 λ -4 λ ^2)+1-2 λ/1-λ)-1/2log(√(1 + 4 λ -4 λ ^2)+1).
  
  This result is much smaller than the (log 2/2) L_A scaling in the (4,2) model with V=0. 
  
  The entanglement entropy takes a different form when V≠ 0. When V=J=1, numerics (as shown in Fig. <ref>) suggests that the steady state entropy has the form
  
      S^(2)=αlog[L/πsin(π L_A/L)],
  
  which is the same as that for the ground state of a 1D critical system with periodic boundary condition. Furthermore, starting from a product state, S^(2) increases as 
  
      S^(2)=αlog t
  
  at early time with the same coefficient α. Such scaling behavior implies an emergent two-dimensional conformal symmetry with dynamical exponent z=1. Similar behavior has also been found in the random non-unitary free fermion dynamics in the small N limit <cit.>. 
  
  We then tune away from the ratio V/J=1 and find that the critical scaling behavior of S^(2) is retained in a large parameter range. More interestingly, we observe that the coefficient α is linearly proportional to the ratio J/V (see inset of Fig. <ref>), indicating that this model remains critical as long as V/J is finite.
  
  In addition, we compute the second Rényi mutual information I^(2) of the steady state defined as
  
      I^(2) = S^(2)_[x_1,x_2]+S^(2)_[x_3,x_4] - S^(2)_[x_1,x_2]∪ [x_3,x_4].
  
  Here we take periodic boundary condition and partition the system into four connected intervals with end points x_1,2,3,4. We present I^(2) in Fig. <ref> and find that it is only a function of cross-ratio, which is defined as 
      η = sinπ/L|x_1-x_2| sinπ/L|x_3-x_4|/sinπ/L|x_1-x_3|sinπ/L|x_2-x_4|.
  
  This numerical result provides further evidence that this model has conformal symmetry.
  In particular, we find that when η→ 0, I^(2)∼η, the same as that for non-unitary random free fermion dynamics in the small N limit <cit.>. 
  
  
  The very short time behavior of S^(2) in the limit J≪ V can be treated perturbatively, see Appendix <ref>. 
  
  
  
  § DISCUSSION AND CONCLUSION
  
  
  In this work we study from the entanglement perspective the non-unitary dynamics of a 1D SYK chain H=J∑_x H_x,x+1^p-i V∑_x H_x^q. We derive the large N self-consistent saddle point solution to the nth Rényi entropy using the path integral method, and numerically study the second Rényi entropy as a function of time, subsystem size, and 
  the dimensionless coupling strength. We find that for the (p,q)=(4,2) model, as we vary the ratio V/J, the steady state can be either in a volume-law phase or an area-law phase. The steady states of these two phases correspond to distinct saddle point solutions, and the two phases are separated by a first order transition. For the non-interacting (p,q)=(2,2) model, the steady state exhibits critical behavior for any finite ratio V/J, indicating the emergent two-dimensional conformal symmetry in this model. 
  
  There are a lot of possible interesting extensions of our work, which we briefly mention. Firstly, in the (4,2) model, the transition is first order when inter-cluster coupling is described by the regular SYK or Brownian SYK interaction. It would be interesting to study the 1/N correction in this model and check if it can round this first order transition to a second order transition. Secondly, the emergent conformal symmetry observed in (2,2) model may need further investigation. We plan to analytically study the critical saddle point in this large N model and find the connection with the critical phenomena observed in the small N model. Lastly, we can use the method developed in this paper to construct other non-unitary dynamics and explore the exotic phases in them.
  
  
  
  
  Acknowledgment We acknowledge helpful discussions with Leon Balents, Yaodong Li, Andy Lucas and Marcin Szyniszewski. PZ acknowledges support from the Walter Burke Institute for Theoretical Physics at Caltech. CL is supported by the NSF CMMT program under Grants No. DMR-1818533. Use was made of computational facilities purchased with funds from the National Science Foundation
  (CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the
  California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF
  DMR-1720256) at UC Santa Barbara.
  
  
  
  
  
  
  
  
  
  § DERIVATION OF EQ. (<REF>)
  
  
  
  In this section we provide the detailed derivation for the bilocal field action (<ref>). By definition, we have
  
      Z_n[{J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1...i_q}](t,L_A)    = Tr_A( Tr_B [ e^-itH|{0}⟩⟨{0}| e^itH^†] ⋯Tr_B [e^-itH|{0}⟩⟨{0}| e^itH^†] _n copies of Tr_B)
         =
      ∑_ϕ_A⟨ϕ_A| ∑_ϕ^1_B⟨ϕ^1_B| e^ - i t J H^p-tV H^q|{0}⟩⟨{0}| e^i t J H^p-tV H^q|ϕ^1_B⟩⋯
                     ∑_ϕ^n_B⟨ϕ^n_B|  e^- i t J H^p-tV H^q|{0}⟩⟨{0}| e^i t J H^p-tV H^q |ϕ^n_B⟩ |ϕ_A⟩,
  
  where 1_A/B=∑_ϕ_A/B|ϕ_A/B⟩⟨ϕ_A/B| is the completeness relation for the Hilbert subspaces ℋ_A or ℋ_B. Note that when using coherent states, additional minus signs lead to the anti-periodic boundary condition of fermions as for thermal ensembles. We now conduct disorder average over the disorder fields J^x,x+1_i_1...i_p/2j_1...j_p/2 and V^x_i_1i_2...i_q and get
  
      Z_n = ∫∏ dJ^x,x+1_i_1...i_p/2j_1...j_p/2 d V^x_i_1...i_q P(J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1...i_q) Z_n[{J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1...i_q}]=∫ D[χ^x_i] e^-S_n[{χ^x_i}],
  
  here P(J^x,x+1_i_1...i_p/2j_1...j_p/2,V^x_i_1...i_q) is the Gaussian distribution of random parameters and 
  
      S_n    =  ∑_x {∑_i ∫_0^2nt1/2χ^x_i∂_τχ^x_i
      -  .
             .1/2∫_[0,2nt]^2dτ dτ'[
      f(τ)f(τ')J^2/p( ∑_i χ^x_i(τ)χ^x_i(τ'))^p/2( ∑_i χ^x+1_i(τ)χ^x+1_i(τ'))^p/2+V^2/q( ∑_i χ^x_i(τ)χ^x_i(τ'))^q]},
  
  the factor f(τ) = ± i has been introduced in the main text (see below Eq. (<ref>)). The boundary conditions for χ^x_i occur at times τ=0,2t,4t,...,2nt due to the trace over A or B and at times τ=t,3t,...,(2n-1)t due to the existence of |{0}⟩⟨{0}|. The latter conditions can be incorporated by redefinition of the Majorana fields:
  
      χ^x_i(s) = {[            χ^x_2j(2mt+(s-4mt)),                 s∈[4mt,(4m+1)t]; -i χ^x_2j-1((2m+1)t+(4m+1)t-s),           s ∈ [(4m+1)t,(4m+3)t];     -χ^x_2j((2m-1)t+s-(4m+3)t),           s ∈ [(4m+3)t,(4m+4)t] ].
  
  for m=1,2,...,n, with n+1≡ 1. The path defined for χ(s) on s ∈[ 0,4nt) is shown in Fig. <ref> (where the “o” is removed). We see that the only boundary conditions needed to specify for χ are
  
      x∈ A     χ^x_i(4mt^+) = -  χ^x_i((4m+4)t^-),
      
      x∈ B     χ^x_i((4m+2)t^-)= χ^x_i((4m+6)t^+),
      
         χ^x_i(4mt^+)=- χ^x_i((4m+8)t^-).
  
  Eqs. (<ref>) and (<ref>) together define the contour 𝒞. Note that the partial operator will also inherit this boundary condition, which we denote as ∂_s,x. On contour 𝒞 the direction of real-time evolution has the explicit form
  
      f(s) = {[ + i   for  s ∈[(2m+1)t,(2m+2)t),;     -i   for  s ∈ [2mt,(2m+1)t). ].
  
  
  Next we introduce bilocal fields
  
      1 = ∫_𝒞 [D G_x] [DΣ_x] e^-  1/2∑_x ∫ d τ d τ' Σ_x(τ,τ')(N/2G_x(τ,τ') - ∑_i^N/2χ^x_i(τ) χ^x_i(τ') ),
  
  and plug it into the action (<ref>). This allows us to integrate out the Majorana fields. Note that the projector P(s,s') (see Eq. (<ref>)) needs to be introduced during this procedure. This finally leads to Eq. (<ref>) in the main text. Note that in the main text and following appendices we remove all the “o” to keep the notation concise.
  
  
  
  § REPLICA SYMMETRY, TWIST FIELDS, AND PERTURBATIVE LIMIT 
  
  
  
  One observation from Fig. (<ref>) and the boundary conditions (<ref>) for χ^x_i is that if we redefine the fields by
  
      χ^x_i(s) = {[                     χ^x_i(s+4t), x ∈ A and s ∈ [(4m+2)t,(4m+4)t),;                        χ^x_i(s),                        otherwise ].
  
  then all the χ^x_i carries the same boundary conditions that were original carried only by χ^x_i(s) with x∈ B. Formally χ^x_i(s) is related to χ^x_i(s') by some orthogonal matrix 𝒜(s,s') when x ∈ A[This orthogonal matrix is a version of the twist field; or in more rigorous sense, it implements the temporal branch cut that connect two twist fields at the boundary x^*.
  Viewing the n Majorana fields from the n replicas as separate fields, we can introduce the twist field 𝒯 to incorporate the boundary condition between the n fields, defined by
  
      𝒯_x^*(t)χ^x_i,m(τ)→χ^x_i,m+1(τ),   x >x^*, τ =t,
  
  while 𝒯_x^* acts trivially on χ^x_i,m(τ) otherwise. Then we have
  
      𝒯^†_x^*(t)χ^x_i,m(τ)→χ^x_i,m-1(τ),   x> x^*,τ = t.
  
  This way we have
  
      Z_n(t)   =∫_𝒞[D χ^x_i,m]𝒯_x^*+L_A(4t)𝒯^†_x^*+L_A(2t)𝒯^†_x^*(4t)𝒯_x^*(2t) e^-NS
         = Tr[𝒯_x^*+L_A𝒯^†_x^*ρ^⊗ n𝒯^†_x^*+L_A𝒯_x^*ρ^⊗ n].
  
  ]. This redefinition reveals some sort of replica “symmetry” in the partition function which allows to calculate the entropy in the small J limit perturbatively.
  In the small J limit, the saddle point solution for G^n_x is replica diagonal and is to leading order of J well approximated by the J=0 solution: G^n_x = 1_n× n⊗ G_q, where G_q is the Green's function for a single-site SYK_q cluster.  These guarantee that Z_n(t,L_A) looks almost exactly like n decoupled copies of Z(t), except for local terms at subsystem boundaries
  
      δ S    =  S_n,saddle -  n S_1,saddle
         =J^2/4p∑_x^*=0,L_A{-∫_[0,4nt]^2dsds' f(s)f(s')
      ((𝒜^TG^n_x^*𝒜)(s,s') G^n_x^*+1(s,s'))^p/2P + n∫_[0,4t]^2 dsds' f(s)f(s')
      (G^1_x^* G^1_x^*+1)^p/2P},
  
  where x^* are summing over ∂ A = {0,L_A}, and we denoted the saddle point solution of Z(t) by G^1_x.
  
   Using the expressions of 𝒜 and P, δ S can be simplified to
  
      S^(2) = δ S/n-1
      = 2× 4n/n-1J^2/4p∫^t_0 ds ∫^2t_t ds' f(s)f(s')G_q^p(s-s'),
  
  where the prefactor 4n is the number of disjoint sectors (blocks) in the Green's functions on which the nonzero entries of the projected Green's functions G_AP = 𝒜^TG_x^n 𝒜P and G_BP=G_x^nP do not overlap. 
  
  
  
  
  
  
  § PERTURBATIVE RESULT FOR THE (4,4) AND (2,2) MODELS
  
  
  
  
  In this Appendix we outline the main result for the (p,q)=(4,4) model and additional perturbative result for the (2,2) model. For the (4,4) model, the strong inter-cluster coupling limit J≫ V shows similar behavior as the (4,2) model: both are described by Eq. (<ref>). As one moves away from this limit, S^(2) decays with increasing V just as in the (4,2) model, however the decay is much slower and no drastic change of slope of S^(2)(V/J) is observed in our numerics. 
  
  The weak coupling limit J≪ V is more interesting, in which S^(2) shows an early time behavior that scales linearly with log t, see Fig. <ref>.
  This scaling can be understood from Eq. (<ref>). Using the conformal Green's function for q=4
  
      G_c(s,s') =  1/(2V √(π)2t/πsinπ(s-s')/2t)^1/2
  
  and plug it into Eq. (<ref>), we have
  
      S^(2) =  J^2/2 π V^2log2t/πϵ
      - π  J^2 ϵ^2/24 V^2 t^2  +  O((ϵ/t)^3),
  
  The logarithmic t behavior observed in Fig. <ref> is then identified with the first term.
  
  We now apply the same method to the (2,2) model in the weak coupling limit J≪ V, in which an early time S^(2)∝log t behavior is also observed (see Fig. <ref>). The scaling behavior can again be obtained by plugging the conformal Green's function (<ref>) to Eq. (<ref>). Note that since for generic q we have G_c∝(sinπ(s-s')/2t)^-2/q while the integrand in Eq. (<ref>) is p-power of G_c, we expect that all the p=q model exhibits the same behavior (up to a numerical constant) for any value of p=q. For q=p=2 we get
  
      S^(2)
      =
      4 J^2/π^2 V^2log2t/πϵ
      - J^2 ϵ^2/3 V^2 t^2
      +O((ϵ/t)^3),
  
  which indeed is the same as the (4,4) result Eq. (<ref>) up to a numerical prefactor. 
  
  To conclude, the early time behavior of S^(2) at J/V≪ 1 for both models is in good agreement with the perturbative results (<ref>) and (<ref>). We emphasize that these perturbative results should only work for V^-1≪ t≪Λ, where a large time cutoff Λ needs to be introduced, since we know that S^(2) must be bounded at long time. We also note that although early time log t behavior appears for both J/V≪ 1 and J/V ∼ 1 regimes in the (2,2) model, they are still qualitatively different since the prefactor α∝ J^2/V^2 in one regime while α∝ J/V in the other.
  
  
  
  
  
  
  
  
  
  