Barrier from chaos: operator entanglement dynamics of the reduced density matrix

TL;DR

This paper investigates the dynamics of operator entanglement in different quantum systems, revealing a three-phase behavior and identifying barriers that relate to quantum chaoticity across rational CFTs, random circuits, and holographic CFTs.

Contribution

It introduces a unified framework to analyze operator entanglement dynamics and identifies distinct entanglement barriers in various models, linking them to quantum chaoticity.

Findings

Operator entanglement exhibits three phases: growth, plateau, and decay.

The duration and nature of the plateau phase vary among models.

Entanglement barriers may serve as indicators of quantum chaoticity.

Abstract

It is believed that thermalization drives the reduced density matrix of a subsystem to approach a short-range entangled operator. If the initial state is also short-range entangled, it is possible that the reduced density matrix remains low-entangled throughout thermalization; or there could exist a barrier with high operator entanglement between the initial and thermalized reduced density matrix. In this paper, we study such dynamics in three classes of models: the rational CFTs, the random unitary circuit, and the holographic CFTs, representing systems of increasing quantum chaoticity. We show that in all three classes of models, the operator entanglement (or variant of) exhibits three phases, a linear growth phase, a plateau phase, and a decay phase. The plateau phase characterized by volume-law operator entanglement corresponds to the barrier in operator entanglement. While it is present in all three models, its persistence and exit show interesting distinctions among them. The rational CFTs have the shortest plateau phase, followed by the slowest decay phase; the holographic CFTs mark the opposite end, i.e. having the longest plateau phase followed by a discontinuous drop; and the random unitary circuit shows the intermediate behavior. We discuss the mechanisms underlying these behaviors in operator entanglement barriers, whose persistence might serve as another measure for quantum chaoticity.

Equations

$$|U⟩=\frac{1}{\sqrt{\text{tr}(U{U}^{†})}}\sum _{ij}{U}_{ij}|i⟩|j⟩$$

$${\rho }_{AB}(t)={\text{Tr}}_{{\mathcal{H}}_C}|\psi (t)⟩⟨\psi (t)|$$

$$|{\rho }_{AB}(t)⟩=\sum _{ij}({\rho }_{AB}(t){)}_{ij}|i⟩|j⟩$$

$${\Vert {\rho }_{AB}(t)\Vert }^2=⟨{\rho }_{AB}(t)|{\rho }_{AB}(t)⟩=\text{tr}({\rho }_{AB}^2(t))$$

$$S_n^{\mathrm{op}}\left(A,{\rho }_{AB}(t)\right)=\frac{1}{1-n}\log {\text{tr}}_{{\mathcal{H}}_A\otimes {\mathcal{H}}_A}{\left[{\text{tr}}_{{\mathcal{H}}_B\otimes {\mathcal{H}}_B}(|{\rho }_{AB}(t)⟩⟨{\rho }_{AB}(t)|/{\Vert {\rho }_{AB}(t)\Vert }^2)\right]}^n$$

$$S_n^{\mathrm{op}}(A,{\rho }_{AB}^{\alpha }(t))$$

$${\mathcal{T}}_3^{-1}=(\alpha +1,\cdots ,3\alpha )(3\alpha +1,\cdots ,5\alpha )\cdots \left((2n-1)\alpha +1,\cdots ,2n\alpha ,1,\cdots ,\alpha \right)$$

$${\mathcal{T}}_2={\left(\alpha +1,3\alpha +1,\cdots ,(2n-1)\alpha +1\right)}^{-1}\left(1,2\alpha +1,\cdots ,(2n-2)\alpha +1\right)$$

$${\mathcal{T}}_1=(123\dots 2\alpha )(2\alpha +1\cdots 4\alpha )\cdots \left((2n-2)\alpha +1\cdots ,2n\alpha \right)$$

$$h_1=h_3=n{\mathfrak{h}}_{2\alpha }{h}_2=2{\mathfrak{h}}_n$$

$$|\psi ⟩=e^{-H\frac{\beta }{4}}|B⟩$$

$$|\psi (\tau )⟩=e^{-H\tau -H\frac{\beta }{4}}|B⟩,⟨\psi (\tau )|=⟨B|e^{H\tau -H\frac{\beta }{4}}$$

$$e^{-(n-1)S_n}=\frac{⟨{\mathcal{T}}_1(w_1,{\bar{w}}_1){\mathcal{T}}_2(w_2,{\bar{w}}_2){\mathcal{T}}_3(w_3,{\bar{w}}_3){⟩}_{\text{strip}}}{⟨{\mathcal{T}}_1(w_1,{\bar{w}}_1){\mathcal{T}}_1(w_3,{\bar{w}}_3){⟩}_{\text{strip}}}$$

$$w_1=i\tau +{\ell }_B,w_2=i\tau ,w_3=i\tau -{\ell }_A$$

$$z=i\exp \left[-\frac{2\pi }{\beta }w\right]$$

$$e^{-(n-1)S_n}$$

$$\sim {\left(\frac{\beta }{2\pi }\right)}^{2h_2}\frac{⟨{\mathcal{T}}_1(z_1){\mathcal{T}}_2(z_2){\mathcal{T}}_3(z_3){\mathcal{T}}_1(z_4){\mathcal{T}}_2(z_5){\mathcal{T}}_3(z_6){⟩}_{\mathbb{C}}}{⟨{\mathcal{T}}_1(z_1){\mathcal{T}}_1(z_4){\mathcal{T}}_1(z_3){\mathcal{T}}_1(z_6){⟩}_{\mathbb{C}}}$$

$$z_1$$

$$z_4$$

$$\left(\left(\left((z_3\leftarrow z_6)\leftarrow z_2\right)\leftarrow z_5\right)\leftarrow z_1\right)\leftarrow z_4$$

$${\mathcal{T}}_1(z_3){\mathcal{T}}_1(z_6)\to (z_6{)}^{-2h_1}\text{mathds}1,{\mathcal{T}}_2(z_2){\mathcal{T}}_2(z_5)\to (z_5{)}^{-2h_2}\text{mathds}1,{\mathcal{T}}_3(z_1){\mathcal{T}}_3(z_4)\to (z_4{)}^{-2h_1}\text{mathds}1$$

$${\mathcal{T}}_1(z_3){\mathcal{T}}_1(z_6)\to (z_6{)}^{-2h_1}\text{mathds}1,{\mathcal{T}}_1(z_1){\mathcal{T}}_1(z_4)\to (z_4{)}^{-2h_1}\text{mathds}1$$

$$\approx$$

$$e^{-(n-1)S_n}\approx {\left(\frac{\beta }{2\pi }\right)}^{2h_2}(z_5{)}^{-2h_2}={\left(\frac{\beta }{2\pi }\right)}^{2h_2}\exp \left[-2h_2\left(\frac{2\pi }{\beta }t\right)\right]$$

$$S_n^{\mathrm{op}}(A,{\rho }_{AB}^{\alpha }(t))$$

$$\left(\left(\left((z_3\leftarrow z_2)\leftarrow z_6\right)\leftarrow z_5\right)\leftarrow z_1\right)\leftarrow z_4$$

$${\mathcal{T}}_3(z_3){\mathcal{T}}_2(z_2)\to (z_2{)}^{-h_2}{\mathcal{T}}_1(z_3),{\mathcal{T}}_1(z_3){\mathcal{T}}_3(z_6)\to (z_6{)}^{h_2-2h_1}{\mathcal{T}}_2(z_3)$$

$${\mathcal{T}}_2(z_3){\mathcal{T}}_2(z_5)\to (z_5{)}^{-2h_2}\text{mathds}1,{\mathcal{T}}_1(z_1){\mathcal{T}}_1(z_4)\to (z_4{)}^{-2h_1}\text{mathds}1$$

$$\text{6-point function}\approx (z_4{)}^{-2h_1}(z_2{)}^{-h_2}(z_6{)}^{h_2-2h_1}(z_5{)}^{-2h_2}$$

$$e^{-(n-1)S_n}\approx {\left(\frac{\beta }{2\pi }\right)}^{2h_2}{\left(\frac{z_6}{z_2{z}_5^2}\right)}^{h_2}={\left(\frac{\beta }{2\pi }\right)}^{2h_2}\exp \left(-\frac{2\pi }{\beta }{\ell }_A{h}_2\right)$$

$$S_n^{\mathrm{op}}(A,{\rho }_{AB}^{\alpha }(t))\approx \frac{4{\mathfrak{h}}_n}{n-1}\left(\ln \left(\frac{2\pi }{\beta }\right)+\frac{2\pi }{\beta }\frac{{\ell }_A}{2}\right)$$

$$\left(\left(\left((z_3\leftarrow z_2)\leftarrow z_6\right)\leftarrow z_1\right)\leftarrow z_5\right)\leftarrow z_4$$

Full text

1 Introduction

Quantum entanglement has played a central role in elucidating many recent progresses of theoretical physics, ranging from being an order parameter for topological orders topo_order_1 ; topo_order_2 to an explicit probe for bulk geometries in AdS/CFT correspondence RT ; HRT . Away from the static settings, the dynamics of quantum entanglement is also a crucial and universal ingredient in understanding systems out of equilibrium nandkishore_many_2015 ; calabrese_entanglement_2009 . It encodes the information about thermalization and entropy generation in the process of regaining equilibrium, and probes some surprising aspects of the systems such as pre-thermalizationkollar_generalized_2011 ; berges_prethermalization_2004 ; abanin_rigorous_2017 ; gring_relaxation_2012 , many body quantum chaoshosur_chaos_2016 ; srednicki_chaos_1994 , quantum scars turner_weak_2018 ; ho_periodic_2018 ; khemani_signatures_2018 , and in the context of AdS/CFT even the black hole interior behind the horizon hartman:2013 .

While the entanglement properties are usually defined and studied with respect to a state $|\psi\rangle$, one can also study the entanglement properties of an operator $U:\mathcal{H}\to\mathcal{H}$ that maps within a Hilbert space $\mathcal{H}$. Operationally we can pick a set of orthonormal basis in $\mathcal{H}$ and its dual $\overline{H}$ : $\{|i\rangle\}$ and $\{\langle j|\}$, and write $U=\sum_{ij}U_{ij}|i\rangle\langle j|$. The operator state can then be obtained by 111Strictly speaking, one needs to specify a mapping $f:\bar{H}\rightarrow H$ in order to define the operator state. Different choices of $f$ could lead to an ambiguity in the form of a unitary transformation $V$ acting on the second copy of the operator state defined. This could affect the operator entanglement. Here we are picking a particular $f$ by requiring $(\langle i|,f(\langle j|))=\delta_{ij}$. To remove such ambiguity, we could restrict the basis $\{|i\rangle\}$ and $\{\langle i|\}$ to have the form of tensor product states on A and B, in which the ambiguous unitary transformation factorize $V=V_{A}\otimes V_{B}$ and does not affect entanglement.

