11institutetext: Condensed Matter Theory Center and Joint Quantum Institute,
Department of Physics, University of Maryland, College Park, Maryland 20742,
USA22institutetext: Department of Physics, Brandeis University, Waltham,
Massachusetts 02453, USA

# Note on entropy dynamics in the Brownian SYK model

Shao-Kai Jian 2,1 Brian Swingle<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the time evolution of Rényi entropy in a system of two coupled
Brownian SYK clusters evolving from an initial product state. The Rényi
entropy of one cluster grows linearly and then saturates to the coarse grained
entropy. This Page curve is obtained by two different methods, a path integral
saddle point analysis and an operator dynamics analysis. Using the Brownian
character of the dynamics, we derive a master equation which controls the
operator dynamics and gives the Page curve for purity. Insight into the
physics of this complicated master equation is provided by a complementary
path integral method: replica diagonal and non-diagonal saddles are
responsible for the linear growth and saturation of Rényi entropy,
respectively.

##  1 Introduction

Interest in the entropy dynamics of quantum many-body systems has increased
dramatically in the past few years. In the context of the black hole
information problem Hawking:1974particle , recent works have derived the Page
curve Page:1993average  of an evaporating black hole from holographic
calculations Penington:2019entanglement ; Almheiri:2019the ;
Penington:2019replica ; Almheiri:2019replica . These calculations have been
quickly generalized to various related situations Almheiri:2020page ;
Almheiri:2019islands ; Almheiri:2019entanglement ; Rozali:2019information ;
Gautason:2020page ; Hashimoto:2020islands ; Hartman:2020flat ;
Hollowood:2020islands ; Krishnan:2020Page ; Chen:2020evaporating ;
Chen:2020bra ; Hartman:2020islands ; Anegawa:2020notes ; Akers:2020leading ;
Balasubramanian:2020islands ; Balasubramanian:2020entanglement ;
Ling:2020island ; Bhattacharya:2020topological ; Marolf:2020observation . In
particular, using a semiclassical saddle point analysis of the gravitational
path integral, it was shown that replica wormhole configurations, though
exponentially suppressed, become dominant after the Page time
Penington:2019replica ; Almheiri:2019replica . Meanwhile, entropy dynamics are
also much explored in various circuit models Dahlsten:2007emergence ;
Znidaric:2007exact ; Lashkari:2011towards ; Nahum:2016quantum ;
vonKeyserlingk:2017operator ; Nahum:2018operator ; Rakovszky:2018diffusive ;
Vedika2018:operator ; Chan:2018solution ; Chan:2018spectral ;
Bertini:2018exact ; Zhou:2019operator ; Zhuang:2019scrambling . Notably, Page-
like behavior has been obtained in random circuit models Czech:2011black ;
Mathur:2011correlations ; Bradler:2015one ; Tokusumi:2018quantum ;
Piroli:2020a ; Liu:2020a , including the Brownian SYK model
Sunderhauf:2019quantum . It is thus interesting to make a connection between
these two types of methods in a single model.

Since the Brownian SYK model is amenable to both saddle point methods and
circuit techniques, in this note we report a calculation of the Page curve in
a model consisting of two coupled Brownian SYK clusters using both saddle
point methods and operator dynamics. The quantity we are interested in is the
Rény entropy of one cluster after tracing out the other. To formulate a path
integral representation of the Rényi entropy, the initial state is taken to be
a tensor product of thermofield double (TFD) states in each subsystem obtained
by doubling the Hilbert space to left and right sides111The Brownian model
does not have a fixed Hamiltonian with which we can define a finite
temperature, so we consider the infinite-temperature TFD state in each
subsystem, i.e., a maximally entangled state.. This is the Brownian version of
the setup in Penington:2019replica : the entanglement between left and right
sides is maximal, whilst it is initially zero between the two subsystems. We
find analytic solutions to the saddle point equations and show that replica
diagonal and non-diagonal solutions are responsible for early time growth and
the late time saturation of the Rényi entropy, respectively.

On the other hand, the density matrix dynamics can be directly analyzed using
an operator dynamics approach. For simplicity we take a tensor product of
Kourkoulou-Maldacena states Kourkoulou:2017pure  in each subsystem as initial
state. Again, there is no entanglement between the two subsystems initially.
Then the system is evolved under the full Hamiltonian, and the Page curve is
obtained from a corresponding master equation. We compare the results from
these two methods, and find excellent agreement even though the two approaches
use slight different initial states. Complementing the saddle point method,
the master equation knows about the microstate from the perspective of
symmetry, i.e., the Fermi parity in our case. The Hilbert space factorizes
into different Fermi parity sectors, leading to an order one correction to the
coarse grained entropy. Note that while the operator dynamics approach gives
access only to the second Rényi entropy, in the saddle point method, we are
able to get solutions for Rényi entropy for arbitrary Rényi index $n$.

The rest of the paper is organized as follows. In Section 2 we discuss a
saddle point analysis of the path integral representation of Rényi entropy.
Both replica diagonal and replica non-diagonal solutions are obtained
analytically and checked numerically. The Page curve is obtained using these
solutions, with the replica diagonal solution being responsible for the linear
growth and the replica non-diagonal solution leading to the saturation to the
coarse grained entropy. In section 3 we study the Page curve using the
operator dynamics of the Brownian SYK model. We derive the master equation
governing the operator size distribution function. The initial linear growth
and late time saturation can be obtained analytically from the master
equation, which shows exact agreement with the saddle point analysis.

##  2 Rényi entropy dynamics from saddle points

###  2.1 Coupled Brownian SYK clusters

The time-dependent Hamiltonian of two coupled Brownian SYK clusters labelled
by $a=1,2$ is given by

$\displaystyle
H(t)=\sum_{|A|=q,a=1,2}J_{A}^{a}(t)[\psi_{a}]_{A}+\sum_{|A|=|B|=q/2}V_{A,B}(t)[\psi_{1}]_{A}[\psi_{2}]_{B},$
(1)

where $A=j_{1}...j_{|A|}$ denotes an ascending list of length $|A|$, $q$ is an
even integer, and

$\displaystyle[\psi_{a}]_{A}\equiv
i^{|A|/2}\psi_{j_{1},a}\psi_{j_{2},a}...\psi_{j_{|A|},a},$ (2)

is a short-hand notation for an $|A|$-body interaction. The $\psi_{j,a}$,
$j=1,...,N_{a}$, are Majorana fermions in subsystem $a$, and satisfy
$\\{\psi_{j,a},\psi_{j^{\prime},a^{\prime}}\\}=\delta_{jj^{\prime}}\delta_{aa^{\prime}}$.
The summations in Eq. (1) are over all possible lists with the indicated
number of fermions. $J_{A}^{a}(t)$ and $V_{A,B}(t)$ are Brownian random
interactions within and between the two subsystems, respectively. The
interaction strength is drawn from a Gaussian distribution with mean zero and
variance given by

$\displaystyle\overline{J_{A}^{a}(t)J_{A^{\prime}}^{a^{\prime}}(t^{\prime})}$
$\displaystyle=$
$\displaystyle\frac{2^{q-1}q!}{q^{2}N_{a}^{q-1}}\mathcal{J}\delta(t-t^{\prime})\delta_{A,A^{\prime}}\delta_{a,a^{\prime}},$
(3) $\displaystyle\overline{V_{A,B}(t)V_{A^{\prime},B^{\prime}}(t^{\prime})}$
$\displaystyle=$
$\displaystyle\frac{2^{q}(q/2)!^{2}}{q^{2}N_{1}^{(q-1)/2}N_{2}^{(q-1)/2}}\mathcal{V}\delta(t-t^{\prime})\delta_{A,A^{\prime}}\delta_{B,B^{\prime}},$
(4) $\displaystyle\delta_{A,A^{\prime}}$ $\displaystyle\equiv$
$\displaystyle\delta_{j_{1},j_{1}^{\prime}}...\delta_{j_{|A|},j_{|A^{\prime}|}^{\prime}}.$
(5)

The over line denotes an average over the Gaussian distribution of couplings.
The interaction strength has dimension one (the dimension of energy), so
$\mathcal{J}$ and $\mathcal{V}$ also have dimension one, while the
$\delta$-function makes up another dimension one and also indicates that the
couplings are Brownian variables, uncorrelated in time. Regarding the
prefactor, the dependence on $N_{a}$ is chosen to facilitate the large-$N$
limit and the dependence on $q$ is chosen to facilitate the large-$q$
expansion in Appendix D. In general, the coupling between two subsystem does
not have to be the same $q$-body interaction as the interaction within each
subsystems, but we make such a choice for simplicity.

###  2.2 Setup

To investigate the entropy dynamics, we consider a similar setup to Ref.
Penington:2019replica : starting from the tensor product of two thermofield
double (TFD) states in each of the subsystems ($a=1,2$), we focus on the Rényi
entropy of subsystem $a=1$ by tracing out subsystem $a=2$. Because Brownian
random interactions do not conserve energy, we simply consider an infinite
temperature TFD state, which is a maximally entangled state. To prepare such a
state, we double the Hilbert space by introducting left ($L$) and right ($R$)
copies of the fermions, $\psi_{j,a,L}$ and $\psi_{j,a,R}$, for both subsystems
$a=1,2$. Then the maximally entangled state and the initial density matrix are
given by

$\displaystyle(\psi_{j,a,L}+i\psi_{j,a,R})|\infty\rangle=0,\quad\forall
a=1,2,\quad\forall j=1,...,N,\quad\rho_{0}=|\infty\rangle\langle\infty|.$ (6)

Consider a time evolution generated by the sum of left and right Hamiltonians.
The random couplings are identical between the two sides, up to an overall
coefficient, with $H_{L}(t)=H(t;\psi_{j,a,L})$ and
$H_{R}(t)=(-1)^{q/2}H(t;\psi_{j,a,R})$. This choice implies that
$H_{R}|\infty\rangle=H_{L}|\infty\rangle$. Hence, the reduced density matrix
$\rho_{1}$ of the subsystem $a=1$ (including both $L$ and $R$ pieces) at time
$t$ is

$\displaystyle\rho(t)=U(t)\rho_{0}U^{\dagger}(t),\quad
U(t)=\mathcal{T}e^{-i\int_{0}^{t}dt^{\prime}(H_{L}(t^{\prime})+H_{R}(t^{\prime}))},\quad\rho_{1}(t)=\text{Tr}_{2}[\rho(t)].$
(7)

where $\mathcal{T}$ denotes time ordering, and $\text{Tr}_{a}$ denotes the
trace over subsystem $a$. The $n$-th Rényi entropy is
$e^{-(n-1)S_{n}}=\text{Tr}_{1}[\rho_{1}(t)^{n}]$. This joint left-right
evolution is equivalent to a single sided evolution for twice the time, i.e.,
$U(t)=\mathcal{T}e^{-i\int_{0}^{2t}dt^{\prime}H(t^{\prime})}$.

A path integral representation of the trace of the $n$-th power of the reduced
density matrix is obtained via a standard replica trick using $n$ copies of
the system and twist fields to implement modified boundary conditions on
subsystem $1$. This formulation gives
$\text{Tr}_{1}[\rho_{1}^{n}]=\frac{Z_{(n)}}{Z_{(1)}^{n}}$, where $Z_{(n)}$,
$n\geq 2$ is the replicated partition function with twist operators inserted,
and $Z_{(1)}\equiv Z$ is the partition function of a single replica. The
replicated partition function $Z_{(n)}=\int[D\psi]e^{-I}$ can be implemented
in a Keldysh contour with two twist operators at $t=0$ and $t=T$, respectively
Penington:2019replica . The insertion of twist operators in the contour is
shown in Fig. 1. The effective action for the path integral of the replicated
systems is

$\displaystyle
I=\sum_{s=\pm}s\int_{0}^{T}dt\Big{(}\frac{1}{2}\sum_{j,a,\alpha}\psi_{j,a,s}^{\alpha}\partial_{t}\psi_{j,a,s}^{\alpha}+i\big{(}\sum_{A,a,\alpha}J_{A}^{a}[\psi_{a,s}^{\alpha}]_{A}+\sum_{A,B,\alpha,\beta}g_{s}^{\alpha\beta}V_{A,B}[\psi_{1,s}^{\alpha}]_{A}[\psi_{2,s}^{\beta}]_{B}\big{)}\Big{)},$
(8)

where $s=\pm$ stands for the forward and backward contour,
$\alpha,\beta=1,2,...,n$ are the replica indices, and
$g^{\alpha\beta}_{+}=\delta^{\alpha\beta}$,
$g^{\alpha\beta}_{-}=\delta^{\alpha+1,\beta}\equiv\epsilon^{\alpha\beta}$ is
due to the twist operator.

Figure 1: The schematic plot of the Keldysh contour and twist operator. The
blue and red contours represent subsystem $a=1$ and $a=2$, respectively. The
black dots at two ends of contour represent twist operators. The red contour
is the $\alpha$ replica while in $s=-$ of the blue contour it changes to
$\alpha-1$ due to the twist operator.

At this point, the Brownian random interactions are integrated out. Strictly
speaking, what we calculate is the logarithm of the average of the replicated
partition function, i.e.,
$\frac{1}{1-n}\log\overline{\text{Tr}_{1}[\rho_{1}^{n}]}$, instead of the
average Rényi entropy
$\frac{1}{1-n}\overline{\log\text{Tr}_{1}[\rho_{1}^{n}]}$ which involves an
additional replica trick to obtain the averaged logarithm. Nevertheless, due
to the large-$N$ structure of the Brownian models, we expect the circuit-to-
circuit fluctuation are suppressed vonKeyserlingk:2017operator ;
Zhou:2018emergent  so that both quantities agree with each other at large $N$.
After integrating out the Brownian variables and introducing the bilocal
fields $G$ and $\Sigma$,