[TABLE]

A convenient example to familiarize the set up is to imagine the operator being the thermal density matrix $\rho_{\beta}$, and the corresponding operator state of $\rho_{\beta}^{\frac{1}{2}}$ being the thermal field double (TFD) state. The operator entanglement then refers to the entanglement of the operator state $|U\rangle$, defined with respect to a subsystem $A\subset\mathcal{H}\otimes\mathcal{H}$ prosen_operator_2007 ; bandyopadhyay_entangling_2005 ; prosen_chaos_2007 ; pizorn_operator_2009 .

In this paper, we study a particular class of dynamical operator entanglement: those associated with the reduced density matrices obtained from quenched states222There are related works recently on the operator entanglement nie_signature_2018 and negativity kudler-flam_quantum_2019 of the unitary evolution matrix, which are different from the operator entanglement of reduced density matrix studied in this work. 333Shortly before the post of this work, we became aware of a new preprint kusuki_dynamics_2019 , which contained the result of the operator entanglement in a local operator quench.. Let us lay out the set up step by step. Given a quenched state $|\psi(t)\rangle$ that starts with only short-range entanglement at $t=0$, and is defined on the total system $ABC$, we first obtain the reduced density matrix on the subsystem $AB$ by tracing out $\mathcal{H}_{C}$:

[TABLE]

For the interest of this paper, we take $C$ to be semi-infinite that surrounds the finite subsystem $AB$ like a heat-bath. The wave-functional of the corresponding (unnormalized) operator state $|\rho_{AB}(t)\rangle$ on the doubled subsystem Hilbert space $\mathcal{H}_{AB}\otimes\mathcal{H}_{AB}$ (see Fig. 1) is then given by the density matrix element $(\rho_{AB}(t))_{ij}=\langle i|\rho_{AB}(t)|j\rangle$ as choi_completely_1975 ; jamiolkowski_linear_1972 :

[TABLE]

where $\{|i\rangle\}$ is an orthonormal basis in $\mathcal{H}_{AB}$.

To normalize the state, we need the norm of the $|\rho_{AB}(t)\rangle$, which is the purity of the state

[TABLE]

It is in general less than $1$ and $t$-dependent when region AB is entangled with C.

We are interested in computing the operator Rényi entropy of the normalized operator state $|\rho_{AB}(t)\rangle/\left\lVert\rho_{AB}(t)\right\rVert$ on the subsystem $A\otimes A$:

[TABLE]

where we use $S^{\text{op}}_{n}\left(A,\rho_{AB}\right)$ to denote the operator Rényi entropy of $\rho_{AB}$, and $S_{n}\left(AB,\psi\right)$ as the state entanglement of the quenched state $|\psi\rangle$. As $|\psi(t)\rangle$ evolves in time $t$, so does $\rho_{AB}(t)$ and consequently $S_{n}^{\rm op}\left(A,\rho_{AB}(t)\right)$. Finding out the $t$-dependent dynamics of $S^{\text{op}}_{n}\left(A,\rho_{AB}(t)\right)$ is the main task of this paper.

An important motivation to study the dynamics of such operator entanglement comes from the following question that is both theoretically interesting and of potential practical significance: is it possible to simulate the density matrix efficiently in the matrix product operator (MPO) framework throughout the thermalization process? Since large operator entanglement indicates large bond dimension in a MPO (Fig. 1), the question can be quantified as weather or not to have low operator entanglement all the way to thermalization.

Various proposals pizorn_operator_2009 ; haegeman_time-dependent_2011 ; haegeman_unifying_2016 ; leviatan_quantum_2017 have been made for directly simulating the reduced density matrix to study the evolution to the thermal states. The operator entanglement, crudely speaking, measures the logarithm the bond dimension to represent the reduced density matrix by a matrix product operator (Fig. 1). In the process of quenched thermalization, the initial state is short-range entangled, and as a result the operator state is also short-range entangled with a low operator entanglement. On the other hand, if the system thermalizes, the final reduced density matrix will be close to an identity operator at high temperature. Thus while the state entanglement increases and saturates to volume-law, the operator entanglement is again very small. Assuming this feature of low operator entanglement persists in between, numerical algorithms that use variational approaches haegeman_time-dependent_2011 ; haegeman_unifying_2016 ; leviatan_quantum_2017 at the low bond dimension manifold should be enough to reproduce all the correlation functions faithfully.

The question is whether the feature of low operator entanglement persists or not. In general, the time-dependence of operator entanglement may result in intermediate “bumps”. If the values at the bumps are comparable to system size, then we should think of them as barriers in operator entanglement, i.e. barriers for efficient simulation of the density matrix.

In this paper we address this question in three classes of models, the two dimensional rational CFTs, the random unitary circuit, and the holographic CFTs. They represent systems with increasing quantum chaoticity maldacena_bound_2015 . All of the three models give a “growth-plateau-drop” pattern as shown in Fig. 1. The plateau after the initial growth is of volume-law value, and thus is a barrier that prevents efficient simulation using the MPO representation. The results are briefly summarized in the following outline of the paper.

In section 2, we study the problem in 2-dimensional CFTs. The computation was performed in Ref. dubail_entanglement_2017 under the implicit assumption that the theory is rational. We redo the computation and allow the operator to be the powers of the reduced density matrix, i.e. $|\rho^{\alpha}_{AB}\rangle,\alpha\in\mathbb{N}$. We find that the story is similar to what is raised up in Ref. asplund_entanglement_2015 : the rational CFT assumption entails additional OPE regimes among wist operators, and a quasi-particle interpretation which gives the result in Fig. 1. For irrational CFTs, the corresponding OPE regime disappears and explicit calculations become difficult.

In section 3, we perform the computation in a random unitary circuit using tools developed in Refs. zhou_emergent_2018 ; jonay_coarse-grained_2018 ; nahum_quantum_2017 . The random unitary circuit is a recently proposed quantum circuit model that is believed to capture the universal features of chaotic dynamics with local interactions. The circuit averaged entanglement in this model is mapped to a statistical mechanical problem of domain walls. At different stages of the evolution, three different types of the domain wall configurations dominate, which gives rise to the three continuous segments in Fig. 1.

The holographic CFTs occupy the maximally chaotic corner of irrational CFTs, in which explicit calculations become plausible again via AdS/CFT. In section 4, we study the dynamics of entanglement entropy in holographic CFTs for the operator state $|\rho^{1/2}_{AB}\rangle$. This coincides with the reflected entropy considered in Ref. tom:2019 , where it was found to closely relate to the entanglement wedge cross section EoP in the static case. Using the picture derived in tom:2019 , the dynamics of operator entanglement in this case is driven by the interplay between two HRT surfaces homologous to $AB$ and $A$ respectively in the bulk of the original state. The three different dominant configurations of the HRT surfaces then give rise to the three phases for the operator entanglement. In particular, the holographic plot in Fig. 1 is obtained from explicit calculations for $\text{AdS}_{3}/\text{CFT}_{2}$ in appendix A, which extends the results of section 2 into regimes inaccessible by OPE analysis.

While in all three models there exists a plateau barrier phase characterized by volume-law operator entanglement, thus denying the possibility of efficiently simulating these density matrices throughout thermalization, the length/persistence of these barriers seem to be correlated with how chaotic the systems are, suggesting it as another measure for quantum chaoticity. We elaborate this point further in the discussion section 5.

2 Two Dimensional CFTs

In this section, we investigate the operator Rényi entropy in two dimensional conformal field theories defined on an infinite line. The partition into subsystems ABC is illustrated in Fig. 2 .

2.1 Replica structure and twist operators

When computing state Rényi entropy, one can work in an orbifold theory by introducing pairs of twist operators that is equivalent to the replica structure calabrese_entanglement_2004 ; calabrese_quantum_2016 . In two dimensional CFTs the computation becomes especially tractable since the twist operators become quasi-local with known conformal dimensions. This has resulted in many progresses for computing state entanglements in 2d CFTs both in the ground state and quench scenarioscalabrese_entanglement_2009 ; calabrese_quantum_2016 .

The standard procedure to compute the state entropy is to represent the state by a Euclidean path-integral with an open boundary. The reduced density $\rho_{A}$ is then obtained by sewing along the compliment $A^{c}$. Quantities such as $\text{tr}(\rho_{A}^{n})$ then connect neighboring copies of $\rho_{A}$ by gluing the upper and lower rims of the open slits along $A$, leading to a partition function on a branched manifold of the form a staircase geometry. This is equivalent to the insertion of twist fields at the entanglement cut.

The operator Rényi entropy works in an analogous picture albeit with more complicated permutation structures. We consider the more general state $|\rho_{AB}^{\alpha}\rangle/\left\lVert\rho_{AB}^{\alpha}\right\rVert,\alpha\in\mathbb{N}$, i.e. operator state associated with the $\alpha$-th power of the reduced density matrix.

The reduced density matrix $\rho_{AB}$ is represented by a path-integral with a slit along $A\cup B$. The power of the reduced density matrix $\rho_{AB}^{\alpha}$ is given by gluing the upper and lower rims of the slits in a staggered manner out of $\alpha$ copies of $\rho_{AB}$, thus leaving two rims, one from the first copy and one from the $\alpha$-th copy. Placing its complex conjugate produces the unnormalized operator density matrix, see Fig. 3.

To compute the corresponding $n$-th operator Rényi entropy, we generalize Eq. (5):

[TABLE]

It proceeds in a few steps. The first term in Eq. (6) involves partial tracing in Fig. 3 over $B$, i.e. gluing along B within each operator density matrix. The resulting path-integral is then replicated $n$ times, whose slits along $A$ are then glued staggeringly, see Fig. 3. The second term in Eq. (6) comes from the norm of the operator state, see Fig. 3, where the trace is performed in $A\cup B$.

According to the descriptions above, one can work out the equivalent twist fields to be inserted (more details can be referred to Ref. tom:2019 ):

[TABLE]

Fig. 3 and Fig. 3 demonstrate the twist field insertions that reproduce Fig. 3 and Fig. 3 respectively. Our set up corresponds to the limit A and B being adjacent to each other, thus the two twist operators in the middle are fused via OPE to give $\mathcal{T}^{-1}_{3}\mathcal{T}^{-1}_{1}\to\mathcal{T}_{2},\;\mathcal{T}_{1}\mathcal{T}^{-1}_{1}\to\mathds{1}$. Define the scaling dimension of size $n$ cyclic permutation to be $\mathfrak{h}_{n}$, then the dimensions of the twist fields are:

[TABLE]

2.2 Boundary state set up

We set up the quench dynamics by considering the following initial state:

[TABLE]

where $|B\rangle$ is a conformal boundary state. The imaginary evolution makes the state normalizable with short-range entanglement at length of order $\beta$. We follow Cardy and Calabrese calabrese_quantum_2016 to set up the calculation for operator entanglement. This has been done for the operator state $|\rho_{AB}\rangle/\left\lVert\rho_{AB}\right\rVert$ in dubail_entanglement_2017 . We repeat the analysis of dubail_entanglement_2017 for general $\alpha\in\mathbb{N}$, and compute $S^{\rm op}_{n}(A,\rho^{\alpha}_{AB}(t))$.

We shall emphasize that the results are only valid for rational CFTs. A similar problem also occurs in the two interval entanglement asplund_entanglement_2015 , where the behaviors of the four point function distinguish the integrable and chaotic models. We will clarify this as we proceed.

For CFTs, it is easier to first wick rotate to Euclidean time $t\to-i\tau$:

[TABLE]

and analytically continue back later. Both $|\psi(\tau)\rangle$ and $\langle\psi(\tau)|$ can be constructed from Euclidean path integrals with boundaries. The procedure for computing the operator Rényi entropy discussed in subsection 2.1 simply amounts to inserting the appropriate twist operators when gluing the two path integrals for $|\psi(\tau)\rangle$ and $\langle\psi(\tau)|$. The setup is shown in Fig. 4: we have a strip geometry separated by the thermal length and three twists fields inserted for the unnormalized state $|\rho_{AB}\rangle$ and two twist fields for the normalization $\left\lVert\rho_{AB}^{\alpha}\right\rVert$ (Fig. 4).

The operator entanglement we are interested is the ratio of correlation functions

[TABLE]

where

[TABLE]

We proceed by conformally mapping to the upper half plane using

[TABLE]

The conformal boundary condition of $|B\rangle$ then means that the correlation functions satisfy the same ward-identities as holomorphic correlation functions on the full complex plane, but with doubled number of points coming from the images:

[TABLE]

where $z_{4}=\bar{z}_{1},z_{5}=\bar{z}_{1},z_{6}=\bar{z}_{3}$ are the images $z_{1},z_{2}$ and $z_{3}$. When continued to real time, their coordinates become:

[TABLE]

2.3 Evaluation via OPE analysis

In taking the high temperature limit $\beta\rightarrow 0$, the holomorphic conformal invariants of the correlation functions are pushed towards the boundary of moduli space for generic $t$. As a result the evaluation of correlation functions are controlled by particular OPE channels and the corresponding singularities, provided that they exist. In this subsection we proceed by assuming that all singularities arising from “chiral” OPEs exist, a point we shall come back to comment in the next subsection.

As $t$ evolves, the conformal invariants bounce from one boundary of moduli space to another, and the dominant OPE channels switch accordingly. As we summarize in the following, there are 4 phases of Eq. (14) with different OPE structure (see Fig. 5), they correspond to the segments of the ”growth-plateau-drop” pattern in Fig. 1, also see review in Ref. asplund_entanglement_2015

2.3.1 $0<t<\ell_{A}/2$

During this regime, the OPE structure is as follows:

[TABLE]

The fusion proceeds in the following steps:

[TABLE]

for the numerator and

[TABLE]

for the denominator. Thus we can evaluate the numerator and denominator respectively as:

[TABLE]

under which Eq. 14 evaluates to

[TABLE]

Thus

[TABLE]

It gives the linear growth with the rate same as that of the state entanglement.

2.3.2 $\ell_{A}/2<t<\ell_{B}/2$

In this regime, the OPE structure switches to:

[TABLE]

Now the fusion proceeds in the following steps:

[TABLE]

for the numerator. The denominator stays the same. The numerator changes to:

[TABLE]

Eq. 14 now evaluates to

[TABLE]

This gives the operator entanglement:

[TABLE]

It contains a time-independent volume-law term, and corresponds to the plateau barrier.

2.3.3 $\ell_{B}/2<t<\ell/2$

In this regime, the OPE structure switches to:

[TABLE]

Now the fusion proceeds in the following steps:

[TABLE]

for the numerator, while the denominator remains unchanged. So the numerator further changes to:

[TABLE]

Eq. 14 now evaluates to

[TABLE]

This gives the operator entanglement:

[TABLE]

which corresponds to the drop phase.

2.3.4 $t>\ell/2$

In this last regime, the OPE structure settles down to:

[TABLE]

The fusion process decompose into two clusters:

[TABLE]

In this phase, the denominator also undergoes a phase transition, whose fusion process now takes the form:

[TABLE]

The numerator and denominator is then evaluated to be:

[TABLE]

Eq. 14 now evaluates to

[TABLE]

This gives the short-range operator entanglement:

[TABLE]

The results for all four regimes are thus summarized as below (see Fig. 6)

[TABLE]

2.4 OPE singularities and quasi-particle picture

The 4 segments of the results are fixed by the 4 “holomorphic” fusion channels of the 6-point correlation function Eq. (14), see Fig. 5, and the corresponding singular behaviors therein. However, as was pointed out in asplund_entanglement_2015 , in generic CFTs the class of singular behaviors is only a subset of those produced by chiral OPEs as in the last subsection. The holomorphic correlation function with images inserted only captures the conformal transformation property of the boundary state; the OPE occurs in the UHP and is controlled by both the holomorphic and anti-holomorphic coordinates $(z,\bar{z})$. In particular, there are only two types of OPEs in the UHP:

[TABLE]

The first line Eq. (33) corresponds to ordinary OPE $\left(z_{1}\to z_{2},\;\bar{z}_{1}\to\bar{z}_{2}\right)$ in the interior of UHP; the second line Eq. (34) corresponds to the dominant channel when $\mathcal{O}_{1}$ approaches the boundary $z_{1}\leftrightarrow\bar{z}_{1}$, and fuses into a boundary operator (identity in this case).

One can check that of the 4 channels shown in Fig. 5, three of them consists solely of legitimate OPEs in the UHP and therefore whose singularities exist in generic CFTs:

• •

$0<t<\ell_{A}/2$: $z_{i}\leftrightarrow\bar{z}_{i},\;i=1,2,3$;

• •

$\ell_{A}/2<t<\ell_{B}/2$: $\left(z_{2}\to z_{3},\bar{z}_{2}\to\bar{z}_{3}\right),\;z_{3}\leftrightarrow\bar{z}_{3},\;z_{1}\leftrightarrow\bar{z}_{1}$;

• •

$t>\ell/2$: $\left(z_{2}\to z_{3},\bar{z}_{2}\to\bar{z}_{3}\right),\;\left(z_{3}\to z_{1},\bar{z}_{3}\to\bar{z}_{1}\right)$.

On the other hand, the fusion channel in the third phase $\left(\ell_{B}/2<t<\ell/2\right)$ of Fig. 5 cannot be decomposed into a series of legitimate OPEs in the UHP. Therefore the corresponding singularities giving rise to the linear decline phase in Fig. 6 do not exist in generic CFTs. It is only in rational CFTs that the existence of these singularities is necessary in order to satisfy the crossing symmetry of correlation function, see asplund_entanglement_2015 . We thus conclude that the linear decline behavior for the phase $\left(\ell_{B}/2<t<\ell/2\right)$ is a consequence of working with rational CFTs.

The nature of these singularities characterizing rational CFTs becomes more transparent if instead of the boundary state $|B\rangle$, we consider the thermal field double (TFD) state defined in the doubled CFT${}_{1}\times$CFT2:

[TABLE]

where $|n\rangle$ are eigenstates of $H$. The reduced density matrix is thus defined on (AB)${}_{1}\times$(AB)2, and the corresponding operator state is defined on a further double copy of (AB)${}_{1}\times$(AB)2. The stripe geometry in Fig. 4 becomes a cylinder of circumference $\beta$, with another set of twist operators inserted at the images. The operator Rényi entropy is then related to the full 6-point correlation function on the complex plane:

[TABLE]

As conjectured and argued in asplund_entanglement_2015 , the TFD state can be used to reproduce the singularities in the boundary state. In this case, the singularities in $\left(\ell_{B}/2<t<\ell/2\right)$ takes the form of two “tangled” light cone singularities. To see what this means, we compute the following two conformal invariants:

[TABLE]

The first line Eq. (37) indicates that $\left(z_{2},\bar{z}_{2}\right)$ is approaching the light cone “tip” of $1$ and $3$, i.e. simultaneously approaching the light cones of both $1$ and $3$; the second line Eq. (38) indicates that $\left(z_{6},\bar{z}_{6}\right)$ is approaching the light cone tip of $3$ and $5$. Each one corresponds to a light cone singularity discussed in asplund_entanglement_2015 , see Fig. 7. However, we cannot isolate them into two separate “clusters” as in asplund_entanglement_2015 ; they are tangled together, and this makes it impossible to factorize the 6-point function into a product of lower-point functions. In principle, such light cone limits of correlation functions in chaotic 2d CFTs, including finite $c$ corrections, can be explicitly studied via the structure of fusion matrix lightcone1 ; lightcone2 ; lightcone3 . We leave this for future studies.

The time dependence in Eq. (32) describes rational CFTs. As a consistency check, it can be exactly reproduced by the quasi-particle picture, which is believed to model the essential physics regarding entanglement evolution in integrable models. We present this in the thermal field double set up, see Fig. 8. The initial state is represented by the middle line. After the quench, there are coherent pairs of quasi-particles generated and they move to the opposite directions. As time goes on, the overlap regions on the line, which are the area that emitted quasi-particles at time $t$ that now entangled between region $A$ and $B^{\prime}$, is proportional to the operator entanglement. It gives rise to the 4 segments in Eq. (32) and the symmetric trapezoid in Fig. 6.

3 Random Unitary Circuit

In this section, we introduce the random unitary circuit model and the corresponding domain wall picture for the entanglement. Using this, we compute and interpret the change of the quench operator entanglement as the succession of three dominant domain wall configurations.

Random unitary circuit is a tensor network that models chaotic evolution with local interactions. The structure of the network is shown in Fig. 9, where the horizontal direction represents the system: a one-dimensional chain with $q$ bits of local degrees of freedom. The dynamics evolves along the vertical direction. Each block is an $q\times q$ independent Haar random unitary matrix. The network is constructed to involve only local interactions that entangle (almost maximally) nearby sites in one time step. And the locations switch by one lattice spacing in alternating steps as the system evolves.