$\displaystyle\int
d\hat{G}_{a}\delta\Big{(}G_{a,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2})-\frac{1}{N_{a}}\sum_{j}\psi_{j,a,s}^{\alpha}(t_{1})\psi_{j,a,s^{\prime}}^{\beta}(t_{2})\Big{)}$
(9) $\displaystyle=\int
d\hat{G}_{a}d\hat{\Sigma}_{a}\exp\Big{[}-\frac{N_{a}}{2}\Sigma_{a,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2})\Big{(}G_{a,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2})-\frac{1}{N_{a}}\sum_{j}\psi_{j,a,s}^{\alpha}(t_{1})\psi_{j,a,s^{\prime}}^{\beta}(t_{2})\Big{)}\Big{]},$
(10)

we arrive at the following effective action,

$\displaystyle-I$ $\displaystyle=$
$\displaystyle\sum_{a}N_{a}\Big{[}\log\text{Pf}[\partial_{t}\hat{\sigma}^{z}-\hat{\Sigma}_{a}]+\int
dt_{1}dt_{2}\big{(}-\frac{1}{2}\Sigma_{a,ss^{\prime}}^{\alpha\beta}G_{a,ss^{\prime}}^{\alpha\beta}+\frac{\mathcal{J}}{4q^{2}}\delta(t_{12})c_{ss^{\prime}}(2G_{a,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2}))^{q}\big{)}\Big{]}$
(11) $\displaystyle+\sqrt{N_{1}N_{2}}\frac{\mathcal{V}}{2q^{2}}\int
dt_{1}dt_{2}\delta(t_{12})c_{ss^{\prime}}(2G_{1,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2}))^{q/2}g^{\alpha\gamma}_{s}g^{\beta\delta}_{s^{\prime}}(2G_{2,ss^{\prime}}^{\gamma\delta}(t_{1},t_{2}))^{q/2}.$

where $t_{12}\equiv t_{1}-t_{2}$, and $c_{++}=c_{--}=-1$, $c_{+-}=c_{-+}=1$ is
due to the Keldysh evolution. The summation over the replica indices and the
contour indices is implicit. $\hat{\sigma}^{z}$ denotes the Pauli matrix
acting on contour space.

From the above effective action the Schwinger-Dyson equations are

$\displaystyle\hat{G}^{-1}_{a}$ $\displaystyle=$
$\displaystyle\hat{\sigma}^{z}\partial_{t}-\hat{\Sigma}_{a},$ (12)
$\displaystyle\Sigma_{a,ss^{\prime}}^{\alpha\beta}$ $\displaystyle=$
$\displaystyle
c_{ss^{\prime}}\delta(t_{12})\Big{[}\frac{\mathcal{J}}{q}(2G_{a,ss^{\prime}}^{\alpha\beta})^{q-1}+\sum_{\gamma\delta}\sqrt{\frac{N_{\bar{a}}}{N_{a}}}\frac{\mathcal{V}}{q}(2G_{a,ss^{\prime}}^{\alpha\beta})^{q/2-1}[g_{(a)}]^{\alpha\gamma}_{s}[g_{(a)}]^{\beta\delta}_{s^{\prime}}(2G_{\bar{a},ss^{\prime}}^{\gamma\delta})^{q/2}\Big{]},$

where we have defined $g_{(1)}=g$, $g_{(2)}=g^{T}$, and $\bar{1}=2$,
$\bar{2}=1$.

###  2.3 Saddle point solutions

For simplicity, we assume the two clusters have equal numbers of Majorana
fermions $N_{1}=N_{2}=N$. To look for a replica diagonal solution, we can
start by looking for a solution when the inter-cluster coupling
$\mathcal{V}=0$. In this case, the problem reduces to $n$ independent
replicas. Moreover, because of the Brownian nature of the problem, the self-
energy is local in time (2.2). Starting from the Green’s functions ansatz

$G_{a}^{\alpha\beta}(t_{1},t_{2})=\delta^{\alpha\beta}\frac{f(t_{12})}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right),$ (14)

one uses (2.2) to get

$\displaystyle\Sigma_{a}^{\alpha\beta}(t_{1},t_{2})=\delta^{\alpha\beta}f(0)^{q-1}\frac{\mathcal{J}}{q}\delta(t_{12})\left(\begin{array}[]{cccc}0&-1\\\
1&0\end{array}\right),\quad\Sigma_{a}^{\alpha\beta}(\omega)=\delta^{\alpha\beta}f(0)^{q-1}\frac{\mathcal{J}}{q}\left(\begin{array}[]{cccc}0&-1\\\
1&0\end{array}\right),$ (19)

where $\hat{\Sigma}_{a}(\omega)=\int dt_{12}\hat{\Sigma}_{a}(t_{12})e^{i\omega
t_{12}}$. Here the limit $\mathcal{J}T\gg 1$ is implicit, as we are interested
in the long-time behaviors of Rényi entropy, so the Fourier transform becomes
an integral. We also solve the Schwinger-Dyson equation (12, 2.2) numerically
in Appendix B for finite $T$, and find excellent agreement with the analytic
solution we give in the following.

Plugging the self-energy into (12) gives

$\displaystyle G_{a}^{\alpha\beta}(\omega)$ $\displaystyle=$
$\displaystyle\delta^{\alpha\beta}\frac{1}{\omega^{2}+f(0)^{2q-2}(\frac{\mathcal{J}}{q})^{2}}\left(\begin{array}[]{cccc}i\omega&-f(0)^{q-1}\frac{\mathcal{J}}{q}\\\
f(0)^{q-1}\frac{\mathcal{J}}{q}&-i\omega\end{array}\right),$ (22)
$\displaystyle G_{a}^{\alpha\beta}(t_{1},t_{2})$ $\displaystyle=$
$\displaystyle\delta^{\alpha\beta}\frac{e^{-f(0)^{q-1}\frac{\mathcal{J}}{q}|t_{12}|}}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right).$ (25)

Comparing it to the ansatz, we find
$f(t_{12})=e^{-\frac{\mathcal{J}}{q}|t_{12}|}$, so the replica diagonal
solution is

$\displaystyle
G_{a}^{\alpha\beta}(t_{1},t_{2})=\frac{e^{-\frac{\mathcal{J}}{q}|t_{12}|}}{2}\delta^{\alpha\beta}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right),\quad\Sigma_{a}^{\alpha\beta}(t_{1},t_{2})=\frac{\mathcal{J}}{q}\delta^{\alpha\beta}\delta(t_{12})\left(\begin{array}[]{cccc}0&-1\\\
1&0\end{array}\right).$ (30)

Now we claim that the above function (30) is still a solution to the
Schwinger-Dyson equation with finite coupling $\mathcal{V}>0$. The reason is
two-fold: (1) for the replica-diagonal solution, the $\mathcal{V}$ term in
(2.2) is only non-vanishing on intra Keldysh contour $s=s^{\prime}$ term due
to the twist operator, and (2) the self-energy only depends on the Green’s
function at $t_{12}=0$ which is vanishing on the intra Keldysh contour
$s=s^{\prime}$ term. So if one plugs (30) into the Schwinger-Dyson equation,
the $\mathcal{V}$ term in (2.2) vanishes and the it solves the equation.

Besides the replica diagonal solution, the twist operator induces new replica
non-diagonal solutions that are the analog of the wormhole solutions found in
Penington:2019replica . We assume the subsystem $a=2$ still hosts the replica
diagonal solution,
$G_{2}^{\alpha\beta},\Sigma_{2}^{\alpha\beta}\propto\delta^{\alpha\beta}$, and
find that the subsystem $a=1$ supports a replica non-diagonal solution. To get
a replica non-diagonal solution, the nontrivial part must come from the twist
operator in (2.2), corresponding to the inter Keldysh contour correlation
function that crosses the twist operator as shown in Fig. 1. The self-energy
of subsystem $a=1$ is given by

$\displaystyle\Sigma_{1,-+}^{\alpha\beta}=\frac{\delta(t_{12})}{q}\Big{[}\mathcal{J}(2G_{1,-+}^{\alpha\beta})^{q-1}+\mathcal{V}(2G_{1,-+}^{\alpha\beta})^{q/2-1}\sum_{\gamma}\epsilon^{\alpha\gamma}(2G_{2,-+}^{\gamma\beta})^{q/2}\Big{]}.$
(31)

For a diagonal solution
$\Sigma_{1,-+}^{\alpha\beta}\propto\delta^{\alpha\beta}$, the second term in
(31) vanishes, and the equation reduces to the diagonal solution (30) as we
have shown. However, the second term in (31) also suggests a replica non-
diagonal solution, $\Sigma_{1,-+}^{\alpha\beta}\propto\epsilon^{\alpha\beta}$.

This leads us to consider a replica non-diagonal ansatz for subsystem $a=1$
(and a replica diagonal ansatz for subsystem $a=2$),

$\displaystyle\hat{G}_{1}=\frac{f_{1}(t_{12})}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-\tilde{\epsilon}^{T}\\\
\tilde{\epsilon}&-\text{sgn}(t_{12})\end{array}\right),\quad\hat{G}_{2}=\frac{f_{2}(t_{12})}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right),$ (36)

where the replica indices are implicit, and
$\tilde{\epsilon}^{\alpha\beta}\equiv\text{sgn}(\alpha-\beta)\delta^{\alpha+1,\beta}$.
Note that $\tilde{\epsilon}^{T}\epsilon=1$. The sign prefactor in
$\tilde{\epsilon}$ is due to an emergent time ordering between different
replicas when different replicas develop nonvanishing correlations. This is
also confirmed by the numerical solutions (see Fig. 2). A detailed derivation
is given in Appendix A, where one finds that $f_{1}$ and $f_{2}$ are
$f_{1}(t_{12})=f_{2}(t_{12})=e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|t_{12}|}$.
The replica non-diagonal solutions read

$\displaystyle\hat{G}_{1}=\frac{e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|t_{12}|}}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-\tilde{\epsilon}^{T}\\\
\tilde{\epsilon}&-\text{sgn}(t_{12})\end{array}\right),\quad\hat{G}_{2}=\frac{e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|t_{12}|}}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right).$ (41)

Thus, we have found both replica diagonal and non-diagonal solutions. We check
these solutions by numerically iterating the Schwinger-Dyson equation (12,
2.2). In the numerical calculations, we focus on $n=2,3$. As shown in Fig. 2,
the agreement between the analytics and the numerical solutions is quite good.
Note the small contributions near the twist operators in the replica non-
diagonal case. The analytical result above corresponds to the limit of large
$T$ where these boundary contributions can be neglected, but they do give
important contributions to the on-shell action as we discuss below.

Figure 2: The numerical replica diagonal and non-diagonal solutions for $n=2$
(a,b), and for $n=3$ (c,d). The coordinates are arranged such that
$(2kT,(2k+2)T)$ belongs to the $k$-th replica, and within each replica,
$(2kT,(2k+1)T)$ [$((2k+1)T,(2k+2)T)$] is the forward $s=+$ (backward $s=-$)
contour. Small nonvanishing values appear at the locations of the twist
operators. This has important contributions to the onshell action, and we
discuss it in Appendix C.

###  2.4 Page curve from saddle points

We now show that the replica diagonal and non-diagonal solutions lead to the
linear increase and the saturation of the Rényi entropy, respectively. We
first give analytic results for Rényi entropy from the two saddle point
solutions, and then numerically evaluate the onshell action to verify our
analytic results.

Recall that the replica diagonal solution (30) is actually the same as the
solution of Schwinger-Dyson equation (2.2) when $\mathcal{V}=0$. This means
the replica diagonal solution is a solution of replicated action without
inserting the twist operator because the twist operator effect is proportional
to $\mathcal{V}$. As a result, the first line in (11) counts the total Hilbert
space dimension. For each replica system, we have two subsystems $a=1,2$ each
hosting $N$ Majorana fermions (remember we have doubled the Hilbert space to
prepare the maximally entangled initial state within each subsystems), thus
the total Hilbert dimension is $2^{nN}$. In the second line of (11), for
replica diagonal solution the inter Kelysh contour component $s\neq
s^{\prime}$ is zero. The intra Keldysh contour component $s=s^{\prime}$ leads
to a linear increase of Renyi-$n$ entropy222The factor
$\delta(x)\text{sgn}(x)^{q}$ seems to give a vanishing result. But we can
consider a smeared out $\delta$-function, which leads to a nonvanishing
result. We verify numerically in Appendix C that this smearing procedure
indeed gives the correct linear growth., i.e.,

$\displaystyle-\frac{I^{(1)}}{N}$ $\displaystyle=$ $\displaystyle n\log
2-\frac{n\mathcal{V}}{q^{2}}\int
dt_{1}dt_{2}\delta(t_{12})\big{(}e^{-\frac{\mathcal{J}}{q}|t_{12}|}\text{sgn}(t_{12})\big{)}^{q}=n\log
2-\frac{n\mathcal{V}T}{q^{2}},$ (42)
$\displaystyle\quad\frac{e^{-I^{(1)}}}{Z^{n}}$ $\displaystyle=$ $\displaystyle
e^{-\frac{nN\mathcal{V}T}{q^{2}}},$ (43)

where the first term in $I^{(1)}$ is canceled by the denominator which is the
Hilbert space of $n$ replicated systems $Z^{n}=2^{nN}$.