The circuit retains minimally the features of evolution by a local Hamiltonian – the local generation and propagation of large entanglement – yet dispenses particular forms of interactions. The average over quenched randomness is what provides analytic handles on the dynamics of chaos. Recently, there have been fruitful results based on this model and its variants, about the scaling behaviors of out-of-time-order correlator nahum_operator_2017 ; von_keyserlingk_operator_2017 ; khemani_operator_2017 , entanglement growth jonay_coarse-grained_2018 ; nahum_dynamics_2017 ; nahum_quantum_2017 ; zhou_emergent_2018 ; chan_solution_2017 , hydrodynamic long tails in $U(1)$ conserved systems khemani_operator_2017 , etc. Similar non-unitary random circuit were also proposed to reconstruct the holographic space time as an arguably more explicit tensor network structure. The Ryu-Takayanagi (RT) surface of entanglement, transition of entanglement scalings, etc. were considered qi_space-time_2018 ; hayden_holographic_2016 ; vasseur_entanglement_2018 in this setting.

Among these results, the random unitary circuit enables a “domain wall” picture to understand the universal behavior of entanglement. It is motivated by the minimal cut picture of a general tensor network. In this picture, the number of bonds the minimal cut traverses through the tensor network gives an upper bound of the bond dimension to represent the states, and therefore is also an upper bound for the entanglement swingle2012entanglement . The domain wall picture is a refinement of the minimal cut by taking unitarity into account jonay_coarse-grained_2018 . They can be viewed as world-lines of objects carrying permutation data. Unitarity puts causality constraints on their trajectories and possible forms of interactions. More explicitly, computing the circuit averaged entanglement is mapped to a statistical mechanical problem of domain walls, where the minimal free energy gives the entanglement.

The objects carrying permutation data can be understood from the underlying replica structure defined by $\text{tr}(\rho_{A}^{n})$, see Fig. 9 for $n=2$. They perform permutation operations at the entanglement cut, which is analogous to the twist fields in the field theory definition of $\text{tr}(\rho_{A}^{n})$ using path integrals. In random unitary circuit, the random average yields an effective statistical mechanics model, whose solutions extend the permutations from the entanglement cut into the interior of the tensor network. In this sense, the domain wall gives a geometric picture to the entanglement of the random unitary circuit, in much the same way as the RT surface ryu_aspects_2006 in holography. In fact, domain walls of permutation in non-unitary random circuit were formulated to be the tensor network analog of the RT surface in attempts qi_space-time_2018 ; hayden_holographic_2016 to understand the geometric structure of the space time. In a different language, the formulation has also been studied in the context of machine learning you_entanglement_2018 .

3.1 The effective model

We now review the effective statistical mechanical model in the random unitary circuit.

Let us begin by looking at the interior of the tensor network of the quantities $\text{tr}(\rho_{A}^{n})$ shown in Fig. 10, for illustration we pick $n=2$. The overlapped blue and red squares are $n=2$ copies of the Haar random matrix $U_{q\times q}$ and its complex conjugate $U^{*}_{q\times q}$. Each block with $U$ (blue in Fig. 10) and $U^{*}$ (red in Fig. 10) is replicated $n$ times in the network. Random average gives rise to “contractions” between $U$ and $U^{*}$ among the copies. Each contraction can be represented as a permutation element in the $n$-th symmetric group, see Fig. 10. The permutation elements are what emerge as the dynamical degrees of freedom after random average, and we call them “spins”. The effective model is represented by the statistical mechanics of these spins444There are actually negative weights for some of the spin configurations. But if we integrate out half of the spins, all the spin configurations in leading order in $q$ will have positive weight. We can then treat the negative weights as perturbative corrections to the leading order configurations, thus regaining the statistical model interpretation.. Boundary conditions are imposed at both the top and bottom boundaries of the tensor network. They play different roles: the boundary condition on the top comes from the replica structure for computing the Rényi entropy; the boundary condition on the bottom encodes the initial state, see Fig. 9.

For the purpose of studying entanglement, it is more convenient to describe the configurations in terms of domain walls rather than spins, because the former fits better with the boundary conditions to be imposed and unitarity constraints. Different spin domains are characterized by different channels of contractions between the $n$ copies of $U$ and $U^{*}$, the domain walls can be understood as the twists of those contractions. For example, in order to evaluate the $n$-th Rényi entropy, the contractions operation on the top boundary of Fig. 9 are identities acting on the subsystem A, and permutations $(12\cdots n)$ acting on the subsystem B respectively. The domain wall is a twist between them, inserted at the entanglement cut. If the spins in the left domain are all $\sigma_{1}$ and those in the right domain $\sigma_{2}$, then the domain wall is a permutation $\tau=\sigma^{-1}_{1}\sigma_{2}$. In a general tensor network, these permutations may become other tensors when propagating from the top boundary into the interior. But in random unitary circuits, the domain walls (or permutations) are the only degrees of freedom that survive the random average. For instance the average purity $\overline{\text{tr}\left(\rho_{A}^{2}\right)}$ will only contain wandering domain walls as in Fig. 11. The Rényi entropy thus is given by the freedom energy of the domain walls subject to the boundary conditions at both the top and bottom boundaries. For the product initial state we are considering, we have free boundary condition at the bottom.

This picture can be extended to more generic entropic quantities such as the operator Rényi entropy we are interested, which corresponds to imposing more complicated boundary conditions, see Fig. 12. In this case, the top boundary can host more generic domain walls. They can decomposed into elementary domain walls, which are defined to consist only of single transpositions. The elementary domain walls are the building blocks for the underlying dynamics. For example, the domain wall $(123)$ consists of two elementary domain walls. The locality of the interaction and unitarity of the evolution restricts the allowed behaviors of elementary domain walls (see Fig. 13). When propagating downward, an elementary domain wall must remain inside the light cone and can not go back. When two elementary domain walls meet, they combine according to the group multiplication rule; and the same applies when a composite domain wall split.

As mentioned, the domain walls need to minimize the free energy while satisfying the boundary conditions imposed at the top and bottom of the circuit. An elementary domain wall will cost $s_{\rm eq}=\ln q$ amount of free energy for each discrete time step it takes. Therefore at the leading order in large $q$ expansion, the minimizing configuration is simply given by free elementary domain walls taking the shortest paths555Entropic corrections are of order $\frac{1}{\log q}$. . At higher orders in large $q$, elementary domain walls can have two types of interactions. One is the statistical interaction between two non-commutative domain walls, such as those with transpositions $(12)$ and $(13)$. It corrects the entanglement velocity in order $\mathcal{O}(1/\log q)$. For the special case of Rényi index $n=2$, this effect is absent. Another weaker interaction at $\mathcal{O}(1/q^{4}\ln q)$ order appears between commutative domain walls like $(12)$ and $(34)$. The suppression of domain wall interaction at $q=\infty$ effectively reduce the quench disorder to annealed disorder. We can thus equating $\overline{e^{-(n-1)S_{n}}}$ and $e^{-(n-1)\overline{S_{n}}}$.

3.2 Domain walls in operator entanglement

In this section, we study time evolution for the operator entanglement of the random unitary circuit. We proceed by listing the domain wall configurations that minimize the free energy at different stages of the evolution. We take $n=2$ to simplify the domain wall structure, and work at leading order in $q\to\infty$ so that the elementary domain walls can be treated independently. This brings in the reduction

[TABLE]

where the two terms corresponds have the initial domain wall configurations shown in Fig. 14.

To warm up, let us first consider the evolution of the 2nd Rényi entropy of the state, which is subtracted in Eq. (39). When $t$ is small, the two domain walls at the entanglement cut cost the same amount of energy when going down in any direction within the light cone. Due to entropy consideration jonay_coarse-grained_2018 ; zhou_emergent_2018 , macroscopically they go down vertically. The free energy increases linearly with rate $v_{E}\ln q\approx\ln q$. At time $t>\frac{\ell}{2}$, the two domain walls can meet and annihilate at the intersection points of their light cones. The free energy of this configuration is $\ell\ln q$, which will be the minimal when $t\gtrsim\frac{\ell}{2}$. This switch from the configuration 1 to configuration 2 in Fig. 14 gives rise to the linear growth and saturation behaviors of entanglement in a quenched state.

On the other hand, the first term in Eq. (39) has six elementary domain walls to begin with. Three types of (macroscopic) configurations will occur in their life time. Configuration 1 is the beginning stage, where all domain walls going vertically down666Free random wall will predominantly go vertically, with fluctuation of size $\sqrt{t}$. when $t$ is small. The growth rate is $6\ln q$. In configuration 2, the four domain walls from the boundaries of region A first meet and fuse into two domain wall. These two go down vertically together with the other two domain walls, giving rise to a growth rate $4\ln q$. In configuration 3 the middle domain walls split into two domain walls on the left and right, which then annihilate the remaining domain walls, giving the saturation value $2\ell\ln q=2S_{2}(AB)$. By equating their free energies, we can determine the respective transition time:

[TABLE]

In Fig. 15 we draw the time evolution of both terms in Eq. (39). The 1st term grows faster than the 2nd term in configuration 1 and then at the same rate in configuration 2. It eventually come to a tie with the state entanglement when reaching configuration 3. The operator entanglement is the difference between these two curves, which we draw in Fig. 15. It contains three stages of development: linear growth, plateau and linear decrease.

We can generalize these results to integer $n$ and $\alpha$ assuming the $q\rightarrow\infty$ limit. The domain wall at the top boundary becomes identical to the permutation data in the twist fields in Eq. (7). There are again two terms

[TABLE]

The first term is a generalization of Fig. 14: on the top boundary, both the cuts between A,C and B,C host $(2\alpha-1)n$ elementary domain walls, the cut between A, B host $2(n-1)$ elementary domain walls. Its dynamics also contain the 3 configurations in Fig. 14. The second term is a generalization of Fig. 14, the difference is only that the number of elementary domain walls is $(2\alpha-1)n$ rather than $2$ at the entanglement cut. From the geometry, one can see immediately that the transition time scales do not dependent on $n$ and $\alpha$, as long as the order of the configurations to appear is the same as $n=2,\alpha=1$. This can be confirmed by an explicit computation

[TABLE]

Further computation shows that the plateau value of $S_{n}^{\rm op}(A,\rho_{AB}^{\alpha}(t))$ is also $n$ and $\alpha$ independent.

We then claim that for a given integer $n\geq 2$ and $\alpha\geq 1$, the $q\rightarrow\infty$ collapses all the operator entanglement $S_{n}^{\rm op}(A,\rho_{AB}^{\alpha}(t))$ to the curve in Fig. 15.

The results above are derived for the large $q$ random unitary circuit. In a realistic system or even random circuit with finite $q$, we argue that such a generalized domain wall picture still exist, and the we can incorporate the non-universal feature into a coarse grained line tension function – the free energy of the domain wall per unit length.

We can not justify this assumption to the full extent, even when restricted to chaotic systems. Instead we motivate this assumption from the following aspects.

Microscopically, the domain walls are the twist of the permutations. These permutations, as we draw in Fig. 10, represents different ways to contract among copies of $U$ and $U^{*}$. These pairwise contraction gives a real deterministic value given in this case by formula shown in Fig. 10. There will contributions other than these contractions if we take a unitary gate in Fig. 10. But these will usually be complex numbers, can consequently dephase in time in chaotic systems. This is also the reason why the domain wall picture works so well in the random average of the random unitary circuit: the random average makes the dephasing exact.

Macroscopically, we can generally assume that the entanglement growth of a given cut can be written as a integral of a local quantity in space time. In the examples where the minimal cut works, it would be the integral of the increment of entanglement along the cut. Because of the translational invariance in space and time (true in the random average sense in the random circuit), we imagine that such minimal cut should be a line (or membrane in higher dimension), and the growth rate should depend only on the slope of this line. This motivates the line tension function $\mathcal{E}_{n}(v)$ for Rényi index $n$ and slope $v$.

Ref. jonay_coarse-grained_2018 motivates the existence of such line tension function by first assuming the entanglement growth rate $\partial_{t}S$ to be dependent only the a function of the slope $\Gamma(\partial_{x}S)$ at the leading order in coarse graining. The assumption only uses the local information of the entanglement (the derivative) which is reasonable in a system with only local interactions. The authors then deduce that the line tension function is the Legendre transform of $\Gamma(x)$. By constraints of entanglement (e.g. subadditivity imposed on $\Gamma(x)$), they further derive that $\mathcal{E}_{n}(v)$ has to satisfy several properties

[TABLE]

where $v_{B}$ is the butterfly velocity that characterizes the speed of operator spreading and many-body chaos jonay_coarse-grained_2018 . When entanglement saturates in Fig. 15, the free energy cost is $2\times\ell\mathcal{E}_{n}(v)/v$, which is minimized at slope $v=v_{B}$ rather than $v=1$, though the saturation value is still the maximal entanglement.

In the configurations we considered in Fig. 14, the domain walls either go vertically or at slope $v_{B}$ (which maximizes the line tension), we can then only use $\mathcal{E}(0)=v_{E}^{(n)}$ and $\mathcal{E}(v_{B})=v_{B}$ to estimate the transition points.

The random unitary circuit at the $n=2$ case is a known example that the line tension function works well in practice even for $q=2$. This is because after random average, the finite $q$ correction only introduces weak interactions between the domain walls in Fig. 14. In leading order, we have nahum_dynamics_2017 ; nahum_operator_2017 ; von_keyserlingk_operator_2017 ; zhou_emergent_2018

[TABLE]

where $v_{E}$ is a perturbative series containing higher $\mathcal{O}(1/q^{8})$ term due to the interactions of the domain walls.

Using these data, we can determine the transition points.

[TABLE]

We look at the symmetric case $\ell_{A}=\ell_{B}$. When $q=2$, we have $t_{1}>t_{2}$, which means that configuration 2 will be skipped. There is only one transition time from configuration 1 to configuration 3

[TABLE]

This is verified in the numerical results in Fig. 16, where $\ell=28$, $\ell_{A}=\ell_{B}=6$ and $q=2$. As we can see the plateau phase found in large $q$ is absent, instead only a peak occurs at about $t_{3}\approx 12.5$.

It therefore suggests that the persistence of the plateau regime is a signature of the strongly chaotic system, and likely will only occur with large on-site degrees freedom.

4 Holography

In this section, we look at holographic CFTs which represent systems that are maximally chaotic in terms of saturation of the chaos bound maldacena_bound_2015 . In AdS/CFT, the entanglement entropy is conveniently captured geometrically by the area of the RT/HRT surfaces RT ; HRT . To study the time evolution for operator entanglement in holographic CFTs, our strategy is to start from a bulk description for operator states, whose asymptotic boundaries are two-copied subsystems, and then apply the HRT prescription. The results are obtained using related pictures derived in Ref. tom:2019 . We first review the bulk picture for the stationary operator states, and then generalize to the dynamical case for generic $\alpha$. After that we study the operator entanglement entropy for $\alpha=1/2$, which is the case considered in tom:2019 , and turns out to be the natural case for constructing the bulk operator state. We will see that the operator entanglement of the reduced density matrix undergoes three phases in the bulk, which is partly driven by the background geometry representing the density matrix operator state.

4.1 Stationary states

We begin by introducing the prescription derived in Ref. tom:2019 using the language of Euclidean path-integral, which works explicitly for stationary states, or more generally states that are time reflection symmetric under: $t\to-t$.

The goal is to construct the bulk geometry dual to the operator state $|\rho_{AB}\rangle$, associated with the reduced density matrix:

[TABLE]

We assume that the global state $|\psi\rangle$ in the boundary CFT has a smooth bulk geometry. Furthermore, it is stationary with a reflection symmetry $t\to-t$ and thus admits a Euclidean path-integral definition. The reduced density matrix element between two field configurations $\phi_{1,2}$ on $AB$ is then given by:

[TABLE]

where $\Phi$ denotes collectively the relevant quantum fields, and $f(\Phi)$ denotes the state-creating sources for $|\psi\rangle$. The path integral is over configurations that has a branch-cut along $AB$, and the field configurations along the branch-cut are fixed to be $\phi_{1}$ and $\phi_{2}$ from above and below respectively.

Alternatively, this defines the (un-normalized) wave-functional of the operator state $|\rho_{AB}\rangle$ on the doubled subsystem $\{AB,A^{\prime}B^{\prime}\}$:

[TABLE]

In order to find the bulk dual of $|\rho_{AB}\rangle$, let us consider computing its norm:

[TABLE]

For stationary states this is given by a boundary Euclidean path integral over a boundary manifold $\partial\mathcal{M}$ doubly branched across $AB$ (see Fig. 17). Through the dictionary of AdS/CFT, this fixes the boundary conditions for the corresponding Euclidean path-integral over bulk fields and geometries, which matches the CFT result. Let us further suppose that the bulk Euclidean path integral is dominated by a particular saddle point contribution, denoted by $\left(\mathcal{M},g,\varphi\right)$ where $g,\varphi$ are respectively the metric and the collection of bulk fields:

[TABLE]

For simplicity we omit writing $\varphi$ in specifying bulk solutions from now on. On the other hand, path-integrals over such branched manifolds are related to the refined 2nd Renyi-entropy xi:2016 . For holographic theories and assume replica symmetry, the dominant saddle $\mathcal{B}$ can be constructed as in Ref. xi:2016 and takes the form of two identical “wedges” (see Fig. 17):

[TABLE]

where each wedge $\mathcal{\tilde{B}}_{1,2}$ can be obtained by solving under the original single-copied boundary condition, the back-reacted geometry upon inserting a “cosmic brane” $\gamma^{AB}_{2}$ that anchors at the boundary along the entangling surface $\partial{(AB)}$, and has the tension

[TABLE]

This creates a conical deficit of angle: $\Delta\phi=\pi$, so each wedge extends an angular range of $\pi$ around $\gamma^{AB}_{2}$. The two wedges are glued at a bulk Cauchy surface $\Sigma(h)$ with induced metric $h$:

[TABLE]

$\Sigma(h)$ also contains the cosmic brane $\gamma^{AB}_{2}$. The reflection symmetry characterizing the stationary state acts by reflecting across $\Sigma$ and interchanging the two wedges: $\mathcal{\tilde{B}}_{1}\leftrightarrow\mathcal{\tilde{B}}_{2}$. Semi-classically, we can interpret $\Sigma(h)$ as the bulk spatial geometry at $t=0$ dual to the operator state $|\rho_{AB}\rangle$.

Let us now describe the prescription operationally (see Fig. 18). To find the bulk geometry dual to the stationary operator state $|\rho_{AB}\rangle$, we start from the bulk dual of the original global state $|\psi\rangle$, insert a cosmic brane $\gamma^{AB}_{2}$ of tension $T=\frac{1}{4G_{N}}$ anchoring towards the boundary at $\partial(AB)$, let it back-react and settle down. The spatial geometry dual to $|\rho_{AB}\rangle$ is then obtained by identifying two copies of the back-reacted spatial geometries across the settled-down position of the cosmic brane $\gamma^{AB}_{2}$.

4.2 Dynamical states

We can extend the prescription to the time-dependent case, i.e. $|\psi\rangle$ is a dynamical state. In this case, the (un-normalized) density matrix elements, i.e. the wave-functional of the operator state, need to be constructed using the Schwinger-Keldysh formalism SK1 ; SK2 . At the instant $t=t_{0}$, it involves a path-integral over two Lorentzian time-folds $\mathcal{L}^{+}$ and $\mathcal{L}^{-}$, ranging from $t=(-\infty,t_{0})$ and $t=(t_{0},-\infty)$ respectively and glued at $t=t_{0}$ across $C$ (see Fig. 19) :

[TABLE]

where $f$ denotes the state-creating sources. Now let us follow the previous route and compute the norm:

[TABLE]

where $\partial\mathcal{M}$ is a boundary (Lorentzian) manifold branched over two Schwinger-Keldysh contours Mukund:2016 , i.e. four time-folds $\mathcal{L}^{\pm}_{1,2}$, and glued at $t=t_{0}$ (see Fig. 19). The gluing creates a local singularity on the branched manifold, in the form of a non-standard Rindler temperature: $T_{\text{Rindler}}=1/4\pi$.

Similar to before, we can push this path-integral into the bulk via AdS/CFT, and suppose further that there exists a dominant contribution to this from the Lorentzian saddle-point geometry $\left(\mathcal{B}(t_{0}),g\right)$ Mukund:2016 ; Skenderis:2009 :

[TABLE]

Strictly speaking, the existence of such a dominant Lorentzian saddle stands on shakier ground compared to the Euclidean counterpart. We shall not dwell on this too much in this paper other than acknowledging the potential subtlety, and proceed assuming such a saddle exists. The cosmic-brane construction has a natural generalization to the dynamical case. In particular, let us still assume that the replica symmetry on the boundary, which in the dynamical case interchanges the two Schwinger-Keldysh two-folds, remains a symmetry of the bulk saddle. There exists a co-dimension 2 bulk hypersurface $\gamma^{AB}_{2}$ of fixed points under this symmetry, which ends on the boundary at $\partial(AB)$ at $t=t_{0}$. The bulk saddle $\mathcal{B}(t_{0})$ takes the form of two identical bulk Schwinger-Keldysh two-folds, glued at a bulk Cauchy surface $\Sigma(t_{0})$ that contains $\gamma^{AB}_{2}$ (see Fig. 20):