For the replica non-diagonal solution, we first show that it leads to a time-
independent on-shell action333For finite $T$, we need to consider the effects
of the twist operators which serve as boundary conditions at the end of the
Keldysh contour. The on-shell action from the replica non-diagonal solution is
then not time independent and will receive important corrections. We discuss
this correction in Sec. 2.5 and in Appendix D.. The on-shell action is a
function of $\mathcal{J}T$ and $\mathcal{V}T$,
$I=I[\mathcal{J}T,\mathcal{V}T]$, so its time derivative reads

$\displaystyle\frac{dI}{dT}$ $\displaystyle=$
$\displaystyle\frac{\mathcal{J}}{T}\frac{\partial
I}{\partial\mathcal{J}}+\frac{\mathcal{V}}{T}\frac{\partial
I}{\partial\mathcal{V}}$ (44) $\displaystyle=$ $\displaystyle-\frac{N}{T}\int
dt_{1}dt_{2}\delta(t_{12})c_{ss^{\prime}}\Big{[}\frac{\mathcal{J}}{4q^{2}}\sum_{a}(2G_{a,ss^{\prime}}^{\alpha\beta})^{q}+\frac{\mathcal{V}}{2q^{2}}(2G_{1,ss^{\prime}}^{\alpha\beta})^{q/2}g^{\alpha\gamma}_{s}g^{\beta\delta}_{s^{\prime}}(2G_{2,ss^{\prime}}^{\gamma\delta})^{q/2}\Big{]}.$
(45)

Plugging the replica non-diagonal solution (41) into the time derivative, we
have

$\displaystyle\frac{dI^{(2)}}{dT}$ $\displaystyle=$
$\displaystyle\frac{N}{T}\int
dt_{1}dt_{2}\delta(t_{12})(1+1-1-1)\Big{[}\frac{\mathcal{J}+\mathcal{V}}{2q^{2}}\big{(}e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|t_{12}|}\big{)}^{q}\Big{]}=0,$
(46)

where the vanishing prefactor is coming from summing over $c_{ss^{\prime}}$.
The difference between replica diagonal and non-diagonal solutions originates
from the twist operator, and the non-diagonal solution has a nontrivial
contribution from the inter Keldysh contour component, so the time derivative
vanishes. Because the bulk of the on-shell action vanishes, its value is
determined by boundary effects near the twist operators.

Explicitly, the onshell action for the replica non-diagonal solution is given
by

$\displaystyle-\frac{I^{(2)}}{N}=\sum_{a}\log\text{Pf}[\hat{\sigma}^{z}\partial_{t}-\hat{\Sigma}_{a}]$
$\displaystyle+(1-q)\int\delta(t_{12})c_{ss^{\prime}}\Big{[}\sum_{a}\frac{\mathcal{J}}{4q^{2}}(2G_{a,ss^{\prime}}^{\alpha\beta})^{q}+\frac{\mathcal{V}}{2q^{2}}(2G_{1,ss^{\prime}}^{\alpha\beta})^{q/2}g^{\alpha\gamma}_{s}g^{\beta\delta}_{s^{\prime}}(2G_{2,ss^{\prime}}^{\gamma\delta})^{q/2}\Big{]},$
(47) $\displaystyle\qquad\quad=\sum_{a}\log\text{Pf}[\hat{G}_{a}^{-1}],$ (48)

where the second line vanishes for the same reason as in (46), and we use (12)
to get the third line. Hence, the on-shell action is determined by the
Pfaffian of the Green’s function.

It is convenient to approach the calculation of the Pfaffian by recalling that
the $a=2$ subsystem still hosts the replica diagonal solution, which is also a
solution when there is no twist operator in (11). Consider the action without
twist operators,

$\displaystyle-\frac{I^{(0)}}{N}$ $\displaystyle=$
$\displaystyle\sum_{a}\Big{[}\log\text{Pf}[\partial_{t}\hat{\sigma}^{z}-\hat{\Sigma}_{a}]+\int
dt_{1}dt_{2}\big{(}-\frac{1}{2}\Sigma_{a,ss^{\prime}}^{\alpha\beta}G_{a,ss^{\prime}}^{\alpha\beta}+\frac{\mathcal{J}}{4q^{2}}\delta(t_{12})c_{ss^{\prime}}(2G_{a,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2}))^{q}\big{)}\Big{]}$
(49) $\displaystyle+\frac{\mathcal{V}}{2q^{2}}\int
dt_{1}dt_{2}\delta(t_{12})c_{ss^{\prime}}(2G_{1,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2}))^{q/2}(2G_{2,ss^{\prime}}^{\alpha\beta}(t_{1},t_{2}))^{q/2},$

where the replica indices are trivial because there are no twist operators.
Similar to the previous calculation, one can show that the large-$N$ solution
of (49) for both subsystems is the same as the diagonal solution in subsystem
$a=2$ given in (36). We copy it here using the same symbol for convenience,

$\displaystyle\hat{G}_{2}=\frac{e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|t_{12}|}}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right).$ (52)

The onshell action of (49) is simply
$-\frac{I^{(0)}}{N}=2\log\text{Pf}[\hat{G}_{2}^{-1}]$ which is actually the
logarithm of Hilbert space dimension, i.e., $e^{-I^{(0)}}=Z^{n}$. Using this
action as denominator, the Rényi entropy from replica non-diagonal solution
can be written as

$\displaystyle\frac{e^{-I^{(2)}}}{Z^{n}}=e^{-(I^{(2)}-I^{(0)})}=\exp
N\Big{(}\log\text{Pf}[\hat{G}_{1}^{-1}]-\log\text{Pf}[\hat{G}_{2}^{-1}]\Big{)},$
(53)

where $\hat{G}_{1}$ and $\hat{G}_{2}$ are given by (41). However, as
discussed, (41) neglects the effect of twist operators near the boundary,
which must be included to get the correct answer. Thus we need to calculate
the Pfaffian using the full numerical solution where the effect of twist
operators is automatically included. We evaluate the Pfaffian for
$n=2,3,...,10$ in Appendix C and find that

$\displaystyle\log\text{Pf}[\hat{G}_{1}^{-1}]-\log\text{Pf}[\hat{G}_{2}^{-1}]=(1-n)\log
2.$ (54)

We expect that this result holds true for any integer $n$ because it gives the
correct answer for Brownian evolutions where the Rényi entropy is maximized at
late time.

From the on-shell actions of the replica diagonal solution (42) and replica
non-diagonal solution (53), we find that the $n$-th Rényi entropy at time $T$
is

$\displaystyle e^{-(n-1)S_{n}(T)}$ $\displaystyle=$
$\displaystyle\frac{e^{-I^{(1)}}+e^{-I^{(2)}}}{Z^{n}}=e^{-\frac{nN\mathcal{V}T}{q^{2}}}+2^{-N(n-1)},$
(55) $\displaystyle S_{n}(T)$ $\displaystyle=$
$\displaystyle\frac{1}{1-n}\log\Big{[}e^{-\frac{nN\mathcal{V}T}{q^{2}}}+2^{-N(n-1)}\Big{]}\rightarrow\begin{cases}\frac{n}{n-1}\frac{N\mathcal{V}T}{q^{2}},\quad&T\ll
T^{\ast}\\\ N\log 2,\quad&T\gg T^{\ast}\end{cases}$ (56)

which is the Page curve in the coupled Brownian SYK models and
$T^{\ast}=\frac{n-1}{n}\frac{q^{2}}{\mathcal{V}}\log 2+\mathcal{O}(1/N)$. To
convert the result to entropy per Majorana, we notice that the number of
Majorana fermion in the doubled Hilbert space of subsystem $a=1$ is $2N$, so
we have

$\displaystyle\mathcal{S}_{n}(t)$ $\displaystyle\equiv$
$\displaystyle\frac{S_{n}(T=2t)}{2N}=\begin{cases}\frac{n}{n-1}\frac{\mathcal{V}}{q^{2}}t,\quad&t\ll
t^{\ast}\\\ \frac{1}{2}\log 2,\quad&t\gg t^{\ast}\end{cases}$ (57)

where $t^{*}=T^{*}/2$ is the Page time. It is interesting to note the Page
time increases for increasing $n$, but remains finite for
$n\rightarrow\infty$. We also evaluate the Rényi entropy of the replica
diagonal and non-diagonal solution numerically. The results for $n=2,3$ are
shown in Fig. 3.

Figure 3: (a) The second Rényi entropy and (b) the third Rényi entropy per
Majorana fermion. The solid line is from the replica diagonal solution, which
grows linearly. The dashed line is from the replica non-diagonal solution that
takes over after Page time. The value of the dashed line is $\frac{1}{2}\log
2$. The dotted line plots the first line in (57) for comparison.

So far we have mainly discussed the Rényi entropy. It is an interesting
question to understand the entanglement entropy which can be obtained in
principle by analytically continue $n\rightarrow 1$ of the Rényi entropy.
There are two ways of doing the analytical continuation depending on when to
take the limit. One can keep the action off-shell and take the limit near
$n\approx 1$ to get the corresponding saddle-point equation
Almheiri:2019replica ; Penington:2019replica . In our case, since $n$ is the
dimension of various matrices $G,\Sigma,g$, etc, it is not clear how to do it.
Another way is to take the analytically continuation of the on-shell saddle-
point action. This amounts to evaluate all possible saddle-point solutions at
general $n$. We have evaluated the fully connected replica non-diagonal
solution which is responsible for the saturation of the Rényi entropy at long
times. It is reasonable to expect that there are lower symmetric saddle-point
solutions that also contribute to the action. This is also the reason that a
naive $n\rightarrow 1$ limit of (55) is ill-defined. For $n=3$ case, we
explicitly evaluate different replica non-diagonal saddle-point solutions (not
shown in the paper). We hope to generalize the calculation to arbitrary $n$,
and then analytically continue the result to $n\rightarrow 1$ in the future
work.

###  2.5 Finite time effects

So far we showed that the replica diagonal solution gives rise to linear
entropy growth while the replica non-diagonal solution gives rise to entropy
saturation at long times. When the non-diagonal saddle dominates, there were
small contributions localized near the twist operator that lead to important
effects. Here we discuss the effects of these contributions on the timescale
for entropy saturation.

In fact, there are two important times in the problem. The first is the time
at which the two saddles exchange dominance. We refer to this as the Page time
and note that, at large $N$, it is independent of $N$. The second is the time
for the entropy to reach within a few bits of its maximal value. We refer to
this time as the strong scrambling time. As we now show, this time scales like
$\log N$ at large $N$.

The key point is that in the replica non-diagonal saddle, boundary effects
contribute a term in the on-shell action of the form

$(1-n)Nb_{n}(T)e^{-\lambda T}$ (58)

where $b_{n}$ grows no faster than polynomial in $T$, and $\lambda$ is a local
interaction scale, which in our case is given by
$\lambda=2(\mathcal{J}+\mathcal{V})/q$ at large $q$ limit. After the Page
time, the Renyi entropy is thus

$S_{n}=N\log 2-Nb_{n}(T)e^{-\lambda T}.$ (59)

Hence, as $N$ goes to infinity at fixed $T$, the entropy differs from its
saturation value by an amount extensive in $N$. On the other hand, by taking
$T\sim\frac{1}{\lambda}\log N$, the correction term can be made order unity
instead of order $N$. At these times, the entropy is therefore within a few
bits of its saturation value. This form of the on-shell action follows from
the exponential decay of correlations. In the limit of large time, the twist
operator contribution is local and cannot depend on the temporal extent of the
system. Since correlations are exponentially decaying in time, it follows that
finite time effects must vanish exponentially fast, up to a polynomial
prefactor. A detailed analysis of this physics is possible in the large $q$
limit as discussed in Appendix D.

##  3 Purity from operator dynamics

Here we study the entropy dynamics from the perspective of operator dynamics,
focusing on the purity, $e^{-S_{2}}$. This approach is complementary to the
path integral approach, including offering easier access to finite $N$
corrections.

###  3.1 Operator dynamics of the coupled Brownian SYK models

As before, we consider a Hilbert space made up of $N_{1}+N_{2}$ Majorana
fermions. The Majorana fermions satisfy
$\\{\psi_{i,a},\psi_{j,b}\\}=\delta_{ab}\delta_{ij}$, and form an orthonormal
basis

$\displaystyle\Gamma_{A,B}=i^{[(|A|+|B|)/2]}2^{(|A|+|B|)/2}[\psi_{1}]_{A}[\psi_{2}]_{B},\quad[\psi_{a}]_{A}=\psi_{j_{1},a}...\psi_{j_{|A|},a},$
(60)

where $A=j_{1}...j_{|A|}$ is an ascending list of length $|A|$ and $[x]$ is
the largest integer less than or equal to $x$. In particular, we use
$\Gamma_{A,0}=\Gamma_{A,\emptyset}$ and $\Gamma_{0,A}=\Gamma_{\emptyset,A}$ to
denote the basis locating solely in subsystem $a=1,2$.

We again take $N_{1}=N_{2}$ for simplicity, so the Hilbert space dimension is
$2^{N}$. The inner product of operators is
$\text{Tr}[\mathcal{O}^{\dagger}\mathcal{O}^{\prime}]$ for any two operators
$\mathcal{O}$ and $\mathcal{O}^{\prime}$. We can decompose any operator by the
basis (60)

$\displaystyle\mathcal{O}=\sum_{A,B}c_{A,B}\Gamma_{A,B},\quad
c_{A,B}=\frac{1}{2^{N}}\text{Tr}[\mathcal{O}\Gamma_{A,B}].$ (61)

If the operator is normalized according to
$2^{-N}\text{Tr}[\mathcal{O}^{\dagger}\mathcal{O}]=1$, implying
$\sum_{A,B}|c_{A,B}|^{2}=1$, then because unitary evolution preserves the
normalization, we can interpret $|c_{A,B}|^{2}$ as the probability of finding
the operator $\mathcal{O}$ in basis operator $\Gamma_{A,B}$. Since the
disorder averaged theory has an emergent $SO(N)\times SO(N)$ symmetry for each
of the subsystem, one expects the dynamics depends only on the length of
operators not on the specific list $A$. So we define the probability
distribution of an operator in subsystem $a=1,2$ with length $m,m^{\prime}$ to
be

$\displaystyle
p_{m,m^{\prime}}(t)=\sum_{|A|=m,|B|=m^{\prime}}|c_{A,B}(t)|^{2}=\sum_{|A|=m,|B|=m^{\prime}}2^{-2N}|\text{Tr}[\mathcal{O}(t)\Gamma_{A,B}]|^{2}.$
(62)

We are interested in finding the master equation governing the time
development of this probability distribution.

In terms of the basis, the Hamiltonian (1) can be rewritten as

$\displaystyle H=\sum_{A,B}J_{A,B}\Gamma_{A,B}^{(q)},\quad J_{A,0}\equiv
2^{-q/2}J_{A}^{1},\quad J_{0,A}\equiv 2^{-q/2}J_{A}^{2},\quad J_{A,B}\equiv
2^{-q/2}V_{A,B},$ (63)

where the superscript $\Gamma^{(q)}$ implies the length of basis is $q$. The
short-hand summation denotes the same summation as in (1). Due to the Brownian
nature of the interactions, we can view the Hamiltonian as a random circuit.
At each time step, the evolution operator is generated by
$U(dt)=e^{-iH(t)dt}$. When we discretize the time interval by a tiny time step
$dt\ll T$, it is easy to check the average over Brownian variables takes the
following rule

$\displaystyle\overline{J_{A,B}(t)dtJ_{A^{\prime},B^{\prime}}(t^{\prime})dt}=\delta_{AA^{\prime}}\delta_{BB^{\prime}}\delta_{tt^{\prime}}\sigma_{A,B}dt,$
(64)
$\displaystyle\sigma_{A,0}=\sigma_{0,A}=\frac{q!\mathcal{J}}{2q^{2}N^{q-1}}\equiv\sigma_{0},\quad\sigma_{A,B}=\frac{(q/2)!^{2}\mathcal{V}}{q^{2}N^{q-1}}\equiv\sigma_{1}.$
(65)

In the following, the over line denoting the disorder average is omitted for
notational simplicity.

For a generic operator $\mathcal{O}$, the infinitesimal evolution
$\mathcal{O}(t+dt)=U^{\dagger}(dt)\mathcal{O}(t)U(dt)$ is

$\displaystyle\mathcal{O}(t+dt)=e^{iHdt}\mathcal{O}(t)e^{-iHdt}$ (66)
$\displaystyle=$
$\displaystyle\mathcal{O}(t)+i[H,\mathcal{O}(t)]dt-\frac{1}{2}\\{H^{2}dt^{2},\mathcal{O}(t)\\}+H\mathcal{O}(t)Hdt^{2}$
(67) $\displaystyle=$
$\displaystyle\left(1-(2C_{N}^{q}\sigma_{0}+(C_{N}^{q/2})^{2}\sigma_{1})dt\right)\mathcal{O}(t)+i[H,\mathcal{O}(t)]dt+\sum_{A,B}\sigma_{A,B}dt\Gamma_{A,B}^{(q)}\mathcal{O}(t)\Gamma_{A,B}^{(q)}.$
(68)

where $C_{n}^{m}\equiv\frac{n!}{m!(n-m)!}$ denotes the number of
$m$-combinations of set with $n$ elements. In the second line, we expand the
exponential function in the first line and keep up to the second order in
$dt$. In the third line we have performed the disorder average of the $dt^{2}$
terms by using (64), and due to the Ito calculus, they become linear in $dt$.
Assuming the operator $\mathcal{O}$ is Hermitian, the distribution at time
$t+dt$ via the definition (62) is given by

$\displaystyle p_{m,m^{\prime}}(t+dt)$ $\displaystyle=$ $\displaystyle
2^{-2N}\sum_{|A|=m,|B|=m^{\prime}}|\text{Tr}[O(t+dt)\Gamma_{A,B}]|^{2}$ (70)
$\displaystyle=$
$\displaystyle\left(1-2(2C_{N}^{q}\sigma_{0}+(C_{N}^{q/2})^{2}\sigma_{1})\right)p_{m,m^{\prime}}(t)$
$\displaystyle+2^{-2N}\sum_{|A|=m,|B|=m^{\prime}}\Big{(}2\text{Tr}[\mathcal{O}(t)\Gamma_{A,B}]\sum_{C,D}\sigma_{C,D}\text{Tr}[\mathcal{O}(t)\Gamma_{C,D}^{(q)}\Gamma_{A,B}\Gamma_{C,D}^{(q)}]$
$\displaystyle-\sum_{C,D}\sigma_{C,D}\text{Tr}^{2}(\mathcal{O}(t)[\Gamma_{A,B},\Gamma_{C,D}^{(q)}])\Big{)}dt,$

where the summation over $C,D$ is the same as in (63). We have performed
disorder average in the second line.

The master equation of $p_{m,m^{\prime}}(t)$ can be derived straightforwardly.
We leave the detailed derivation in Appendix E, and the result is

$\displaystyle\frac{dp_{m,m^{\prime}}(t)}{dt}$ $\displaystyle=$
$\displaystyle-4\Big{[}\sigma_{0}\sum_{k=1,\text{odd}}^{\min(q,m)}C_{N-m}^{q-k}C_{m}^{k}+\sigma_{0}\sum_{k^{\prime}=1,\text{odd}}^{\min(q,m^{\prime})}C_{N-m^{\prime}}^{q-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}$
(71)
$\displaystyle+\sigma_{1}\sum_{k=0}^{\min(q/2,m)}\sum_{k^{\prime}=0}^{\min(q/2,m^{\prime})}\frac{1-(-1)^{k+k^{\prime}}}{2}C_{N-m}^{q/2-k}C_{m}^{k}C_{N-m^{\prime}}^{q/2-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}\Big{]}p_{m,m^{\prime}}(t),$
$\displaystyle+4\Big{[}\sigma_{0}\sum_{k=1,\text{odd}}^{\min(q,m)}C_{N-(m+q-2k)}^{k}C_{m+q-2k}^{m-k}p_{m+q-2k,m^{\prime}}(t)$
$\displaystyle+\sigma_{0}\sum_{k^{\prime}=1,\text{odd}}^{\min(q,m^{\prime})}C_{N-(m^{\prime}+q-2k^{\prime})}^{k^{\prime}}C_{m^{\prime}+q-2k^{\prime}}^{m^{\prime}-k^{\prime}}p_{m,m^{\prime}+q-2k^{\prime}}(t)$
$\displaystyle+\sigma_{1}\sum_{k=0}^{\min(q/2,m)}\sum_{k^{\prime}=0}^{\min(q/2,m^{\prime})}\frac{1-(-1)^{k+k^{\prime}}}{2}C_{N-(m+q/2-2k)}^{k}C_{m+q/2-2k}^{m-k}$
$\displaystyle\times
C_{N-(m^{\prime}+q/2-2k^{\prime})}^{k^{\prime}}C_{m^{\prime}+q/2-2k^{\prime}}^{m^{\prime}-k^{\prime}}p_{m+q/2-2k,m^{\prime}+q/2-2k^{\prime}}(t)\Big{]},$

where the first two lines is the out-going rate and the rest is the in-coming
rate. It is straightforward but tedious to show that the following
distribution is a stationary solution to the master equation

$\displaystyle
p^{(\text{st})}_{m,m^{\prime}}=2^{-2N}C_{N}^{m}C_{N}^{m^{\prime}},$ (72)

which means the probability in any basis $\Gamma_{A,B}$ is the same, i.e.,
$|c_{A,B}^{(\text{st})}|^{2}=2^{-2N}$. This is consistent with the expectation
of approaching an infinite temperature state in the Brownian evolution, where
no any specific basis is preferred.

It is instructive also to check the symmetry of the master equation. First,
after the disorder average the model has an $SO(N)\times SO(N)$ symmetry, and
this is the reason that the master equation can be reduced $p_{A,B}\rightarrow
p_{m,m^{\prime}}$ depending only on the length of the basis.

Second, depending on the parity of $q/2$, i.e., the interactions between two
Brownian SYK models, the model has $Z_{2}^{f}\times Z_{2}^{f}$ for even $q/2$,
i.e., the Fermi parity is conserved separately in two subsystems, or
$Z_{2}^{f}$ for odd $q/2$, i.e., only the total Fermi parity is conserved.
This leads to the result that $p_{m,m^{\prime}}$ couples only to
$p_{m+q/2-2k,m^{\prime}+q/2-2k^{\prime}}$ in (71). If $q/2$ is even, then the
master equations for distribution $p_{m,m^{\prime}}$,
$(m,m^{\prime})\in(\text{even},\text{even})$, $(\text{even},\text{odd})$,
$(\text{odd},\text{even})$, $(\text{odd},\text{odd})$ decouple. If $q/2$ is
odd, then the master equations for distribution $p_{m,m^{\prime}}$,
$m+m^{\prime}\in\text{even},\text{odd}$ decouple.

Finally, some special operators are conserved due to the Fermi parity
symmetry. For even $q/2$, four operators $\Gamma_{0,0}$,
$\Gamma_{\\{12...N\\},0}$, $\Gamma_{0,\\{12...N\\}}$,
$\Gamma_{\\{12...N\\},\\{12...N\\}}$ are conserved, while for odd $q/2$, only
two operators $\Gamma_{0,0}$, $\Gamma_{\\{12...N\\},\\{12...N\\}}$ are
conserved.

###  3.2 Purity evolution of a pure state

Now we relate the purity evolution to the operator dynamics. The discussion in
the following is general for any initial density matrix $\rho^{\text{(in)}}$
in the Hilbert space span by $2N$ Majorana operators. The density matrix
evolves

$\displaystyle\rho(t)$ $\displaystyle=$ $\displaystyle
U(t)\rho^{\text{(in)}}U^{\dagger}(t)\equiv\sum_{A,B}c_{A,B}(t)\Gamma_{A,B},\quad
c_{A,B}(0)=c^{(\text{in})}_{A,B},$ (73)

so the reduced density matrix is

$\displaystyle\rho_{\text{1}}(t)$ $\displaystyle=$
$\displaystyle\text{Tr}_{2}[\sum_{A,B}c_{A,B}(t)\Gamma_{A,B}]=2^{N/2}\sum_{A}c_{A,0}(t)\Gamma_{A},$
(74) $\displaystyle e^{-S_{2}(t)}$ $\displaystyle=$
$\displaystyle\text{Tr}_{1}[\rho_{1}(t)^{2}]=2^{3N/2}\sum_{A}|c_{A,0}(t)|^{2}.$
(75)

Here $\Gamma_{A}\equiv i^{[|A|/2]}2^{|A|/2}[\psi_{1}]_{A}$ is the basis in
subsystem $a=1$. In the second equation, we have used the orthonormal property
of the basis $\Gamma_{A}$, namely,
$\text{Tr}_{1}[\Gamma_{A}\Gamma_{B}]=2^{N/2}\delta_{AB}$.

The evolution of the operator wavefunction $c_{A,B}$ is captured by the master
equation (71), so the purity dynamics is also dictated by the master equation.
We should notice that the density matrix is not normalized with respect to
$2^{-N}\text{Tr}[\mathcal{O}^{\dagger}\mathcal{O}]=1$. For simplicity, let us
consider pure initial state so that $\text{Tr}[\rho^{2}]=1$. We can look at
normalized operator,

$\displaystyle\mathcal{O}_{\rho}=2^{N/2}\rho=\sum_{A,B}2^{N/2}c_{A,B}(t)\Gamma_{A,B},\quad
p_{m,m^{\prime}}^{\mathcal{O}}=\sum_{|A|=m,|B|=m^{\prime}}2^{N}|c_{A,B}|^{2}.$
(76)

So the purity is

$\displaystyle e^{-S_{2}(t)}$ $\displaystyle=$ $\displaystyle
2^{N/2}\sum_{m}p_{m,0}^{\mathcal{O}}(t),$ (77)

which is $2^{N/2}$ times the probability of finding the operator in subsystem
$a=1$.

###  3.3 Setup

We consider the following state (not be confused with the TFD state considered
in the previous section), $|\infty\rangle$, that is similar to the Kourkoulou-
Maldacena state Kourkoulou:2017pure

$\displaystyle(\psi_{2j-1,a}+i\psi_{2j,a})|\infty\rangle=0,\quad\langle\infty|(\psi_{2j-1,a}-i\psi_{2j,a})=0,\quad\forall
a=1,2,\quad\forall j=1,...,N.$ (78)

So the initial density matrix is

$\displaystyle\rho^{(\text{in})}$ $\displaystyle=$
$\displaystyle|\infty\rangle\langle\infty|=\prod_{j=1}^{N/2}(\frac{1}{2}-i\psi_{2j-1,1}\psi_{2j,1})\prod_{j=1}^{N/2}(\frac{1}{2}-i\psi_{2j-1,2}\psi_{j,2})$
(79) $\displaystyle=$ $\displaystyle
2^{-N}\prod_{j=1}^{N/2}(1-\Gamma_{\\{2j-1,2j\\},0})(1-\Gamma_{0,\\{2j-1,2j\\}})\equiv\sum_{A,B}c^{(\text{in})}_{A,B}\Gamma_{A,B},$
(80)