[TABLE]

Furthermore, being a saddle point, the bulk geometry near $\gamma^{AB}_{2}$ should be locally smooth. For Lorentzian manifolds this implies that the local Rindler temperature takes the standard value $T_{\text{Rindler}}=1/2\pi$. This refers to the branched four-folds, and as a result, each of the two-folds $\tilde{\mathcal{B}}_{1,2}(t_{0})=\tilde{\mathcal{B}}^{+}_{1,2}(t_{0})\cup\tilde{\mathcal{B}}^{-}_{1,2}(t_{0})$ has a local Rindler temperature near $\gamma^{AB}_{2}$ of $T_{\text{Rindler}}=1/\pi$. This is the back-reacted bulk solution obtained by inserting a cosmic brane of tension $T_{2}=\frac{1}{8G_{N}}$ ending on $\partial(AB)$ at $t=t_{0}$, with the boundary being the Lorentzian two-fold $\mathcal{L}^{+}\cup\mathcal{L}^{-}$ that defines the density matrix in Eq. (55). The settled-down location of the cosmic brane is $\gamma^{AB}_{2}$. The bulk solution $\tilde{\mathcal{B}}_{1}(t_{0})$ obtained this way then defines the leading order description of the bulk state dual to the operator state $|\rho_{AB}(t_{0})\rangle$. We can understand the instantaneous spatial geometry $\Sigma(t_{0})$ as being obtained by identifying two copies of the back-reacted spatial geometries ending on $AB$ at $t=t_{0}$ across the cosmic brane $\gamma^{AB}_{2}$.

It is worthwhile to clarify the notions of dynamics. The operator state $|\rho_{AB}(t_{0})\rangle$ “evolves” as we change $t_{0}$. However, this evolution is not unitary, as manifested by the $t_{0}$-dependent norm $\mathcal{N}(t_{0})$. As a result, such $t_{0}$-dependence is not encoded by a single Lorentzian bulk solution; instead, $t_{0}$ serves as a “label” for the one-parameter family of bulk saddles $\left\{\tilde{\mathcal{B}}_{1}(t_{0}),\;t_{0}\in\mathbb{R}\right\}$. On the other hand, for a fixed $t_{0}$ one can evolve the state $|\rho_{AB}(t_{0},s)\rangle=U(s)|\rho_{AB}(t_{0})\rangle$ using some unitary operator $U(s)$ and study the dynamics in $s$. A canonical candidate is $U(s)=e^{is\left(K^{\psi}_{AB}+K^{\psi}_{A^{\prime}B^{\prime}}\right)}$, where $K^{\psi}_{AB}$ is the half-side operator modular Hamiltonian, and is related to the original state Hamiltonian $K^{\Psi}_{AB}=-\ln{\rho_{AB}}$ by:

[TABLE]

$K^{\psi}_{A^{\prime}B^{\prime}}$ is an identical copy of $K^{\psi}_{AB}$. Notice that $U(s)$ is not the full modular-flow, which acts as a symmetry of the operator state. In the vicinity of $\gamma^{AB}_{2}(t_{0})$, the action of $U(s)$ can be approximated geometrically by two identical rindler boosts about $\gamma^{AB}_{2}(t_{0})$, with $s$ being the rapidity tom:2018 ; xi:2018 . This effectively generates a time-like evolution in the near-$\gamma^{AB}_{2}(t_{0})$ region of the bulk operator state.

4.3 Generalization to $|\rho^{\alpha}_{AB}(t)\rangle$

One can generalize the construction discussed before, and consider operator states of the form considered in section 2:

[TABLE]

where we take $\alpha$ as an integer to begin with. By making the same assumptions (e.g. replica symmetry, existence of a dominant saddle, etc) it is easy to see that the construction for the bulk dual proceeds the same, with the only change being the tension of cosmic-brane:

[TABLE]

Next, we can analytically continue $\alpha$ away from integral values. A special case occurs at the value $\alpha=\frac{1}{2}$. In this case $T_{\alpha}=0$, i.e. the cosmic brane is tension-less, and the location of the brane is given by the HRT surface in the original geometry HRT . More importantly, since the cosmic brane does not generate time-dependent back-reaction, the “parent” bulk geometry for constructing the dual of the operator state $|\rho^{\frac{1}{2}}_{AB}(t)\rangle$ remain the same as we evolve in $t$. The only $t$-dependent ingredient is the HRT surface across which we make the identification.

4.4 Operator Entanglement of $e^{-\beta/2H}|B\rangle$

Let us now consider the holographic picture for the operator entanglement of the generalized operator state $|\rho^{\frac{1}{2}}_{AB}(t)\rangle$, where we take the global state $|\psi\rangle$ to be pure B-state $e^{-H\beta/2}|B\rangle$. In a holographic $\text{CFT}_{d}$, this is dual to half side of a two-sided external black hole with the metric:

[TABLE]

Furthermore it contains an end-of-the-world (EoW) brane cutting through the interior maldacena_eternal_2003 . The EoW brane can be understood as arising from extending the boundary condition on the CFT into the bulk BCFT . A particular type of boundary condition in the bulk is that the HRT suface can end perpendicular onto the EoW brane.

We are interested in the operator entanglement entropy:

[TABLE]

This coincides with the reflective entropy considered in Ref. tom:2019 , where for static cases it was shown to be equal to twice of the entanglement wedge cross section, the latter has been conjectured as a holographic realization of the entanglement of purification EoP . Via the HRT prescription this is obtained by finding the extremal surface $\gamma^{A}(t)$ in bulk geometry dual to the operator state. For $\alpha=\frac{1}{2}$, the bulk dual is given by the original dual (i.e. half-sided eternal black brane with the EoW brane) identified across the the HRT surface $\gamma^{AB}(t)$. This allows us to extract $S^{\rm op}_{A}(t)$ directly in the original bulk geometry by modifying the rule of finding the HRT surface $\gamma^{A}(t)$ as follows:

• •

A HRT surface $\gamma^{A}(t)$ in the original geometry remains a candidate HRT surface

• •

A candidate HRT surface can also end perpendicularly on $\gamma^{AB}(t)$.

• •

$S^{\rm op}_{A}(t)$ is given by twice the smallest area among all candidate HRT surfaces.

In other words, $S^{\rm op}_{A}(t)$ in this case equals twice of the dynamical version of the entanglement wedge cross section for A inside AB. Let us focus on the case where both $AB$ and $A$ are flat stripes, ranging between $(-\ell_{A},\ell_{B})$ and $(-\ell_{A},0)$ respectively along the spatial $y$ direction, where $0<\ell_{A}\leq\ell/2,\;\ell=\ell_{A}+\ell_{B}$.

The first step is to determine the HRT surface $\gamma^{AB}(t)$ in the parent one-sided black brane geometry, on which the HRT surfaces $\gamma^{A}(t)$ of the operator state can end. The time dependence of $\gamma^{AB}(t)$ has been studied in hartman:2013 , and exhibits two phases in the high temperature limit $\beta\ll 1$:

• •

For $v_{d}t<\ell/2$, $\gamma^{AB}(t)$ takes the form of two planes at fixed values of $y=-\ell_{A}$ and $\ell_{B}$ respectively, extending into the bulk and ending perpendicularly on the EoW brane, with the total proper area growing with time $\text{Area}\left(\gamma^{AB}(t)\right)\propto\frac{8\pi}{\beta}v_{d}t$;

• •

For $v_{d}t>\ell/2$, $\gamma^{AB}(t)$ is a connected extremal surface that extends in the $y$ direction, most of which lie very close to the bifurcating surface $\rho=0$ and has a volume-law proper area $\text{Area}\left(\gamma^{AB}(t)\right)\propto\frac{4\pi}{\beta}\ell$.

where the velocity $v_{d}=\frac{\sqrt{d}(d-2)^{1/2-1/d}}{\left[2(d-1)\right]^{1-1/d}}$ depends on the boundary CFT dimension $d$ hartman:2013 . On top of this evolution, we expect three phases for the behavior of the HRT surface $\gamma^{A}(t)$, summarized as in Fig. 21:

• •

For $v_{d}t<f\ell_{A}$ with a proportionality constant $f$ to be determined, $\gamma^{A}(t)$ is similar to each component of $\gamma^{AB}(t)$ at the same time, and takes the form of a plane at fixed $y=0$, extending into the bulk and ending perpendicularly on the EoW brane. During this phase, the entanglement entropy grows linearly in time: $S^{\rm op}_{A}(t)\propto\frac{4\pi}{\beta}v_{d}t$.

• •

For $f\ell_{A}<v_{d}t<\ell/2$, $\gamma^{A}(t)$ latches onto one of the components of $\gamma^{AB}(t)$. During this phase the entanglement entropy is approximately time-independent and exhibits volume-law scaling: $S^{\rm op}_{A}(t)\propto\frac{4\pi}{\beta}f\ell_{A}$.

• •

For $v_{d}t>\ell/2$, due to the transition of $\gamma^{AB}(t)$, $\gamma^{A}(t)$ undergoes a discontinuous transition, changing from volume-scaling to being short-ranged.

As we can see, the decline phase is absent for holographic theories, instead the plateau barrier persists until the very end of thermalization. In appendix (A), we calculate the operator entanglement $S^{\rm op}_{A}(t)$ explicitly in AdS3/CFT2, where $v_{2}=1$. After subtracting the usual UV divergent contribution, the results are summarized below, and plotted in Fig. 22:

[TABLE]

This can be compared with the rational CFT results by sending $n\to 1,\alpha\to 1/2$ in Eq. (32). In particular, we can identify explicitly that $f=1/2$, agreeing with the CFT results. We make a comment about this. The constant $f$ is determined by the HRT surface $\gamma^{A}(t)$ in the second phase, which has a proper area that scales with the volume $\sim\ell_{A}$, see Fig. 21. Naïvely an HRT surface that scales with volume-law would have its most part sticking to the bifurcating surface, like (half of) $\gamma^{AB}(t)$ in the third phase, see Fig. 21. Let us call surfaces like this the thermal type. It is easy to see that the thermal type would give the incorrect answer $f=1$.

So $\gamma^{A}(t)$ has smaller proper area than the thermal type by a factor of $f$, which is $1/2$ in AdS3/CFT2. Geometrically this can be understood as follows. The thermal type for $\gamma^{A}(t)$ would emerge if we extremize all surfaces that end perpendicularly on the full time-like hyperplane $y=-\ell_{A}$; while the actual HRT surface comes from extremizing those that end perpendicularly on $\gamma^{AB}(t)$, which is a particular Cauchy slice of the $y=-\ell_{A}$ hyperplane. The extremizing procedure yielding the thermal type can be represented by a maxmin procedure Wall_2014 : foliate the relevant spacetime region into Cauchy surfaces; minimize on each Cauchy surface; then maximize the minimal values over all Cauchy surfaces. In this picture, the actual HRT surface corresponds to a particular minimal configuration on the Cauchy surface whose boundary is $\gamma^{AB}(t)$. In the maxmin construction, it has smaller proper area than the thermal type.

As a result, $\gamma^{A}(t)$ does not stick to the bifurcating surface, and in fact passes through the horizon. One can then investigate what fraction of $\gamma^{A}(t)$ is in the interior of the black brane. In appendix A.4 we check this explicitly, and found that the entire volume-scaling part of $\gamma^{A}(t)$ is from behind the horizon, see Fig. 23. More explicitly, the portion of of $\gamma^{A}(t)$ that is behind the horizon covers the interval in $y$ coordinate:

[TABLE]

Such an extensive amount of entanglement from the black hole interior is only accessible from the operator entanglement.

5 Discussion

In this paper, we study the operator entanglement of the reduced density matrix of a quenched state in three representative systems: random unitary circuit with local interactions, rational CFT and holography. We find that in all those cases, the operator entanglement grows linearly from short-range value to an extensive value, stays constant for some time thus forming a plateau, and then drops to short-range value. These three phases has been summarized in Fig. 1, which we reproduce below in Fig. 24 for the readers’ convenience.

To answer the question raised in the introduction as motivation, our results suggest that for generic non-integrable systems such as the random unitary circuits or the holographic CFTs, there exists an operator entanglement barrier of volume-law value during the thermalization of a quenched state. Hence algorithms that use a low entanglement representation can not faithfully reproduce the reduced density matrix throughout the process. This is a natural expectation for chaotic systems. On the other hand, the rational CFTs also share the operator entanglement barrier despite the evolution Hamiltonian being integrable. The reason lies in our choice of the initial global quenched state $|\psi\rangle=e^{-\beta H/4}|B\rangle$, which is a conformal invariant boundary state perturbed by the energy density operator $T_{tt}$. Such perturbation always exists in CFTs, and it has been shown in cardy_quantum_2015 ; wen_entanglement_2018 that the reduced density matrix $\rho_{A}(t)$ under the global quench converges to the thermal density matrix $\frac{1}{Z}e^{-\beta H}$ as $t\rightarrow\infty$. 777More precisely, cardy_quantum_2015 shows that their overlap converges to 1; wen_entanglement_2018 derives the convergence in operators for geometries in which the modular Hamiltonian $K_{A}(t)=-\log(\rho_{A}(t))$ can be explicitly calculated. The fate of thermalization for the initial states we considered puts the rational CFTs in equal footing as the other two chaotic models for comparison. However, it is also known that other irrelevant perturbations of the boundary state (for example by the descendants of the stress tensor) could lead to a Generalized Gibbs ensemble (GGE) rigol_relaxation_2007 at late times cardy_quantum_2015 . The breakdown of thermalization could evade the decay phase, unless the chemical potentials are all tuned to be small. It would be interesting to study how the operator entanglement behave by quenching from the general initial states.

Strictly speaking, the explicit computations we performed in different models are for different operator entanglement quantities: in rational CFT and random unitary circuit we obtained the results for general $n$-th operator Rényi entropy, and for operator state with general $\alpha$; in holographic CFTs what we computed is the reflective entropy, i.e. operator entanglement entropy $n\to 1$ with $\alpha=\frac{1}{2}$. That these results plotted on the same figure are showing the same plateau values therefore requires some technical explanation. First of all, the time scales of the systems are calibrated so that the effective speeds of light controlling the causality of local interactions are all set to $1$. Furthermore, we have re-scaled the entanglement quantities by their corresponding $n$-dependent equilibrium Rényi entropy density. More specifically, what we plot is the ratio $S^{\rm op}_{n}(A,\rho^{\alpha}_{AB}(t))/(\ell s^{(n)}_{\text{eq}})$, where $s^{(n)}_{\text{eq}}$ is the equilibrium Rényi entropy density of the quenched state at $t\rightarrow\infty$. More importantly, one might worry that transition time could depend on the Rényi index $n$, and therefore the comparison between different models in Fig. 24 may be meaningless. We acknowledge this potential subtlety. However we point out that for integer $n$ and $\alpha$, both the CFT and random unitary circuit results are explicitly given: they are independent of $\alpha$ and $n$ only enters as a pre-factor $s_{\rm eq}^{(n)}$ (in fact, it is also $n$ independent in random circuit when $q\rightarrow\infty$.888Strictly speaking, when computing operator entanglement entropy, this is true if we take $q\rightarrow\infty$ before the replica limit. ). Therefore one can analytically continue them to $n=1$, $\alpha=half$ and compare with holographic CFTs, and this should agree with Fig. 24. However, an honest comparison for generic $n$ would require considering the cosmic-brane back-reactions in holography, which is beyond our current computational control. It is possible that Fig. 24 fails to accurately capture the case for generic $n$ (for example, the transition to the decay phase in holography maybe continuous due to back-reactions. ). We leave such investigations for the future.

Now let us discuss the physics behind the three segments of the ”growth-plateau-decay” pattern.

During the linear growth phase, we observe that the growth rates (in units of the equilibrium entropy density) are the same among the three models. This is very similar the the linear growth behaviors of quenched state entanglement, where the rates are also equal among the models despite different entanglement dynamics. However, we should caution that the linear growth of the operator entanglement is not solely a consequence of the corresponding state entanglement increase for the subsystem $A$. During this phase, the operator state $|\rho_{AB}^{\alpha}(t)\rangle\in\mathcal{H}\otimes\mathcal{H}$ is building up entanglement between the subsystems $A\cup B$ and $A^{\prime}\cup B^{\prime}$, thus is drifting away from the form of a product state in $\mathcal{H}\otimes\mathcal{H}$. As a result, the operator entanglement growth rate is not twice of, but is equal to, the state rate. The linear growth stops at $t=\ell_{A}/2$. This coincides with the thermalization time for A, after which the operator entanglement barrier emerges. This suggests that the transition between the two phases is possibly related to the thermalization of region A.

The plateau phase represents an intriguing ”equilibrium” in the course towards thermalization. Despite the operator entanglement taking constant value, the operator state is certainly not static, for example the state Rényi entropy of AB is still linearly increasing. This excludes the possibility of describing the reduced density matrix by any kind of static (equilibrium) ensemble. We can ask a few questions regarding this phase. Firstly, is the calibrated plateau value universal? For example, could it represent the saturation to the maximal operator entanglement subject to the given constraints on the density matrix? Secondly, how universal is the plateau phase? For example, does it survive in chaotic systems away from the large $q$ or large $N$ limits we are considering? And finally, what are the implications for the observables? For example, does it implies any properties in the local correlation functions within region A?

The existence of the decay phase for operator entanglement is dictated by the quenched thermalization of these models. In the high temperature limit, the thermalized reduced density matrix is approximately proportional to identity. Viewed as an operator state, it is therefore short-range entangled with correlation length of order $\beta$. However, the specific behaviors are highly model dependent; in fact, the models we study have different mechanisms for triggering the decay phase.

The linear growth phase and the early plateau phase are the same among all three classes of models. It is the point of exit from the plateau and the subsequent decay phase that serve to differentiate among them. A similar story for the two-interval entanglement/mutual information in the quenched state is discussed in asplund_entanglement_2015 , where it is also found to encode the distinction between the chaotic and rational CFTs. However, if we consider the same quantity for the random unitary circuits, at large $q$ in large separation it exhibits the same behavior as the holographic CFTs 999On the other hand, they were shown to differ for the case of $n=2$ Rényi mutual inormation.. In contrast, by probing the later time behavior of the operator entanglement, the three models are fully distinguished.

Phenomenologically we observe that how late the system exits the plateau phase is correlated with how chaotic it is: the rational CFTs exit at $t=\frac{\ell_{B}}{2}$; the random unitary circuits exit at $t=\frac{\ell}{2}-\frac{\ell_{A}}{4}=\frac{\ell_{B}}{2}+\frac{\ell_{A}}{4}$; the holographic CFTs exit at $t=\frac{\ell}{2}=\frac{\ell_{B}}{2}+\frac{\ell_{A}}{2}$. By the underlying quenched dynamics, the reduced density matrix is fully thermalized at $t=\frac{\ell}{2}$, after which the operator state becomes short-range entangled and the operator entanglement has to drop to short-range value. Therefore the exit has to occur before $t=\frac{\ell}{2}$. This is saturated by a class of maximally chaotic systems: the holographic CFTs. It is therefore tempting to suggest that the persistence of the operator entanglement barrier, calibrated using the standard quenched state $e^{-\beta H/4}|B\rangle$, can serve as an alternative measure for the quantum chaoticity of the system.

To better check and understand the possible relation between the duration of the operator entanglement plateau and quantum chaoticity, it could be helpful in the future to study more examples of intermediate models, such as chaotic CFTs that are not holographic. Here based on the two chaotic models we study, it is still interesting to compare their mechanisms of exiting the plateau phase. In both models they are driven by geometrical transitions of the domain walls and the HRT surfaces respectively. The difference is as follows. In the random unitary circuits, the exit corresponds to a first order transition of the effective statistical mechanics, the operator entanglement is the free energy itself and thus is continuous across the exit. In particular, all the domain walls have the same tension in contributing to the free energy. In the holographic CFTs, the exit is driven by a phase transition of the bulk dual of the operator state, i.e. $\gamma^{AB}(t)$. The operator entanglement (or more specifically, reflective entropy) is given by the area of the HRT surface in the bulk, which changes by responding to the bulk dual transition, and thus is not continuous across the exit. In particular, contrary to the case in the random unitary circuits, $\gamma^{AB}(t)$ and $\gamma^{A}(t)$ play highly distinct roles in determining the operator entanglement. This distinction between the two models is not present when computing two-interval mutual information. This explains why the mutual information cannot distinguish between these two models as pointed out before. In the future, it would be interesting to connect this distinction with the different levels of quantum chaoticity between the two models.

Acknowledgements.

We acknowledge the accommodation and interactive environment of the KITP program “The Dynamics of Quantum Information”, where the question of “entanglement barrier” was raised up and stimulated this work. We acknowledge Xiao Chen, Thomas Faulkner, Yuya Kusuki, Andreas W.W. Ludwig, Márk Mezei, Adam Nahum, Frank Pollmann, Shinsei Ryu, and Kotaro Tamaoka for insightful discussion and previous collaborations on related projects. HW is supported by the Gordon and Betty Moore Foundation through Grant GBMF7392; TZ is supported by the Gordon and Betty Moore Foundation, under the EPiQS initiative, Grant GBMF4304, at the Kavli Institute for Theoretical Physics. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. We acknowledge support from the Center for Scientific Computing from the CNSI, MRL: an NSF MRSEC (DMR-1720256) and NSF CNS-1725797.