where in the second line we used the basis defined in (60). The probability
distribution of the normalized density matrix operator at time zero is

$\displaystyle
p_{2m,2m^{\prime}}^{\mathcal{O}}(0)=2^{-N}C_{N/2}^{m}C_{N/2}^{m^{\prime}},$
(81)

which implies the purity is one (equivalently, the second Rényi entropy is
zero),

$\displaystyle
e^{-S_{2}(0)}=2^{N/2}\sum_{m}p_{m,0}^{\mathcal{O}}(0)=2^{-N/2}\sum_{m=0}^{N/2}C_{N/2}^{m}=1.$
(82)

This is consistent with the fact that the initial state is a product state
between the two subsystems.

###  3.4 Page curve from the master equation

The system is prepared in a pure state with the initial distribution given by
(81). Though it is not easy to solve the master equation exactly, we can
obtain the final probability distribution. Considering the symmetry of the
master equation, the final probability distribution depends on the parity of
$q/2$. If $q/2$ is even, the final distribution is

$\displaystyle
p_{m,m^{\prime}}=\begin{cases}\frac{1-2^{2-N}}{2^{2N-2}-4}C_{N}^{m}C_{N}^{m^{\prime}},&\\{m,m^{\prime}\in
2Z\\}\cap\\{m,m^{\prime}\neq 0,N\\}\\\
2^{-N},&(m,m^{\prime})=(0,0),(0,N),(N,0),(N,N)\\\ 0,&m,m^{\prime}\notin
2Z\end{cases}.$ (83)

which leads to the purity

$\displaystyle e^{-S_{2}}=\frac{2^{N/2}}{2^{N-2}+1},\qquad
S_{2}=\frac{N}{2}\log 2-\log 4+O(2^{-N}).$ (84)

The deficit of $\log 4$ is due to $Z_{2}^{f}\times Z_{2}^{f}$ symmetry, since
the Hilbert space has four decoupled sectors. Because the initial state (79)
is prepared in the (even, even) sector, its maximal entropy is given by
$\frac{N}{2}\log 2-\log 4$.

On the other hand, if $q/2$ is odd, the final distribution is

$\displaystyle
p_{m,m^{\prime}}=\begin{cases}\frac{1-2^{1-N}}{2^{2N-1}-2}C_{N}^{m}C_{N}^{m^{\prime}},&\\{m+m^{\prime}\in
2Z\\}\cap\\{m,m^{\prime}\neq 0,N\\}\\\ 2^{-N},&(m,m^{\prime})=(0,0),(N,N)\\\
0,&m+m^{\prime}\notin 2Z\end{cases}.$ (85)

which leads to the purity

$\displaystyle e^{-S_{2}}=\frac{2^{N/2}}{2^{N-1}+1},\qquad
S_{2}=\frac{N}{2}\log 2-\log 2+O(2^{-N}).$ (86)

The shortage of $\log 2$ is due to the $Z_{2}^{f}$ symmetry. The Hilbert space
has two decoupled sectors, leading to a deficit of $-\log 2$.

Figure 4: (a) The Page curve of the second Rényi entropy per Majorana fermion
$\mathcal{S}_{2}$ from the master equation. We start from the initial
distribution given in (81). We choose $N=66$, $q=4$, and
$\mathcal{V}/\mathcal{J}=0.2$. The dotted line is the first line in (93), and
the dashed line is $S_{2}^{th}=\frac{1}{2}\log 2-\frac{1}{N}\log 4$ where the
finite $N$ correction is included. (b) A log linear plot of
$S^{th}_{2}-S_{2}(t)$ after the Page time. It indicates an exponential
$e^{-\lambda t}$ behavior, consistent with the $\log N$ scrambling time
discussed in Sec. 2.5.

Thus we expect that under the Hamiltonian (63) evolution, the second Rényi
entropy will increase from zero to almost the largest value $\frac{N}{2}\log
2$. We can get the increase rate, which is the outgoing rate from
$P(t)=2^{N/2}\sum_{m}p_{m,0}(t)$, i.e., the second line in (71). The initial
outgoing rate is

$\displaystyle\frac{dP(t)}{dt}\Big{|}_{t=0}$ $\displaystyle=$ $\displaystyle
2^{N/2}\sum_{m=0}^{N}\frac{dp_{m,0}}{dt}\Big{|}_{t=0}=-2^{N/2}\sum_{m=0}^{N}4\sigma_{1}C_{N}^{q/2}\sum_{k=1,\text{odd}}^{\min(q/2,m)}C_{N-m}^{q/2-k}C_{m}^{k}p_{m,0}(0),$
(87) $\displaystyle\approx$
$\displaystyle-4\sigma_{1}C_{N}^{q/2}\int_{0}^{1}ds\sum_{k=1,\text{odd}}^{q/2}\frac{(1-s)^{q/2-k}s^{k}}{(q/2-k)!k!}\sqrt{\frac{N}{\pi}}e^{-N(s-1/2)^{2}}$
(88) $\displaystyle\approx$
$\displaystyle-4\sigma_{1}C_{N}^{q/2}N^{q/2}\frac{1}{2(q/2)!}\int_{0}^{1}ds[1-(1-2s)^{q/2}]\sqrt{\frac{N}{\pi}}e^{-N(s-1/2)^{2}}$
(89) $\displaystyle\approx$ $\displaystyle-\frac{2N\mathcal{V}}{q^{2}}.$ (90)

where in the second line we take the large-$N$ limit with $s=m/N$ fixed, and
use the Gaussian distribution to approximate the Binomial distribution, i.e.,

$\displaystyle
2^{N/2}p_{m,0}^{\mathcal{O}}(0)\approx\frac{2}{N}\sqrt{\frac{N}{\pi}}e^{-N(s-1/2)^{2}}.$
(91)

Thus at $t\ll 1/\mathcal{V}$, the second Rényi entropy grows linearly,

$\displaystyle e^{-S_{2}(t)}=P(t)\approx 1-\frac{2N\mathcal{V}}{q^{2}}t\approx
e^{-\frac{2N\mathcal{V}}{q^{2}}t},\quad
S_{2}(t)=\frac{2N\mathcal{V}}{q^{2}}t.$ (92)

At late time, as we know that the second Rényi entropy will saturate at
$\frac{N}{2}\log 2$, we have thus the following Page curve of second Rényi
entropy per Majorana fermion,

$\displaystyle\mathcal{S}_{2}(t)\equiv\frac{S_{2}(t)}{N}=\begin{cases}\frac{2\mathcal{V}}{q^{2}}t,&t<t^{*}\\\
\frac{1}{2}\log 2,&t>t^{*}\end{cases},$ (93)

where $t^{*}=\frac{q^{2}}{4\mathcal{V}}\log 2$ is the Page time. Even though
we start from a different initial density matrix, we can still compare the
result from the saddle point solution. We find that two results match exactly
for second Rényi entropy per Majorana (57). This is the useful quantity since
the Hilbert dimensions are different for the two cases.

We also numerically solve the master equation starting from (81) for $q/2=2$.
The purity evolution as a function of time is plotted as the solid line in
Fig. 4, explicitly showing the Page curve in the Brownian evolution. The
dotted line and dash line are given by (93). The deviation from the linear
increasing dotted line at early time is due to finite $N$ effect. For the
dashed line, we have included the finite $N$ correction from $Z_{2}^{f}\times
Z_{2}^{f}$ symmetry. Fig. 4 shows that after the Page time, the approach to
the thermal value is given by an exponential function. This is consistent with
a $\log N$ scrambling time as discussed in Sec. 2.5, where it is explained by
the effect of twist operators.

##  4 Conclusion

We studied the Rényi entropy dynamics of coupled Brownian SYK clusters using
both path integral and operator dynamics methods. While the Page curve has
been observed in random circuit models before, we showed how the replica
diagonal and non-diagonal saddle points give rise to the entanglement
behavior. This structure is very similar to the replica wormhole scenario
obtained in holographic calculations. We also discussed the scrambling time
and Page time in Section 2.5 and Appendix D, where the crucial effect of the
twist operator was discussed. One interesting future direction is to calculate
the entanglement entropy directly, and consequently to reveal the role of
entanglement islands in more generic quantum mechanical systems.

## Acknowledgements

We would like to thank Meng Cheng and Pengfei Zhang for useful discussions.
S.-K. J. would like to acknowledge helpful discussions with Shenglong Xu in
related collaborations. This work is supported by the Simons Foundation via
the It From Qubit Collaboration.

##  Appendix A Replica non-diagonal solution

According to the ansatz (36) and the Schwinger-Dyson equation (2.2), the self-
energy is

$\displaystyle\hat{\Sigma}_{1}$ $\displaystyle=$
$\displaystyle\tilde{f}_{1}\left(\begin{array}[]{cccc}0&-\tilde{\epsilon}^{T}\\\
\tilde{\epsilon}&0\end{array}\right),\quad\tilde{f}_{1}=\Big{(}\frac{\mathcal{J}}{q}f_{1}(0)^{q-1}+\frac{\mathcal{V}}{q}f_{1}(0)^{q/2-1}f_{2}(0)^{q/2}\Big{)},$
(96) $\displaystyle\hat{\Sigma}_{2}$ $\displaystyle=$
$\displaystyle\tilde{f}_{2}\left(\begin{array}[]{cccc}0&-1\\\
1&0\end{array}\right),\quad\tilde{f}_{2}=\Big{(}\frac{\mathcal{J}}{q}f_{2}(0)^{q-1}+\frac{\mathcal{V}}{q}f_{2}(0)^{q/2-1}f_{1}(0)^{q/2}\Big{)},$
(99)

where $\tilde{f}_{i}$ is a constant. Then (12) leads to

$\displaystyle\hat{G}_{1}(\omega)=\frac{1}{\omega^{2}+\tilde{f}_{1}^{2}}\left(\begin{array}[]{cccc}i\omega&-\tilde{f}_{1}\tilde{\epsilon}^{T}\\\
\tilde{f}_{1}\tilde{\epsilon}&-i\omega\end{array}\right),\quad\hat{G}_{2}(\omega)=\frac{1}{\omega^{2}+\tilde{f}_{2}^{2}}\left(\begin{array}[]{cccc}i\omega&-\tilde{f}_{2}\\\
\tilde{f}_{2}&-i\omega\end{array}\right),$ (104)
$\displaystyle\hat{G}_{1}=\frac{e^{-\tilde{f}_{1}|t_{12}|}}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-\tilde{\epsilon}^{T}\\\
\tilde{\epsilon}&-\text{sgn}(t_{12})\end{array}\right),\hat{G}_{2}=\frac{e^{-\tilde{f}_{2}|t_{12}|}}{2}\left(\begin{array}[]{cccc}\text{sgn}(t_{12})&-1\\\
1&-\text{sgn}(t_{12})\end{array}\right).$ (109)

Comparing this to the ansatz (36), we find the replica non-diagonal solution

$\displaystyle
f_{1}(t_{12})=f_{2}(t_{12})=e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|t_{12}|}.$
(110)

##  Appendix B The solution at finite time

We solve the Schwinger-Dyson equation numerically to compare it with the
analytic solutions (30, 41). Starting from a noninteracting solution as an
input, we iterate the Schwinger-Dyson equation (12, 2.2) until the result
converges. This iteration was used in Maldacena:2016remarks  and also in
Penington:2019replica .

Figure 5: The replica diagonal (a,c) and replica non-diagonal (b,d) solution
of the Schwinger-Dyson equation for $n=2$ Rényi entropy. We plot
$G_{++}^{11}(t_{1},t_{2})$ and $G_{+-}^{12}(t_{1},t_{2})$ for replica diagonal
and non-diagonal solutions, respectively. (c,d) are the values of Green’s
function on the dashed line in (a,b). The red dot (black line) represents the
numerical (analytic) solution. We choose the parameter
$q=4,\mathcal{J}/\mathcal{V}=1,\mathcal{J}T=20$. The number of discretization
is $M=400$.

We show the results in Fig. 5. The analytic solution (30, 41) matches the
numeric solution quite well for $\mathcal{J}T=20$ as shown in Fig. 5 for both
diagonal Fig. 5 and non-diagonal Fig. 5 solutions.

##  Appendix C Numerical calculation of the saddle point solution and the
onshell action

For numerical convenience, here we adopt a different convention for the
labeling of fields. We put both Keldysh contour indices and replica indices
into the time argument $0<s<2nT$ (a similar convention is used in
Chen:2020replica ): The forward contour for the $\alpha$-th replica is
$s\in((2\alpha-2)T,(2\alpha-1)T)$, and the backward contour for the
$\alpha$-th replica is $s\in((2\alpha-1)T,2\alpha T)$. We also introduce a
sign factor to capture the forward and backward contour,

$\displaystyle f(s)=\begin{cases}i,&s\in((2\alpha-2)T,(2\alpha-1)T)\\\
-i,&s\in((2\alpha-1)T,2\alpha T)\end{cases}\qquad\alpha=1,...,n.$ (111)

In this convention we also adapt the action such that the interaction between
two clusters is local at $s$. As an illustration, the contour convention for
$n=2$ is shown in Fig. 6. In this case, as seen from the figure, a replica
diagonal solution in subsystem $a=1$ will have nonvanishing correlation
between $s\in(0,2t)$ and $s\in(6t,8t)$ and between $s\in(2t,4t)$ and
$s\in(4t,6t)$. On the other hand, a nonvanishing correlation between
$s\in(0,2t)$ and $s\in(2t,4t)$ and between $s\in(4t,6t)$ and $s\in(6t,8t)$ for
subsystem $a=1$ implies a replica non-diagonal solution. We will assume $a=2$
has a diagonal solution in the following.