Appendix A AdS3/CFT2 Calculations

In this appendix we carry out the holographic calculations explicitly in AdS3/CFT2, where one can use the fact that BTZ black holes are locally isometric to AdS3 and can be obtained as a quotient of it. This would allow us to explicitly solve for the HRT surfaces/geodesics. To simplify the discussion, we will study the TFD state, dual to the (unwrapped) two-sided BTZ black hole. The answer for the pure B-state is simply given by half of the TFD results.

The right exterior of BTZ black hole geometry is:

[TABLE]

The AdS scale is set to 1, and the temperature is given by $\beta=2\pi/r_{0}$. The left exterior has the same geometry and can be obtained by analytically continue $r_{0}t\to r_{0}t+i\pi$. In addition, the future interior of the BTZ black hole can be obtained by $r_{0}t\to r_{0}t+i\pi/2,\;\rho\to i\alpha$ and has the geometry:

[TABLE]

Had we wrapped the $y$ direction to be compact, there will be a compact null circle at $\alpha=\pi/2$ and this corresponds to the black hole singularity. The two exterior and interior regions are all isometric to parts of the poincare patch for AdS3:

[TABLE]

via the identification of coordinates:

[TABLE]

and analytic continuations thereafter. Towards the asymptotic boundary of the right exterior $z\to 0$ ($\rho\to\infty$), the boundary coordinates can be directly identified:

[TABLE]

Similarly, towards the left exterior we have (by sending $t\to-t+i\pi$):

[TABLE]

Denote the end points of the interval $AB$ on the right asymptotic boundary by $(P_{1},P_{3})$, and the other end point of the sub-interval $A$ by $P_{2}$; the corresponding end points on the left asymptotic boundary are denoted by $P_{i=1,2,3}^{\prime}$. Taking the set up as in Sec. 2, they are mapped into the poincare-patch coordinates $(x_{0},x_{1})$ as:

[TABLE]

Via this map we can from now on carry out the analysis in the poincare patch, in which generic space-like geodesics $\left\{x_{0}(z),x_{1}(z)\right\}$ take the following form of “circles”:

[TABLE]

if the trajectories projected onto the $(x_{0},x_{1})$ plane are space-like; otherwise the geodesics take the form:

[TABLE]

A.1 Phase 1

In the first phase, the segments for the HRT surfaces $\gamma^{AB}(t)$ that construct the operator state $|\psi^{1/2}_{AB}(t)\rangle$, and the HRT surface $\gamma^{A}(t)$ that computes the operator entanglement, are half-circles $\gamma^{AB}(t)=\widehat{P_{1}P_{1}^{\prime}}\cup\widehat{P_{3}P_{3}^{\prime}},\;\gamma^{A}(t)=\widehat{P_{2}P_{2}^{\prime}}$ that connect $(P_{1},P_{1}^{\prime}),(P_{3},P_{3}^{\prime})$ and $(P_{2},P_{2}^{\prime})$ respectively, as studied in Ref. hartman:2013 ; shinsei:2018 . They all pass through the black hole interior, and are described by Eq. (71) with $\eta=0,C_{1}=0$ (see the left of Fig. 25):

[TABLE]

The operator entanglement for the TFD state is given by twice the geodesic length of $\widehat{P_{2}P_{2}^{\prime}}$ multiplied by $\frac{1}{4G_{N}}=\frac{c}{6\ell}$:

[TABLE]

where $\Lambda$ is the UV cut-off we placed towards the asymptotic boundary at $z=1/\Lambda$. The operator entanglement for the pure B-state is

[TABLE]

The same result also holds for the geodesic length of $\widehat{P_{1}P_{1}^{\prime}}$ and $\widehat{P_{3}P_{3}^{\prime}}$, after taking into account of a $y$-dependent UV $z$ cut-off, which corresponds to a uniform UV $\rho$ cut-off in the original geometry via $1/z\sim e^{\rho-r_{0}y}$. These are the linear growth identified in Ref. hartman:2013 .

A.2 Phase 2

In this phase, $\gamma^{AB}(t)$ remains the same phase as before, while $\gamma^{A}(t)$ undergoes a transition from a half-circle $\gamma^{A}(t)=\widehat{P_{2}P_{2}^{\prime}}$ to a union of two segments $\gamma^{A}(t)=\widehat{P_{2}\zeta}\cup\widehat{P_{2}^{\prime}\zeta^{\prime}}$ (see the center of Fig. 25). The two segments end orthogonally on $\widehat{P_{1}P_{1}^{\prime}}$ at $\zeta,\zeta^{\prime}\in\widehat{P_{1}P_{1}^{\prime}}$ respectively. Both segments are of the form Eq. (71) and one is the reflection of another across $x_{1}=0$, so we focus on the one of them $\widehat{P_{2}\zeta}$ and study the orthogonality condition between $\alpha=\widehat{P_{1}P_{1}^{\prime}},\beta=\widehat{P_{2}\zeta}$:

[TABLE]

The $(\pm)$ branches in $\widehat{P_{2}\zeta}$ is determined by whether the segment has passed the turning point at $z=R$. In order to satisfy orthogonality as well as the boundary condition we must solve:

[TABLE]

The last equation is satisfiable only if $\zeta$ is on the $(-)$ branch of $\widehat{P_{2}\zeta}$. This is a set of 5 equations for 5 variables $\zeta_{z},\eta,R,C_{1},C_{2}$, so there are a discrete number of solutions. One can eliminate the irrelevant integration constants $C_{1,2}$ and obtain the following equations:

[TABLE]

To obtain solutions, we define $\epsilon^{2}=e^{-r_{0}\ell_{A}}$ and $\kappa^{2}=e^{r_{0}\ell_{A}}\cosh^{-2}{(r_{0}t)}$. In the high temperature limit, $\epsilon\ll 1$ is a small parameter, it is easy to solve to the first few orders in $\epsilon$ expansion of the solution:

[TABLE]

In the current regime, $t>\ell_{A}/2$ and therefore $\kappa\ll 1$, one can further expand in $\kappa$ the leading order solution in $\epsilon$:

[TABLE]

We can use this result to compute the geodesic length of $\widehat{P_{2}\zeta}$:

[TABLE]

The operator entanglement for the TFD state is equal to 4 times $\frac{\mathcal{L}(P_{2},\zeta)}{4G_{N}}$, from both the left exterior region and the other copy across $\gamma^{AB}(t)$:

[TABLE]

The operator entanglement for the pure B-state is thus:

[TABLE]

A.3 Phase 3

In this phase, $\gamma^{AB}(t)$ undergoes a phase transition from $\gamma^{AB}(t)=\widehat{P_{1}P_{1}^{\prime}}\cup\widehat{P_{3}P_{3}^{\prime}}$ to $\gamma^{AB}(t)=\widehat{P_{1}P_{3}}\cup\widehat{P_{1}^{\prime}P_{3}^{\prime}}$. Responding to this background transition, $\gamma^{A}(t)$ transits from $\gamma^{A}(t)=\widehat{P_{2}\zeta}\cup\widehat{P_{2}^{\prime}\zeta^{\prime}}$ to $\gamma^{A}(t)=\widehat{P_{2}\chi}\cup\widehat{P_{2}^{\prime}\chi^{\prime}}$ that end perpendicularly on $\widehat{P_{1}P_{3}}$ and $\widehat{P_{1}^{\prime}P_{3}^{\prime}}$ at $\chi$ and $\chi^{\prime}$ respectively (see the right of Fig. 25).

We focus on the $\left(\widehat{P_{1}P_{3}},\widehat{P_{2}\chi}\right)$ system. The geodesics are collinear projected onto the $(x_{0},x_{1})$ plane, so the geodesic solutions are:

[TABLE]

The geodesic length of $\widehat{P_{1}P_{3}}$ is:

[TABLE]

Following a similar procedure as before, we can solve for the geodesic solution for $\widehat{P_{2}\chi}$ perturbatively in $e^{-r_{0}\ell_{B}}$ and $e^{-r_{0}\ell_{A}}$:

[TABLE]

The intersection point $\chi$ is on the $(-)$ branch, so the geodesic length of $\widehat{P_{2}\chi}$ is:

[TABLE]

This is exponentially suppressed. To relate these results with the CFT result Eq. (32) more precisely, we need to define the UV cut-off scale with a factor of $r_{0}=2\pi/\beta$. This can be understood as coming from the factor $r_{0}^{2}$ in the $(dt,dy)$ part of the metric in Eq. (65).

A.4 Extensive entanglement from behind the horizon

An interesting question for the second phase A.2 is whether and how much of the HRT surface $\gamma^{A}(t)\sim\widehat{P_{2}\zeta}$ remains in the black hole interior. To answer this, recall that the horizon for the BTZ black hole is mapped into the following null-sheets in the poincare patch of AdS3:

[TABLE]

emanating from the bifurcating surface at $\{x_{0}=x_{1}=0,z\in\mathbb{R}^{+}\}$. The black hole interior corresponds to $\{|x_{0}|>|x_{1}|,z\in\mathbb{R}^{+}\}$. From Eq. (76) and Eq. (80) it is easy to see that the end point $\zeta$ is in the interior, although very close to the horizon from behind:

[TABLE]

On the hand, it is also easy to work out where $\widehat{P_{2}\zeta}$ crosses the horizon by solving $x^{\beta}_{0}(z^{c})=x^{\beta}_{1}(z^{c})$ of Eq. (76):

[TABLE]

This is on the $(+)$ branch of $\widehat{P_{2}\zeta}$, i.e. before the turning point. The contribution to the operator entanglement from the black hole interior can thus be computed from:

[TABLE]

This captures the entire volume-law contribution to $S_{A}(t)$. In this phase, we have extensive operator entanglement from behind the horizon. To have a better picture for what $\gamma^{A}(t)$ looks like in black brane geometry, let us express these in terms of the $y$ coordinate representing the spatial dimension:

[TABLE]

Recall that the projection of $\gamma^{A}(t)$ onto the $y$-coordinate covers the interval $\{y\in\mathbb{R}:-\ell_{A}\leq y\leq 0\}$, the horizon-crossing occurs at:

[TABLE]

which stays fixed as we take $\ell_{A}\gg\beta$. This thus directly demonstrates that the interior segment $\{y\in\mathbb{R}:-\ell_{A}\leq y\leq y_{c}\}$ covers the entire volume contribution. Furthermore, the interior segment sticks very close to the horizon, though from behind.