Figure 6: The contour convention of $n=2$ Rényi entropy. The blue (red) solid
line represents the subsystem $a=1$ ($a=2$). The arrows indicate the direction
of the time parametrization. The black dashed lines indicate the interaction
between two subsystems.

In terms of this convention, the action is

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
ds(\frac{1}{2}\psi\partial_{s}\psi+f(s)H(s)),$ (112) $\displaystyle H(s)$
$\displaystyle=$
$\displaystyle\sum_{|A|=q,a=1,2}J_{A}^{a}(s)[\psi_{a}]_{A}+\sum_{|A|=|B|=q/2}V_{A,B}(s)[\psi_{1}]_{A}[\psi_{2}]_{B}.$
(113)

where $A=j_{1}...j_{|A|}$ denotes an ascending list of length $|A|$, and
$[\psi_{a}]_{A}\equiv
i^{|A|/2}\psi_{j_{1},a}\psi_{j_{2},a}...\psi_{j_{|A|},a}$ is a short-hand
notation for $|A|$-body interaction. The summation is over all possible such
lists from $N_{a}$ Majorana fermions.

In general, the interaction strength is random variable with possible
dependence on time $s$. The distributions of the interactions are defined by
vanishing means and the following variances,

$\displaystyle\overline{J_{A}^{a}(s)J_{A^{\prime}}^{a^{\prime}}(s^{\prime})}$
$\displaystyle=$
$\displaystyle\frac{2^{q-1}q!}{q^{2}N_{a}^{q-1}}\mathcal{J}f_{J}(s-s^{\prime})\delta_{A,A^{\prime}}\delta_{a,a^{\prime}},$
(114)
$\displaystyle\overline{V_{A,B}(s)V_{A^{\prime},B^{\prime}}(s^{\prime})}$
$\displaystyle=$
$\displaystyle\frac{2^{q}(q/2)!^{2}}{q^{2}N_{1}^{(q-1)/2}N_{2}^{(q-1)/2}}\mathcal{V}f_{V}(s-s^{\prime})\delta_{A,A^{\prime}}\delta_{B,B^{\prime}},$
(115) $\displaystyle\delta_{A,A^{\prime}}$ $\displaystyle\equiv$
$\displaystyle\delta_{j_{1},j_{1}^{\prime}}...\delta_{j_{|A|},j_{|A^{\prime}|}^{\prime}}.$
(116)

where the function $f_{J}$ and $f_{V}$ characterize the time dependence of the
variances. For Brownian random variable on the contours,
$f_{J}(s,s^{\prime})=f_{V}(s,s^{\prime})=\sum_{\alpha=0}^{n-1}\delta(|s-s^{\prime}|-2\alpha
T)+\sum_{\alpha=1}^{2n-1}\delta(s+s^{\prime}-2\alpha T)]$. And for regular SYK
model, $f_{J}(s,s^{\prime})=\mathcal{J},f_{V}(s,s^{\prime})=\mathcal{V}$.

The effective action after averaging over random variables and introducing
bilocal fields, i.e. the Green’s function
$G_{a}(s_{1},s_{2})=\frac{1}{N_{a}}\sum_{j=1}^{N_{a}}\psi_{j,a}(s)\psi_{j,a}(s^{\prime})$
and the self-energy $\Sigma_{a}(s,s^{\prime})$, reads

$\displaystyle-I$ $\displaystyle=$
$\displaystyle\sum_{a}N_{a}\Big{[}\frac{1}{2}\text{Tr}\log(G_{0,a}^{-1}-\Sigma_{a})-\frac{1}{2}\int
ds_{1}ds_{2}\Sigma_{a}(s_{1},s_{2})G_{a}(s_{1},s_{2})$ (118)
$\displaystyle+\frac{\mathcal{J}}{4q^{2}}\int
ds_{1}ds_{2}f_{J}(s_{1},s_{2})f(s_{1})f(s_{2})(2G_{a}(s_{1},s_{2}))^{q}\Big{]}$
$\displaystyle+\sqrt{N_{1}N_{2}}\frac{\mathcal{V}}{2q^{2}}\int
ds_{1}ds_{2}f_{V}(s_{1},s_{2})f(s_{1})f(s_{2})(2G_{1}(s_{1},s_{2}))^{q/2}(2G_{2}(s_{1},s_{2}))^{q/2},$

The Schwinger-Dyson equation follows from the effective action is given by

$\displaystyle\hat{G}^{-1}_{a}$ $\displaystyle=$
$\displaystyle\hat{G}_{0,a}^{-1}-\hat{\Sigma}_{a},$ (119)
$\displaystyle\Sigma_{a}(s_{1},s_{2})$ $\displaystyle=$
$\displaystyle\frac{\mathcal{J}}{q}f_{J}(s_{1},s_{2})f(s_{1})f(s_{2})(2G_{a}(s_{1},s_{2}))^{q-1}$
(120)
$\displaystyle+\sqrt{\frac{N_{\bar{a}}}{N_{a}}}\frac{\mathcal{V}}{q}f_{V}(s_{1},s_{2})f(s_{1})f(s_{2})(2G_{a}(s_{1},s_{2}))^{q/2-1}(2G_{\bar{a}}(s_{1},s_{2}))^{q/2}.$

Figure 7: The replica diagonal solution $G_{2}$ (a) and replica non-diagonal
solution $G_{1}$ (b) of Schwinger-Dyson equation for the second Rényi entropy.
(c,d) Comparison between numerical solutions and analytic solutions located on
the dashed line in (b). (d) The nonvanishing correlation sourced by the twist
operator is also given by the profile of $e^{-(\mathcal{J}+\mathcal{V})|t|/q}$
at the dashed line in (b) due to the boundary condition at $t_{1}=2t$. We
choose the parameter $\mathcal{J}T=20$, $\mathcal{V}/\mathcal{J}=0.2$, $q=4$.
Note that $T=2t$. The number of discretization is $M=400$.

Figure 8: The replica diagonal solution $G_{2}$ (a) and replica non-diagonal
solution $G_{1}$ (b) of Schwinger-Dyson equation for the third Rényi entropy.
(b) shows clearly the nonvanishing correlation is sourced by the twist
operator. We choose the parameter $\mathcal{J}T=24$,
$\mathcal{V}/\mathcal{J}=0.2$, $q=4$. The number of discretization is $M=400$.

For simplicity, we consider $N_{1}=N_{2}=N$. We numerically solve the
Schwinger-Dyson equation (119, 120) for $n=2,3$ and look for replica non-
diagonal solution. In doing so, we use
$G_{0,a}(s,s^{\prime})=\frac{1}{2}\text{sgn}(s-s^{\prime})$ for times located
at the same close time path. An illustration of the close time bath for $n=2$
is given in Fig. 6. To get the replica non-diagonal solution, we start from an
initial ansatz with small but non-zero non-diagonal correlations for subsystem
$a=1$ and a diagonal initial ansatz for subsystem $a=2$. The results for $n=2$
are shown in Fig. 7. As we discuss in above, Fig. 7 is a replica non-diagonal
solution. It is intuitive to note from the figures that the only difference
between the replica diagonal solution and the replica non-diagonal solution is
those nonvanishing correlations at $\\{4t<s_{1}<8t\\}\cap\\{0<s_{2}<4t\\}$ and
$\\{0<t_{1}<4t\\}\cap\\{4t<t_{2}<8t\\}$ sourced by the twist operators located
at $s=0,2t,6t,8t$. We also get the results for $n=3$ Rényi entropy, which are
shown in Fig. 8, and there are six twist operators.

We also numerically calculate the Pfaffian in the calculation of onshell
action and check the validity of (54). The result is plotted in Fig. 9 where
we calculate $\log\text{Pf}(G_{1}^{-1}G_{2})$ for the Rényi entropy
$n=2,...,10$. The dashed line is $(1-n)\log 2$, and we find excellent
agreement of the numerical evaluated values and (54).

Figure 9: The Pfaffian in onshell action as a function of $n$. We choose
$\mathcal{J}T=50$, $\mathcal{V}/\mathcal{J}=0.2$, $q=4$. The number of
discretization is $M=40$.

##  Appendix D Non-diagonal solutions and twist operators

In this section, we discuss the effect of the twist operators. We will focus
on the replica non-diagonal solution for the second Rény entropy $n=2$ for
simplicity, while the generalization to other Rényi entropy is
straightforward. The first equation (119) couples functions non-locally in
time domain, while the second equation (120) is local. Using the large-$q$
ansatz, $G_{a}=g_{a0}(1+\frac{g_{a}}{q})$, the Schwinger-Dyson equation
becomes,

$\displaystyle\Sigma_{a}(s_{1},s_{2})$ $\displaystyle=$
$\displaystyle-\frac{1}{q}\partial_{s_{1}}\partial_{s_{2}}g_{a0}g_{a}(s_{1},s_{2}),$
(121) $\displaystyle\Sigma_{a}(s_{1},s_{2})$ $\displaystyle=$
$\displaystyle\frac{\mathcal{J}}{q}f_{J}(s_{1},s_{2})f(s_{1})f(s_{2})(2g_{a0})^{q-1}e^{g_{a}(s_{1},s_{2})}$
(122)
$\displaystyle+\sqrt{\frac{N_{\bar{a}}}{N_{a}}}\frac{\mathcal{V}}{q}f_{V}(s_{1},s_{2})f(s_{1})f(s_{2})(2g_{a0})^{q/2-1}(2g_{\bar{a}0})^{q/2}e^{\frac{1}{2}(g_{a}(s_{1},s_{2})+g_{\bar{a}}(s_{1},s_{2}))}.$

One advantage of the large-$q$ equation of motion is that it becomes local in
time variables. For simplicity, we will assume $N_{1}=N_{2}=N$. We also assume
$g_{20}=\frac{1}{2}\text{sgn}(s_{1}-s_{2})$ when $s_{1},s_{2}$ are located at
the same close time path and zero otherwise, while
$g_{10}=\frac{1}{2}\text{sgn}(s_{1}-s_{2})$ to look for the non-diagonal
solutions.

We can solve it in the regime $0<s_{1}<T$ and $0<s_{2}<T$. For Brownian case,
the large-$q$ saddle point equation reads,

$\displaystyle\partial_{s_{1}}\partial_{s_{2}}g_{a}(s_{1},s_{2})$
$\displaystyle=$ $\displaystyle
2\mathcal{J}\delta(s_{1}-s_{2})e^{g_{a}(s_{1},s_{2})}+2\mathcal{V}\delta(s_{1}-s_{2})e^{\frac{1}{2}(g_{a}(s_{1},s_{2})+g_{\bar{a}}(s_{1},s_{2}))}.$
(123)

The equation can be solved by realizing it is
$\partial_{s_{1}}\partial_{s_{2}}g_{a}(s_{1},s_{2})=0$ when $s_{1}\neq s_{2}$,
and the $\delta$ function leads to a jump in the first derivative at
$s=s_{1}=s_{2}$ that can be solved easily. Supplementing with the boundary
condition $g_{a}(s,s)=0$ and $g_{a}(s_{1},s_{2})=g_{a}(s_{2},s_{1})$, the
solutions are

$\displaystyle
g_{1}(s_{1},s_{2})=g_{2}(s_{1},s_{2})=-(\mathcal{J}+\mathcal{V})|s_{1}-s_{2}|,\quad\\{0<s_{1}<T,0<s_{2}<T\\}.$
(124)

We can extend such calculations to the regime $0<s_{1}<2T$ and $0<s_{2}<2T$,

$\displaystyle
G_{a}(s_{1},s_{2})=\frac{\text{sgn}(s_{1}-s_{2})}{2}\begin{cases}e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|s_{1}-s_{2}|},\quad&\\{s_{1},s_{2}\\}\in\\{(0,T),(0,T)\\}\cup\\{(T,2T),(T,2T)\\}\\\
e^{-\frac{\mathcal{J}+\mathcal{V}}{q}|T-s_{1}-s_{2}|},\quad&\text{otherwise}\end{cases}.$
(125)

This is consistent with the analytic solution (41) and the numeric solution
shown in Fig. (7). More generally, the solutions in
$(\alpha-1)2T<s_{1},s_{2}<2\alpha T$, $\alpha=1,...,n$ will be the same.

To investigate the nonvanishing correlation induced by the twist operator near
the boundary, we first focus on the regime $3T<s_{1}<4T$ and $0<s_{2}<T$. In
this regime, there are two twist operators, where $T_{1}$ locates at $(4T,0)$
and $T_{2}$ locates $(3T,T)$ as also indicated by the nonvanishing
correlations in Fig. 7. So the boundary conditions are $\psi(T)=\psi(3T)$ and
$\psi(0)=-\psi(4T)$. The minus sign is because of the Fermi operator. To
simplify the notation, we shift $s_{1}\rightarrow s_{1}+3T$, so the regime is
$0<s_{1},s_{2}<T$. After the redefinition, the large-$q$ equation of motion in
this regime is

$\displaystyle\partial_{s_{1}}\partial_{s_{2}}g_{1}(s_{1},s_{2})=-2\mathcal{J}\delta(T-s_{1}-s_{2})e^{g_{1}(s_{1},s_{2})}.$
(126)

The absence of the $\mathcal{V}$ term is because the replica diagonal solution
of subsystem $a=2$ vanishes in this regime.

The solution is exponentially suppressed at large $s$, and at the large time,
i.e., $\mathcal{J}T\gg 1$, the two twist operators are separated by a large
distance $T$. So we can further assume in this regime the induced solution is
separable as follows,

$\displaystyle
G_{1}(s_{1},s_{2})=G_{T_{1}}(s_{1},s_{2})+G_{T_{2}}(s_{1},s_{2}),$ (127)

where $G_{T_{i}}$ denotes the induced solution by twist operator $T_{i}$. At
large-$q$ limit, $G_{T_{i}}\approx(1+\frac{g_{T_{i}}}{q})$, and $g_{T_{i}}$
satisfies large-$q$ equation of motion (126). But they satisfy different
boundary conditions because the two twist operators locate at different
places, i.e.,

$\displaystyle g_{T_{1}}(s_{1},0)=-(\mathcal{J}+\mathcal{V})(T-s_{1}),\quad
g_{T_{1}}(T,s_{2})=-(\mathcal{J}+\mathcal{V})s_{2},$ (128) $\displaystyle
g_{T_{2}}(s_{1},T)=-(\mathcal{J}+\mathcal{V})s_{1},\quad
g_{T_{2}}(0,s_{2})=-(\mathcal{J}+\mathcal{V})(T-s_{2}).$ (129)

Let us first look at $g_{T_{1}}$. Owing to the delta function in the right-
hand-side of (126), the solution is not differentiable at $s_{1}+s_{2}=T$, so
we assume

$\displaystyle
g_{T_{1}}(s_{1},s_{2})=\begin{cases}g_{I}(s_{1},s_{2}),&s_{1}+s_{2}\leq T\\\
g_{II}(s_{1},s_{2}),&s_{1}+s_{2}\geq T\end{cases}.$ (130)

Taking the boundary conditions into consideration, the solution has the form,

$\displaystyle
g_{I}(s_{1},s_{2})=-(\mathcal{J}+\mathcal{V})(T-s_{1})+f_{T_{1}}(T-s_{2}),\quad
g_{II}(s_{1},s_{2})=-(\mathcal{J}+\mathcal{V})s_{2}+f_{T_{1}}(s_{1}),$ (131)

and $f_{T_{1}}$ satisfies

$\displaystyle\partial_{s}f_{T_{1}}(s)=\mathcal{J}+\mathcal{V}-2\mathcal{J}e^{-(\mathcal{V}+\mathcal{J})(T-s)+f_{T_{1}}(s)},\quad
f_{T_{1}}(T)=0.$ (132)

It is not hard to solve above differential equation, which leads to the
solution,

$\displaystyle
f_{T_{1}}(s)=-(\mathcal{J}+\mathcal{V})(T-s)-\log\frac{e^{-2(\mathcal{J}+\mathcal{V})(T-s)}\mathcal{J}+\mathcal{V}}{\mathcal{J}+\mathcal{V}}.$
(133)

And consequently, one can get the correlation function
$G_{T_{1}}\approx\frac{1}{2}e^{g_{T_{1}}/q}$ induced by the twist operator
$T_{1}$,

$\displaystyle
g_{T_{1}}(s_{1},s_{2})=\begin{cases}-(\mathcal{J}+\mathcal{V})(T-s_{1}-s_{2})-\log\frac{e^{-2(\mathcal{J}+\mathcal{V})s_{2}}\mathcal{J}+\mathcal{V}}{\mathcal{J}+\mathcal{V}},&s_{1}+s_{2}\leq
T\\\
-(\mathcal{J}+\mathcal{V})(T-s_{1}-s_{2})-\log\frac{e^{-2(\mathcal{J}+\mathcal{V})(T-s_{1})}\mathcal{J}+\mathcal{V}}{\mathcal{J}+\mathcal{V}},&s_{1}+s_{2}\geq
T\end{cases}.$ (134)

To simplify the notation, we define the induced Green’s function as

$\displaystyle
G_{T}(s_{1},s_{2})=\frac{1}{2}\exp\Big{[}-\frac{1}{q}\Big{(}(\mathcal{J}+\mathcal{V})(s_{1}+s_{2})+\log\frac{e^{-2(\mathcal{J}+\mathcal{V})\min(s_{1},s_{2})}\mathcal{J}+\mathcal{V}}{\mathcal{J}+\mathcal{V}}\Big{)}\Big{]}.$
(135)

Intuitively, this function is exponentially suppressed away from
$s_{1}=s_{2}=0$ where the twist operator supposed to be located. So the
solution at the regime is

$\displaystyle G_{1}(s_{1},s_{2})=G_{T}(T-s_{1},s_{2})+G_{T}(s_{1},T-s_{2}).$
(136)

This is an approximate solution with accuracy $O(e^{-\mathcal{J}T})$. The
solutions in other regimes can be obtained in the same way, so we do not have
to detail the calculation.

Now we can discuss the effect of these twist operators to the onshell action.
We mainly discuss the second Rényi entropy, but the results can be extended to
$n$-th Rényi entropy by small modifications. Taking into account the twist
operators, the replica non-diagonal solution is
$G_{1}=G_{2}+G_{T_{1}}+G_{T_{2}}$, where $G_{2}$ denotes the diagonal solution
in subsystem $a=2$ (125), and $G_{T_{1}}$ and $G_{T_{2}}$ are induced
solutions by the twist operators $T_{1}$ and $T_{2}$, respectively. Notice
$G_{2}$ and $G_{T_{i}}$ have different domain of support. The onshell action
is

$\displaystyle\log\frac{e^{-I^{(2)}}}{Z^{2}}$ $\displaystyle=$
$\displaystyle\frac{N}{2}\text{Tr}(\log(G_{1}^{-1}G_{2}))+N\frac{1-q}{q^{2}}\mathcal{J}\int
f(s_{1})f(s_{2})f_{\mathcal{J}}(s_{1},s_{2})(G_{T_{1}}+G_{T_{2}})^{q}$ (137)
$\displaystyle\approx$ $\displaystyle-N\log 2-\frac{N}{2}\int
ds_{1}ds_{2}\partial_{s_{1}}G_{T_{1}}(s_{1},s_{2})\partial_{s_{2}}G_{T_{2}}(s_{2},s_{1})$
(138)

where in the second line, we expand the $\text{Tr}\log$ term and keep the
lowest-order coupling between two twist operators $G_{T_{1}}$ and $G_{T_{2}}$,
because the factorized part contributes to the coarse grained entropy $-N\log
2$ Chen:2020replica , and we neglect other subleading terms in the large-$q$
limit. So including the parts from coupled induced Green’s function, the
contribution from the replica non-diagonal solution is

$\displaystyle\log\frac{e^{-I^{(2)}}}{Z^{2}}\approx-N\log
2+\frac{N}{q^{2}}(\mathcal{J}+\mathcal{V})^{2}T^{2}e^{-\frac{2(\mathcal{J}+\mathcal{V})T}{q}}.$
(139)

In getting above results, we neglect the second term in (135) which will not
change the essential exponential factor. Then the second Rényi entropy from
two saddle points reads

$\displaystyle e^{-S_{2}(T)}$ $\displaystyle=$
$\displaystyle\frac{e^{-I^{(1)}}+e^{-I^{(2)}}}{Z^{2}}=e^{-\frac{2N\mathcal{V}T}{q^{2}}}+e^{-N\log
2+\frac{N}{q^{2}}(\mathcal{J}+\mathcal{V})^{2}T^{2}e^{-\frac{2(\mathcal{J}+\mathcal{V})T}{q}}}.$
(140)

After the Page time when the replica non-diagonal saddle point dominates, the
Rényi entropy is actually not independent of time. The exponentially small
overlaps between two twist operators mean that it takes times proportional to
$\log N$ to fully scramble the information Lashkari:2011towards ;
Gharibyan:2018onset .

The large-$q$ analysis of the twist operator can be extended to the regular
SYK model. Here we calculate it at the infinite temperature for an
illustration. The solutions at diagonal part is simple, yielding the solution

$\displaystyle
G_{a}(s_{1},s_{2})=\begin{cases}\frac{\text{sgn}(s_{1}-s_{2})}{2}\Big{(}\frac{1}{\cosh(\mathcal{J}_{0}|s_{1}-s_{2}|)}\Big{)}^{2/q},\quad&\\{s_{1},s_{2}\\}\in\\{(0,T),(0,T)\\}\cup\\{(T,2T),(T,2T)\\}\\\
\frac{\text{sgn}(s_{1}-s_{2})}{2}\Big{(}\frac{1}{\cosh(\mathcal{J}_{0}|T-s_{1}-s_{2}|)}\Big{)}^{2/q},\quad&\text{otherwise}\end{cases}.$
(141)

where $\mathcal{J}_{0}=\sqrt{\mathcal{J}^{2}+\mathcal{V}^{2}}$. So let us
focus again on the regime $3T<s_{1}<4T$ and $0<s_{2}<T$. Redefining
$s_{1}\rightarrow s_{1}+3T$, the regime is $0<s_{1},s_{2}<T$. The equation of
motion now reads

$\displaystyle\partial_{s_{1}}\partial_{s_{2}}g_{1}(s_{1},s_{2})=-2\mathcal{J}^{2}e^{g_{1}(s_{1},s_{2})}.$
(142)

The absence of the $\mathcal{V}$ term is because the replica diagonal solution
of subsystem $a=2$ vanishes in this regime. A general solution to above
Liouville equation is
$g_{1}(s_{1},s_{2})=\log\frac{h_{1}^{\prime}(s_{1})h_{2}^{\prime}(s_{2})}{\mathcal{J}^{2}(h_{1}(s_{1})-h_{2}(s_{2}))}$.
We expect the solution is exponentially suppressed at large $s$, and at the
large time, i.e., $\mathcal{J}T\gg 1$, the two twist operators are separated
by a large distance $T$. So we can further assume in this regime the induced
solution is separable as follows,

$\displaystyle
G_{1}(s_{1},s_{2})=G_{T_{1}}(s_{1},s_{2})+G_{T_{2}}(s_{1},s_{2}),$ (143)

where $G_{T_{i}}$ denotes the induced solution by twist operator $T_{i}$. At
large-$q$ limit, $G_{T_{i}}\approx(1+\frac{g_{T_{i}}}{q})$, and $g_{T_{i}}$
satisfies the Liouville equation. But they satisfy different boundary
conditions because the two twist operators locate at different places, i.e.,

$\displaystyle
g_{T_{1}}(s_{1},0)=2\log\frac{1}{\cosh\mathcal{J}_{0}(T-s_{1})},\quad
g_{T_{1}}(T,s_{2})=2\log\frac{1}{\cosh\mathcal{J}_{0}s_{2}},$ (144)
$\displaystyle
g_{T_{2}}(s_{1},T)=2\log\frac{1}{\cosh\mathcal{J}_{0}s_{1}},\quad
g_{T_{2}}(0,s_{2})=2\log\frac{1}{\cosh\mathcal{J}_{0}(T-s_{2})}.$ (145)

After we take into account the boundary conditions, it is straightforward to
get the following solutions,

$\displaystyle G_{T_{1}}(s_{1},s_{2})=G_{T}(T-s_{1},s_{2}),\quad
G_{T_{2}}=G_{T}(s_{1},T-s_{2}),$ (146) $\displaystyle
G_{T}(s_{1},s_{2})=\frac{1}{2}\frac{1}{(\cosh\mathcal{J}_{0}s_{1}\cosh\mathcal{J}_{0}s_{2}+\frac{\mathcal{J}^{2}}{\mathcal{J}_{0}^{2}}\sinh\mathcal{J}_{0}s_{1}\sinh\mathcal{J}_{0}s_{2})^{2}}.$
(147)

It will interesting to explore the effect of these twist operators in more
details, which we leave as a future work.

##  Appendix E Derivation of the master equation

We derive the master equation in this section. We start from (70). Using the
properties of the Majorana basis, when $\Gamma_{C,D}^{(q)}$ and $\Gamma_{A,0}$
share even (odd) Majorana operators, it leads to a positive (negative) sign in
the following,

$\displaystyle\Gamma_{C,D}^{(q)}\Gamma_{A,B}\Gamma_{C,D}^{(q)}=\Gamma_{A,B},\quad\text{if
the sum of number of common elements in $A$, $C$ and $B$, $D$ is even}$
$\displaystyle\Gamma_{C,D}^{(q)}\Gamma_{A,B}\Gamma_{C,D}^{(q)}=-\Gamma_{A,B},\quad\text{if
the sum of number of common elements in $A$, $C$ and $B$, $D$ is odd}$

For a fixed list $A$ ($B$), if $C$ ($D$) and $A$ ($B$) have $k$ ($k^{\prime}$)
common elements, the number in the summation over $C$ is given by
$C_{N-m}^{q-k}C_{m}^{k}$
($C_{N-m^{\prime}}^{q-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}$) and
$C_{N-m}^{q/2-k}C_{m}^{k}$
($C_{N-m^{\prime}}^{q/2-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}$) for the intra
and inter subsystem interactions, respectively. So the second line in (70)
becomes

$\displaystyle
2^{-2N}\sum_{|A|=m,|B|=m^{\prime}}2\text{Tr}[\mathcal{O}(t)\Gamma_{A,B}]\sum_{C,D}\sigma_{C,D}\text{Tr}[\mathcal{O}(t)\Gamma_{C,D}^{(q)}\Gamma_{A,B}\Gamma_{C,D}^{(q)}]$
(149) $\displaystyle=$ $\displaystyle
2\Big{[}\sum_{k=0}^{\min(q,m)}(-1)^{k}C_{N-m}^{q-k}C_{m}^{k}\sigma_{0}+\sum_{k^{\prime}=0}^{\min(q,m^{\prime})}(-1)^{k^{\prime}}C_{N-m^{\prime}}^{q-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}\sigma_{0}$
$\displaystyle+\sum_{k=0}^{\min{(q/2,m)}}(-1)^{k}C_{N-m}^{q/2-k}C_{m}^{k}\sum_{k^{\prime}=0}^{\min{(q/2,m^{\prime})}}(-1)^{k^{\prime}}C_{N-m^{\prime}}^{q/2-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}\sigma_{1}\Big{]}p_{m,m^{\prime}}(t).$

One can combine it with the first line in (70) to give the outgoing rate,

$\displaystyle-2(2C_{N}^{q}\sigma_{0}+(C_{N}^{q/2})^{2}\sigma_{1})p_{m,m^{\prime}}(t)+2^{-2N}\sum_{A,B,C,D}2\text{Tr}[\mathcal{O}(t)\Gamma_{A,B}]\sigma_{C,D}\text{Tr}[\mathcal{O}(t)\Gamma_{C,D}^{q}\Gamma_{A,B}\Gamma_{C,D}^{q}]$
(150) $\displaystyle=$
$\displaystyle-4\Big{[}\sigma_{0}\sum_{k=1,\text{odd}}^{\min(q,m)}C_{N-m}^{q-k}C_{m}^{k}+\sigma_{0}\sum_{k^{\prime}=1,\text{odd}}^{\min(q,m^{\prime})}C_{N-m^{\prime}}^{q-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}$
$\displaystyle+\sigma_{1}\sum_{k=0}^{\min(q/2,m)}\sum_{k^{\prime}=0}^{\min(q/2,m^{\prime})}\frac{1-(-1)^{k+k^{\prime}}}{2}C_{N-m}^{q/2-k}C_{m}^{k}C_{N-m^{\prime}}^{q/2-k^{\prime}}C_{m^{\prime}}^{k^{\prime}}\Big{]}p_{m,m^{\prime}}(t).$

In deriving the result we have used the combinatorial identity
$\sum_{k=0}^{\min(m,q)}C_{m}^{k}C_{n}^{q-k}=C_{m+n}^{q}$.

We now calculate the third line in (70). The commutator
$[\Gamma_{A,B},\Gamma_{C,D}^{(q)}]$ vanishes unless $\Gamma_{A,B}$ and
$\Gamma_{C,D}^{(q)}$ shares odd common Majorana operators. Assuming lists $A$
and $C$ have $k$ common elements, the length of
$[\Gamma_{A,B},\Gamma_{C,D}^{(q)}]$ is $m+|C|-2k+|D|$. Here, $|C|=q,q/2$ for
intra and inter subsystem interactions, respectively. The summation over
$|A|=m,|C|$ overcounts the number of terms in $\\{E,|E|=m+|C|-2k\\}$. The
overcounting factor comes from the number of ways to decompose $m+|C|-2k$
length list into two lists $A$ and $C$ with $k$ common elements, i.e.,
$C_{N-(m+|C|-2k)}^{k}C_{m+|C|-2k}^{m-k}$. One can make similar analysis to
list $B$ and $D$. Then the third line in (70) gives rise to the incoming rate,

$\displaystyle-2^{-2N}\sum_{|A|=m,|B|=m^{\prime}}\sum_{C,D}\sigma_{C,D}\text{Tr}^{2}(\mathcal{O}(t)[\Gamma_{A,B},\Gamma_{C,D}^{(q)}])$
(152) $\displaystyle=$ $\displaystyle
4\Big{[}\sigma_{0}\sum_{k=1,\text{odd}}^{\min(q,m)}C_{m+q-2k}^{m-k}C_{N-(m+q-2k)}^{k}p_{m+q-2k,m^{\prime}}(t)$
$\displaystyle+\sigma_{0}\sum_{k^{\prime}=1,\text{odd}}^{\min(q,m^{\prime})}C_{m^{\prime}+q-2k^{\prime}}^{m^{\prime}-k^{\prime}}C_{N-(m^{\prime}+q-2k^{\prime})}^{k^{\prime}}p_{m,m^{\prime}+q-2k^{\prime}}(t)$
$\displaystyle+\sigma_{1}\sum_{k=0}^{\min(q/2,m)}\sum_{k^{\prime}=0}^{\min(q/2,m^{\prime})}\frac{1-(-1)^{k+k^{\prime}}}{2}C_{m+q/2-2k}^{m-k}C_{N-(m+q/2-2k)}^{k}$
$\displaystyle\times
C_{m^{\prime}+q/2-2k^{\prime}}^{m^{\prime}-k^{\prime}}C_{N-(m^{\prime}+q/2-2k^{\prime})}^{k^{\prime}}p_{m+q/2-2k,m^{\prime}+q/2-2k^{\prime}}(t)\Big{]}.$

## References

  * (1) S. Hawking, _Particle Creation by Black Holes_ , _Commun. Math. Phys._ 43 (1975) 199–220.
  * (2) D. N. Page, _Average entropy of a subsystem_ , _Phys. Rev. Lett._ 71 (1993) 1291–1294, [gr-qc/9305007].
  * (3) G. Penington, _Entanglement Wedge Reconstruction and the Information Paradox_ , 1905.08255.
  * (4) A. Almheiri, N. Engelhardt, D. Marolf and H. Maxfield, _The entropy of bulk quantum fields and the entanglement wedge of an evaporating black hole_ , _JHEP_ 12 (2019) 063, [1905.08762].
  * (5) G. Penington, S. H. Shenker, D. Stanford and Z. Yang, _Replica wormholes and the black hole interior_ , 1911.11977.
  * (6) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, _Replica Wormholes and the Entropy of Hawking Radiation_ , _JHEP_ 05 (2020) 013, [1911.12333].
  * (7) A. Almheiri, R. Mahajan, J. Maldacena and Y. Zhao, _The Page curve of Hawking radiation from semiclassical geometry_ , _JHEP_ 03 (2020) 149, [1908.10996].
  * (8) A. Almheiri, R. Mahajan and J. Maldacena, _Islands outside the horizon_ , 1910.11077.
  * (9) A. Almheiri, R. Mahajan and J. E. Santos, _Entanglement islands in higher dimensions_ , _SciPost Phys._ 9 (2020) 001, [1911.09666].
  * (10) M. Rozali, J. Sully, M. Van Raamsdonk, C. Waddell and D. Wakeham, _Information radiation in BCFT models of black holes_ , _JHEP_ 05 (2020) 004, [1910.12836].
  * (11) F. F. Gautason, L. Schneiderbauer, W. Sybesma and L. Thorlacius, _Page Curve for an Evaporating Black Hole_ , _JHEP_ 05 (2020) 091, [2004.00598].
  * (12) K. Hashimoto, N. Iizuka and Y. Matsuo, _Islands in Schwarzschild black holes_ , _JHEP_ 06 (2020) 085, [2004.05863].
  * (13) T. Hartman, E. Shaghoulian and A. Strominger, _Islands in Asymptotically Flat 2D Gravity_ , _JHEP_ 07 (2020) 022, [2004.13857].
  * (14) T. J. Hollowood and S. P. Kumar, _Islands and Page Curves for Evaporating Black Holes in JT Gravity_ , _JHEP_ 08 (2020) 094, [2004.14944].
  * (15) C. Krishnan, V. Patil and J. Pereira, _Page Curve and the Information Paradox in Flat Space_ , 2005.02993.
  * (16) H. Z. Chen, Z. Fisher, J. Hernandez, R. C. Myers and S.-M. Ruan, _Evaporating Black Holes Coupled to a Thermal Bath_ , 2007.11658.
  * (17) Y. Chen, V. Gorbenko and J. Maldacena, _Bra-ket wormholes in gravitationally prepared states_ , 2007.16091.
  * (18) T. Hartman, Y. Jiang and E. Shaghoulian, _Islands in cosmology_ , 2008.01022.
  * (19) T. Anegawa and N. Iizuka, _Notes on islands in asymptotically flat 2d dilaton black holes_ , _JHEP_ 07 (2020) 036, [2004.01601].
  * (20) C. Akers and G. Penington, _Leading order corrections to the quantum extremal surface prescription_ , 2008.03319.
  * (21) V. Balasubramanian, A. Kar and T. Ugajin, _Islands in de Sitter space_ , 2008.05275.
  * (22) V. Balasubramanian, A. Kar and T. Ugajin, _Entanglement between two disjoint universes_ , 2008.05274.
  * (23) Y. Ling, Y. Liu and Z.-Y. Xian, _Island in Charged Black Holes_ , 2010.00037.
  * (24) A. Bhattacharya, A. Chanda, S. Maulik, C. Northe and S. Roy, _Topological shadows and complexity of islands in multiboundary wormholes_ , 2010.04134.
  * (25) D. Marolf and H. Maxfield, _Observations of Hawking radiation: the Page curve and baby universes_ , 2010.06602.
  * (26) O. C. Dahlsten, R. Oliveira and M. B. Plenio, _The emergence of typical entanglement in two-party random processes_ , _Journal of Physics A: Mathematical and Theoretical_ 40 (2007) 8081.
  * (27) M. Žnidarič, _Exact convergence times for generation of random bipartite entanglement_ , _Phys. Rev. A_ 78 (Sep, 2008) 032324.
  * (28) N. Lashkari, D. Stanford, M. Hastings, T. Osborne and P. Hayden, _Towards the Fast Scrambling Conjecture_ , _JHEP_ 04 (2013) 022, [1111.6580].
  * (29) A. Nahum, J. Ruhman, S. Vijay and J. Haah, _Quantum Entanglement Growth Under Random Unitary Dynamics_ , _Phys. Rev. X_ 7 (2017) 031016, [1608.06950].
  * (30) C. von Keyserlingk, T. Rakovszky, F. Pollmann and S. Sondhi, _Operator hydrodynamics, OTOCs, and entanglement growth in systems without conservation laws_ , _Phys. Rev. X_ 8 (2018) 021013, [1705.08910].
  * (31) A. Nahum, S. Vijay and J. Haah, _Operator spreading in random unitary circuits_ , _Phys. Rev. X_ 8 (Apr, 2018) 021014.
  * (32) T. Rakovszky, F. Pollmann and C. W. von Keyserlingk, _Diffusive hydrodynamics of out-of-time-ordered correlators with charge conservation_ , _Phys. Rev. X_ 8 (Sep, 2018) 031058.
  * (33) V. Khemani, A. Vishwanath and D. A. Huse, _Operator spreading and the emergence of dissipative hydrodynamics under unitary evolution with conservation laws_ , _Phys. Rev. X_ 8 (Sep, 2018) 031057.
  * (34) A. Chan, A. De Luca and J. T. Chalker, _Solution of a minimal model for many-body quantum chaos_ , _Phys. Rev. X_ 8 (Nov, 2018) 041019.
  * (35) A. Chan, A. De Luca and J. T. Chalker, _Spectral statistics in spatially extended chaotic quantum many-body systems_ , _Phys. Rev. Lett._ 121 (Aug, 2018) 060601.
  * (36) B. Bertini, P. Kos and T. c. v. Prosen, _Exact spectral form factor in a minimal model of many-body quantum chaos_ , _Phys. Rev. Lett._ 121 (Dec, 2018) 264101.
  * (37) T. Zhou and X. Chen, _Operator dynamics in a brownian quantum circuit_ , _Phys. Rev. E_ 99 (May, 2019) 052212.
  * (38) Q. Zhuang, T. Schuster, B. Yoshida and N. Y. Yao, _Scrambling and complexity in phase space_ , _Phys. Rev. A_ 99 (Jun, 2019) 062334.
  * (39) B. Czech, K. Larjo and M. Rozali, _Black Holes as Rubik’s Cubes_ , _JHEP_ 08 (2011) 143, [1106.5229].
  * (40) S. D. Mathur and C. J. Plumberg, _Correlations in Hawking radiation and the infall problem_ , _JHEP_ 09 (2011) 093, [1101.4899].
  * (41) K. Bradler and C. Adami, _One-shot decoupling and Page curves from a dynamical model for black hole evaporation_ , _Phys. Rev. Lett._ 116 (2016) 101301, [1505.02840].
  * (42) T. Tokusumi, A. Matsumura and Y. Nambu, _Quantum Circuit Model of Black Hole Evaporation_ , _Class. Quant. Grav._ 35 (2018) 235013, [1807.07672].
  * (43) L. Piroli, C. Sünderhauf and X.-L. Qi, _A Random Unitary Circuit Model for Black Hole Evaporation_ , _JHEP_ 04 (2020) 063, [2002.09236].
  * (44) H. Liu and S. Vardhan, _A dynamical mechanism for the Page curve from quantum chaos_ , 2002.05734.
  * (45) C. Sünderhauf, L. Piroli, X.-L. Qi, N. Schuch and J. I. Cirac, _Quantum chaos in the Brownian SYK model with large finite $N$: OTOCs and tripartite information_, _JHEP_ 11 (2019) 038, [1908.00775].
  * (46) I. Kourkoulou and J. Maldacena, _Pure states in the SYK model and nearly- $AdS_{2}$ gravity_, 1707.02325.
  * (47) T. Zhou and A. Nahum, _Emergent statistical mechanics of entanglement in random unitary circuits_ , _Phys. Rev. B_ 99 (2019) 174205, [1804.09737].
  * (48) J. Maldacena and D. Stanford, _Remarks on the Sachdev-Ye-Kitaev model_ , _Phys. Rev. D_ 94 (2016) 106002, [1604.07818].
  * (49) Y. Chen, X.-L. Qi and P. Zhang, _Replica wormhole and information retrieval in the SYK model coupled to Majorana chains_ , _JHEP_ 06 (2020) 121, [2003.13147].
  * (50) H. Gharibyan, M. Hanada, S. H. Shenker and M. Tezuka, _Onset of Random Matrix Behavior in Scrambling Systems_ , _JHEP_ 07 (2018) 124, [1803.08050].