Entanglement scaling of operators: a conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1  +  1d

In a Hilbert space with a bipartition $\mathcal{H} = \mathcal{H}_A \otimes \mathcal{H}_B$, any state $\left\vert{\psi}\right\rangle \in \mathcal{H}$ has a Schmidt decomposition

$$\begin{align*} \left\vert{\psi}\right\rangle = \sum_{i=1}^{r} \sqrt{\lambda_i} \left\vert{\psi_{A,i}}\right\rangle \otimes \left\vert{\psi_{B,i}}\right\rangle \nonumber \end{align*}$$

where the $\lambda_i$ are positive real coefficients, ordered as $\lambda_1 \geqslant \lambda_2 \geqslant \dots \geqslant \lambda_r>0$. r is the Schmidt rank, which cannot be larger than ${\rm min} ({\rm dim} \mathcal{H}_A, {\rm dim} \mathcal{H}_B)$. The $\left\vert{\psi_{A, i}}\right\rangle$'s (resp. $\left\vert{\psi_{B, i}}\right\rangle$) are an orthonormal set of states in $\mathcal{H}_A$ ($\mathcal{H}_B$). Although the Schmidt decomposition is not unique, the set of coefficients $\{\lambda_i \}$ is. For a normalized state, $\left\langle{\psi}\vert{\psi}\right\rangle = 1$, the entanglement entropy $S_\alpha (\left\vert{\psi}\right\rangle)$ is defined as

$$\begin{align} \label{eq:EE} S_\alpha (\left\vert{\psi}\right\rangle) = \frac{1}{1-\alpha} \log \left(\sum_{i=1}^r \lambda_i^\alpha \right), \qquad \alpha >0. \nonumber \end{align} \tag{ 1 }$$

(In this paper, we always refer to $S_\alpha$ as the entanglement entropy, without having $\alpha \rightarrow 1$ in mind; we do not give different names to the von Neumann and the Renyi entropies.) Because of the uniqueness of the $\lambda_i$'s, the entanglement entropy is basis independent: if U_A and U_B are two unitary transformations acting on $\mathcal{H}_A$ and $\mathcal{H}_B$ respectively, then $S_\alpha (U_A \otimes U_B \left\vert{\psi}\right\rangle)= S_\alpha (\left\vert{\psi}\right\rangle)$. For integer values of $\alpha \geqslant 2$, the entanglement entropy can be rewritten as the expectation value of an operator, dubbed swap operator by Zanardi, Zalka and Faoro [25], in a replicated system. One takes α identical replicas of the quantum state $\left\vert{\psi}\right\rangle$, thus creating a new state in a larger Hilbert space, $\left\vert{\psi}\right\rangle ^{\otimes \alpha} = \left\vert{\psi}\right\rangle \otimes \left\vert{\psi}\right\rangle \otimes \dots \otimes \left\vert{\psi}\right\rangle \in \mathcal{H}^{\otimes \alpha}$. The swap operator is the operator that cyclically permutes the α replicas of the subsystem A, leaving B untouched; we write it as $\mathcal{S}^A_{(1, 2, 3, \dots, \alpha)}$. The subscript $(1, 2, 3, \dots, \alpha)$ indicates the permutation of α elements, written in cycle notation; notice that here the permutation σ has a single cycle $(1, \sigma (1), \sigma\circ \sigma (1), {\rm etc})$ of total length α. The entanglement entropy is then

$$\begin{align} \label{eq:EE_replicas} (\alpha \; {\rm integer} \; \geqslant 2) \qquad \quad S_\alpha (\left\vert{\psi}\right\rangle) = \frac{1}{1-\alpha} \log \left(\left\langle {\psi} \right\vert ^{\otimes \alpha} \; \mathcal{S}^A_{(1,2,3,\dots,\alpha)} \; \left\vert{\psi}\right\rangle ^{\otimes \alpha} \right). \nonumber \end{align} \tag{ 2 }$$

The main advantage of the latter expression over equation (1) is that it makes it clear that the quantity inside the logarithm is an observable [25]; it allows to calculate the entanglement entropy in Monte-Carlo simulations [26], it is the basic idea of various proposals to measure the entanglement entropy [27–29], and it is at heart of the recent groundbreaking experiment [30] that managed to measure S₂ in a cold atom setup. It is also the key to field-theory calculations based on the replica trick [15], although equation (2) usually appears in a somewhat disguised form in that context, as those field theory calculations are typically formulated as path integrals on a replicated worldsheet that involve twist operators (we will come back to this point in section 2.2).

Now we turn to the case where $\mathcal{H} = (\mathbb{C}^d){\hspace{0pt}}^{\otimes L}$ is the Hilbert space of a spin chain; d is the dimension of the local degree of freedom on each site, and L is the total number of sites. The subsystem A (resp. B) is made of the first L_A sites ($L_B= L-L_A$ sites). Then one says that $\left\vert{\psi}\right\rangle$ obeys the area law (for a given α) if $S_\alpha (\left\vert{\psi}\right\rangle)$ remains bounded when $L_A, L_B \rightarrow \infty$. A famous theorem of Hastings [12] states that when $\left\vert{\psi}\right\rangle$ is the ground state of a gapped local Hamiltonian, it obeys the area law, with $\alpha =1$. Strictly speaking, the case $\alpha=1$ is usually not sufficient to ensure that a state can be approximated by an MPS, see [10]—although, for ground states of gapped Hamiltonians, Hastings [12] did prove the approximability as well—. As shown by Cirac and Verstraete [9], one sufficient condition for approximability, that is valid beyond ground states, is to obey the area law for $0<\alpha <1$. Then $\left\vert \psi \right\rangle$ can be approximated as

$$\begin{align*} \left\vert \psi \right\rangle \simeq \sum_{\{\sigma_1,\dots,\sigma_L \} } \sum_{\{a_1,\dots,a_{L-1} \}} [M^{\sigma_1}]_{1 a_1} [M^{\sigma_2}]_{a_1 a_2} [M^{\sigma_3}]_{a_2 a_3} \dots [M^{\sigma_L}]_{a_{L-1} 1} \left\vert \sigma_1, \dots,\sigma_L \right\rangle \nonumber \end{align*}$$

with finite matrices $M^{\sigma_i}$. The latter area law, for $0<\alpha<1$, can also be proved for the ground state of gapped Hamiltonians, see for instance the discussion in [31].

Notwithstanding these important aspects concerning the role of the (Renyi) parameter α, it is fair to say that in most known situations, the behavior of $S_\alpha (\left\vert{\psi}\right\rangle)$ does not appear to depend dramatically on the precise value of α, and it is often the case in one-dimensional systems that the area law is either obeyed for all $\alpha >0$, or violated for all $\alpha >0$. Another way to say this is that, usually, there is no phase transition as one varies α; notice that this is of course a requirement if one wants to use the replica trick. The analyticity in α is known to hold, in particular, for the (single interval) ground state entanglement entropy at quantum critical points described by CFT, where the area law is always violated logarithmically, no matter the precise value of α. So in this paper, we will proceed as in the seminal work of Calabrese and Cardy [15], and many subsequent papers: we will rely on the replica trick, calculate the OSEE for α integer $\geqslant 2$, and then make the assumption that the scaling one finds (as a function of system size, or of time) is similar for other values of α.

Instead of approximating states by MPSs, we now want to approximate operators by MPOs. Mathematically, this is essentially the same problem though, since operators on $\mathcal{H}$ are nothing but elements of $\mathcal{H}^{\prime} = {\rm End}(\mathcal{H})$, which is again a Hilbert space with the Hilbert–Schmidt inner product $\left\langle{O_1}\vert{O_2}\right\rangle ^{\prime} = {\rm tr [O_1^\dagger O_2]}$, $O_1, O_2 \in {\rm End}(\mathcal{H})$. (To avoid mathematical difficulties, in this introduction we assume that we are dealing with a finite dimensional $\mathcal{H}$, so it is always true that ${\rm End}(\mathcal{H})$ is also a Hilbert space.) So one can simply replace $\mathcal{H}_A$ by $\mathcal{H}_A^{\prime} = {\rm End}(\mathcal{H}_A)$ and $\mathcal{H}_B$ by $\mathcal{H}_B^{\prime} = {\rm End}(\mathcal{H}_B)$ in the previous section; this straightforwardly leads to the notion of a Schmidt decomposition for the operator O, which we chose to write as

$$\begin{align} \label{eq:schmidt_op} \frac{O}{\sqrt{{\rm tr } [O^\dagger O]} } = \sum_{i=1}^{r} \sqrt{\lambda_i} O_{A,i} \otimes O_{B,i}, \nonumber \end{align} \tag{ 3 }$$

where the $\lambda_i$ are positive real coefficients, in decreasing order as above. The operators $O_{A, i} \in {\rm End}(\mathcal{H}_A)$ obey the orthonormality condition ${\rm tr} [O^\dagger_{A, i} O_{A, j}] = \delta_{i, j}$ (same for the O_B,i's). In equation (3) we have explicitly included a normalization factor in the left-hand side, because it should be emphasized that, while it is perfectly standard to work with quantum states that are properly normalized ($\left\langle{\psi}\vert{\psi}\right\rangle =1$), it is unusual to ask that operators on $\mathcal{H}$ should obey a specific normalization condition. Again, the set of coefficients $\{\lambda_i \}$ is unique, and is invariant under changes of basis of $\mathcal{H}_A$ and $\mathcal{H}_B$. Thus, we are naturally led to the notion of operator space entanglement entropy (OSEE), which, to the best of our knowledge, was investigated first by Zanardi in the context of [1], then rediscovered by Prosen and Pižorn [2], with subsequent studies in [32–34],

$$\begin{align} \label{eq:Renyi_alpha_op} S_\alpha (O) = \frac{1}{1-\alpha} \log \left(\sum_{i=1}^r \lambda_i^\alpha \right). \nonumber \end{align} \tag{ 4 }$$

As in the previous section, when α is an integer $\geqslant 2$, one can rewrite $S_\alpha(O)$ in a way that involves the swap operator in a replicated system. The expression involves α replicas of O, $O^{\otimes \alpha} = O \otimes \dots \otimes O$, their complex conjugates, and the swap operator that cyclically permutes the replicas of subsystem A, $\mathcal{S}^A_{(1, 2, 3, \dots, \alpha)}$, as well as its inverse $\mathcal{S}^A_{(1, \alpha, \alpha-1, \dots, 2)}$—recall that we use the cycle notation for the permutations $(1, \sigma(1), \sigma \circ \sigma(1), {\rm etc})$—. The expression reads

$$\begin{align} \label{eq:OSEE_replicas} (\alpha \; {\rm integer}\; \geqslant 2) \qquad \quad S_\alpha (O) = \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} [ (O^\dagger)^{\otimes \alpha} \cdot \mathcal{S}^A_{(1,\alpha,\alpha-1,\dots,2)} \cdot O^{\otimes \alpha} \cdot \mathcal{S}^A_{(1,2,3,\dots,\alpha)} ]}{{\rm tr} \left[ (O^\dagger)^{\otimes \alpha}\cdot O^{\otimes \alpha} \right] } \right). \nonumber \end{align} \tag{ 5 }$$

This is the key expression that we will use throughout this paper to calculate the OSEE; notice that the denominator inside the logarithm could also be written as $\left({\rm tr} [O^\dagger O] \right){\hspace{0pt}}^\alpha$. Notice also that, if one is willing to relate the OSEE of, say, a density matrix ρ, to a physical observable, then equation (5) does not quite do the job, because the traces involve expressions that are quadratic in the replicated density matrix $\rho^{\otimes \alpha}$, while expectation values must be linear in $\rho^{\otimes \alpha}$. It is not difficult to cure this though; but we defer this discussion to the next section.

For now, let us focus on the relevance of the OSEE to the approximability by MPOs. Since the OSEE is exactly the same quantity as the 'usual' entanglement entropy once one has replaced the space of states $\mathcal{H}$ by the space of operators $\mathcal{H}^{\prime}={\rm End}(\mathcal{H})$, one can make exactly the same statements as above. Namely, if we consider again the spin chain $\mathcal{H} = (\mathbb{C}^d){\hspace{0pt}}^{\otimes L}$ with subsystems A and B corresponding to the first L_A sites and the last $L_B= L-L_A$ sites, then one says that O obeys the operator area law (for a given α) if $S_\alpha (O)$ remains bounded when $L_A, L_B, \rightarrow \infty$. Then the result of Cirac and Verstraete [9] shows that if O obeys the area law for $0<\alpha <1$, then it can be efficiently approximated (in Hilbert–Schmidt norm) by an MPO,

$$\begin{align*} O \simeq \sum_{{\{\sigma_1,\dots,\sigma_L \}}\atop{\{\sigma^{\prime}_1,\dots,\sigma^{\prime}_L \}} } \sum_{\{a_1,\dots,a_{L-1} \}} [M^{\sigma_1^{\prime}}_{\sigma_1}]_{1 a_1} [M^{\sigma_2^{\prime}}_{\sigma_2}]_{a_1 a_2} [M^{\sigma_3^{\prime}}_{\sigma_3}]_{a_2 a_3} \dots [M^{\sigma_L^{\prime}}_{\sigma_L}]_{a_{L-1} 1} \left\vert \sigma^{\prime}_1, \dots,\sigma^{\prime}_L \right\rangle \left\langle \sigma_1, \dots,\sigma_L \right\vert, \nonumber \end{align*}$$

with finite matrices $M^{\sigma_i^{\prime}}_{\sigma_i}$. The discussion of [10] applies to other values of α.

Thus, in trying to approximate operators by MPOs, it is important to understand whether or not they violate the operator area law, and if they do violate it, whether this violation is logarithmic, or worse. Prosen and Pižorn and collaborators investigated this question in a few concrete examples [2, 32–34]; in this paper, we revisit and extend their results, and in particular, we set up the framework for calculating the OSEE in CFT.

1.3.1. The OSEE of $\boldsymbol{ \left\vert{\psi}\right\rangle \left\langle {\psi} \right\vert }$ is exactly twice the entanglement entropy of $\boldsymbol{\left\vert{\psi}\right\rangle }$.

Let us emphasize a simple but important point: the usual entanglement entropy is defined for states, while the OSEE is defined for operators. They do not probe the same class of objects, and it is therefore vain to look for a general relation between the two quantities. That being said, there is one case where they are obviously related: when the operator O is the density matrix of a pure state, $\left\vert {\psi} \right\rangle \left\langle {\psi} \right\vert$. In that case, the OSEE is exactly twice the entanglement entropy of $\left\vert {\psi} \right\rangle$,

$$\begin{align} \label{eq:OSEE_pure} S_\alpha(\left\vert {\psi} \right\rangle \left\langle {\psi} \right\vert) = 2 S_\alpha(\left\vert {\psi} \right\rangle). \nonumber \end{align} \tag{ 6 }$$

This simply follows from the fact that a Schmidt decomposition of the state $\left\vert {\psi} \right\rangle$ automatically induces a Schmidt decomposition for the operator $\left\vert {\psi} \right\rangle \left\langle {\psi} \right\vert$. This simple relation will show up again in several places in this paper.

1.3.2. The OSEE of a density matrix ρ is measurable.

For integer $\alpha \geqslant 2$ the OSEE of a density matrix ρ is measurable. Indeed, using the hermiticity of ρ, elementary manipulations allow to rewrite equation (5) in terms of expectation values in a system with $2\alpha$ replicas—instead of α replicas in equation (5)—,

$$\begin{align*} S_\alpha (\rho) = \frac{1}{1-\alpha}\log \left(\frac{{\rm tr} \left[ \rho^{\otimes 2\alpha} \cdot \mathcal{S}^{A\cup B}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)} \cdot \mathcal{S}^{A}_{(1,2,3,\dots,\alpha)(\alpha+1,2\alpha,2\alpha-1,\dots,\alpha+2)} \right]}{{\rm tr} \left[ \rho^{\otimes 2 \alpha} \cdot \mathcal{S}^{A\cup B}_{(1,\alpha+1)(2,\alpha+2)\dots (\alpha,2\alpha)} \right] } \right). \nonumber \end{align*}$$

Notice that the different swap operators involved act on different subsystems: $\mathcal{S}^A$ acts on the replicas of subsystem A only, while $\mathcal{S}^{A \cup B}$ acts on the full system; notice also that the permutations are different. This formula shows that $S_\alpha (\rho)$ is measurable, at least in principle, since both the numerator and the denominator in the logarithm are expectation values of observables. Notice also the identity

$$\begin{align*} & \mathcal{S}^{A \cup B}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)} \mathcal{S}^{A}_{(1,2,3,\dots,\alpha)(\alpha+1,2\alpha,2\alpha-1,\dots,\alpha+2)} \nonumber \\ & \qquad = \mathcal{S}^{A}_{(1,\alpha+2)(2,\alpha+3)\dots (\alpha-1,2\alpha) (\alpha,\alpha+1)} \mathcal{S}^{B}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)}, \nonumber \end{align*}$$

which follows from the composition of the two permutations in part A, and which leads to a more symmetric expression for the OSEE,

$$\begin{align} \label{eq:OSEE_replicas_rho} S_\alpha (\rho) = \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} \left[ \rho^{\otimes 2\alpha} \cdot \mathcal{S}^{A}_{(1,\alpha+2)(2,\alpha+3)\dots (\alpha-1,2\alpha) (\alpha,\alpha+1)} \cdot \mathcal{S}^{B}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)} \right]}{{\rm tr} \left[ \rho^{\otimes 2 \alpha} \cdot \mathcal{S}^{A\cup B}_{(1,\alpha+1)(2,\alpha+2)\dots (\alpha,2\alpha)} \right] } \right). \nonumber \end{align} \tag{ 7 }$$

1.3.3. The OSEE is not a measure of entanglement in a mixed state.

In the next section, we illustrate how the OSEE can be calculated in practical situations, and we introduce the CFT tools that are needed later; in particular, we discuss the case of thermal density matrices. In section 3, we discuss the OSEE of the reduced density matrix after a global quench, the approach to the thermal density matrix or generalized Gibbs ensemble, and the consequences for numerical simulations of the evolution of one-dimensional quantum systems at large times. In section 4, we turn to the evolution operator ${\rm e}^{- {\rm i} H t}$, and argue that its OSEE typically grows linearly in time, which is reminiscent of the ballistic spreading of entanglement. In section 5, we report an attempt at understanding the OSEE of local operators in Heisenberg picture, which is relevant to Heisenberg-picture DMRG. We conclude in section 6. More details about free fermion calculations, permutations and twist operators in replicated CFT, and about the reported attempt of section 5, are included in three appendices.

Consider the operator T that translates all sites by one (modulo L) to the right: it acts on the basis states of $\mathcal{H}=(\mathbb{C}^d){\hspace{0pt}}^{\otimes L}$ as $T \left\vert{i_1, i_2, \dots, i_L}\right\rangle = \left\vert{i_L, i_1, \dots, i_{L-1}}\right\rangle$. Equivalently, the components of T are $[T]_{i_1 i_2 \dots i_L}^{i_1^{\prime} i_2^{\prime} \dots i_L^{\prime}} = \delta_{i_L, i_1^{\prime}} \delta_{i_1, i_2^{\prime}} \dots \delta_{i_{L-1}, i_{L}^{\prime}}$. It is easy to see that this operator is an MPO with bond dimension d: for each pair $(i, i^{\prime})$ define the $d \times d$ matrix $M^{i^{\prime}}_{i}$ that has elements $[M^{i^{\prime}}_{i}]_{a b} = \delta_{a, i} \delta_{b, i^{\prime}}$; then

$$\begin{align*} [T]_{i_1 i_2 \dots i_L}^{i_1^{\prime} i_2^{\prime} \dots i_L^{\prime}} = {\rm tr} \left[ M^{i_1^{\prime}}_{i_1} M^{i_2^{\prime}}_{i_2} \dots M^{i_L^{\prime}}_{i_L} \right]. \nonumber \end{align*}$$

Now take $\mathcal{H}_A$ (resp. $\mathcal{H}_B$) as the first L_A sites (resp. $L_B=L-L_A$ last sites) of the chain, and let us calculate the OSEE. We need to introduce the 'partial translators' T_A,a,b, defined as

$$\begin{align*} [T_{A,a,b}]_{i_1 i_2 \dots i_{L_A}}^{i_1^{\prime} i_2^{\prime} \dots i_{L_A}^{\prime}} = \delta_{a, i_1^{\prime}} \delta_{i_1, i_2^{\prime}} \dots \delta_{i_{L_A},b}, \nonumber \end{align*}$$

as well as partial translators on part B, T_B,b,a, with a similar definition. We have

$$\begin{align*} T = \sum_{a=1}^d \sum_{b=1}^d T_{A,a,b} \otimes T_{B,b,a}. \nonumber \end{align*}$$

This is close to the Schmidt decomposition (3), however we still need to check that the different term are orthonormal. First, notice that in the left-hand side, T has norm ${\rm tr} [T^\dagger T] = d^L$, which follows from the fact that $T^\dagger T$ is the identity on $\mathcal{H}$. Second, the inner product between two partial translators in the right-hand side is

$$\begin{align*} {\rm tr} \left[ T_{A,a,b}^\dagger T_{A,a^{\prime},b^{\prime}} \right] &= \sum_{\{i\},\{i^{\prime}\} } [T^\dagger_{A,a,b}]_{i_1^{\prime} i_2^{\prime} \dots i_{L_A}^{\prime}}^{i_1 i_2 \dots i_{L_A}} \times [T_{A,a^{\prime},b^{\prime}}]_{i_1 i_2 \dots i_{L_A}}^{i_1^{\prime} i_2^{\prime} \dots i_{L_A}^{\prime}} \nonumber \\ &= \sum_{\{i\},\{i^{\prime}\} } \delta_{a,i_1^{\prime}} \delta_{i_1, i_2^{\prime}} \dots \delta_{i_{L_A},b} \times \delta_{a^{\prime},i_1^{\prime}} \delta_{i_1,i_2^{\prime}} \dots \delta_{i_{L_A},b^{\prime}} \nonumber \\ &= d^{L_A-1} \delta_{a,a^{\prime}} \delta_{b,b^{\prime}}, \nonumber \end{align*}$$

so the partial translators are orthogonal, with some fixed normalization factor. Thus, the correct way to write the Schmidt decomposition of T, with properly normalized terms, is

$$\begin{align*} \frac{1}{d^{L/2}} T = \sum_{a=1}^d \sum_{b=1}^d \frac{1}{d} \left(\frac{1}{d^{(L_A-1)/2}} T_{A,a,b} \right) \otimes \left(\frac{1}{d^{(L_B-1)/2}} T_{B,b,a} \right). \nonumber \end{align*}$$

There are exactly d² non-zero Schmidt coefficients, all equal to 1/d. It follows that

$$\begin{align} S_\alpha(T) = \frac{1}{1-\alpha} \log \left(\sum_{a,b} (1/d^2)^\alpha \right) = 2 \log d, \nonumber \end{align} \tag{ 8 }$$

independently of α. Thus, in this simple example, the OSEE is simply the log of the bond dimension of the MPO, counted twice because there are two boundary points that separate A and B.

Alternatlvely, instead of writing the Schmidt decomposition explicitly, we could have used the formula with the replicas (5), leading to the following amusing calculation. Let us represent the operator T graphically, as

where each line corresponds to a Kronecker delta $\delta_{i, i^{\prime}}$. Then, graphically, one just has to count the number of closed loops; each loop gives a factor d. For instance, the squared norm of T is given by

Now let us evaluate ${\rm tr} [ (T^\dagger){\hspace{0pt}}^{\otimes \alpha} \cdot \mathcal{S}_{(1, \alpha, \alpha-1, \dots, 2)} \cdot T^{\otimes \alpha} \cdot \mathcal{S}_{(1, 2, 3, \dots, \alpha)} ]$. For simplicity, we draw the case $\alpha=2$ only. We represent the replicated T as

and the swap operator that permutes the two replicas of part A as

Then we find, counting the number of closed loops on the diagram, that

More generally, for $\alpha \geqslant 2$, one finds ${\rm tr} [ (T^\dagger){\hspace{0pt}}^{\otimes \alpha} \cdot \mathcal{S}_{(1, \alpha, \alpha-1, \dots, 2)} \cdot T^{\otimes \alpha} \cdot \mathcal{S}_{(1, 2, 3, \dots, \alpha)} ]= d^{\alpha (L-1)}$. Plugging this into equation (5), we recover the fact that $S_\alpha(T) = 2 \log d$, independently of α.

Next, consider an infinitely long critical system described by conformal field theory, with a Hamiltonian H and gapless excitations that propagate at velocity v. The system is at inverse temperature β. If we think of this setup as being the effective, large-scale, description of some microscopic model, then we have to assume that $\beta v \gg a_0$ where a₀ is some characteristic UV length scale (e.g. lattice spacing or Fermi wavelength).

The thermal density matrix is

$$\begin{align} \rho_\beta = \frac{1}{Z} {\rm e}^{-\beta H}, \nonumber \end{align} \tag{ 9 }$$

where the normalization factor $Z = {\rm tr} [{\rm e}^{-\beta H}]$ is the partition function of the CFT on an infinitely long cylinder of circumference $\beta v$, see figure 3(a). Strictly speaking, this partition function is IR divergent; it is convenient to introduce an IR cutoff by cutting the cylinder at its two extremities, leaving a finite cylinder of length $\Lambda \gg \beta v$. It is then a standard result that the partition function goes as

$$\begin{align} \label{eq:Zcylinder} Z \simeq {\rm e}^{\frac{\pi c}{6} \frac{\Lambda}{\beta v} }, \nonumber \end{align} \tag{ 10 }$$

where c is the central charge. Now, let us briefly review the calculation of the 'usual' entanglement entropy. This will provide an opportunity to emphasize the well-known but crucial fact that the 'usual' entanglement entropy obeys a volume law, as opposed to the area law obeyed by the OSEE; it will also be a convenient way to introduce twist fields, which we will need later on.

2.2.1. Preliminary: the volume law for the 'usual' entanglement entropy.

Take A  =  [0, L_A] and $B = (-\infty, 0) \cup (L_A, +\infty)$, and consider $S^{({\rm EE})}_\alpha = \frac{1}{1-\alpha} \log \left({\rm tr} [\rho_A^\alpha] \right)$, for α integer. Here we use the superscript (EE) to emphasize that this quantity is not the OSEE $S_\alpha(\rho_\beta)$. As briefly reviewed in the introduction, $S_\alpha^{({\rm EE})}$ may also been written as the expectation value of the swap operator,

$$\begin{align*} S^{({\rm EE})}_\alpha = \frac{1}{1-\alpha} \log \left({\rm tr} [ \mathcal{S}^A_{(1,2,\dots,\alpha)} \cdot (\rho_\beta)^{\otimes \alpha} ] \right) = \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} [ \mathcal{S}^A_{(1,2,\dots,\alpha)} \cdot ({\rm e}^{- \beta H})^{\otimes \alpha} ] }{Z^\alpha } \right). \nonumber \end{align*}$$

The key point from [15] is that $Z^{\rm twist}_\alpha \equiv {\rm tr} [ \mathcal{S}^A_{(1, 2, \dots, \alpha)} \cdot ({\rm e}^{- \beta H}){\hspace{0pt}}^{\otimes \alpha}]$ is again the partition function of the CFT, but this time on a more complicated surface, see figure 3(b): it is a cyclinder with α identical sheets, where one can go continuously from one sheet to the next by turning around the two points (x,y)  =  (0,0) and (0,L_A)—namely, the two end points of the interval A—. It is possible to calculate this partition function $Z^{\rm twist}_\alpha$ directly, but it is more efficient to view the ratio of partition functions $Z^{\rm twist}_\alpha / Z^\alpha$ as a two-point correlator of twist fields $\mathcal{T}_{(1, \alpha, \alpha-1, \dots, 2)}$ and $\mathcal{T}_{(1, 2, \dots, \alpha)}$ located at (0,0) and (0,L_A) respectively (by convention, the subscript σ in $\mathcal{T}_\sigma$ is the permutation of the replicas one picks when one turns clockwise around the point; see also the appendix for a discussion of twist fields). The twist fields are primary operators with scaling dimension $\Delta_\alpha = \frac{c}{12} (\alpha - 1/\alpha)$. The two-point function of primary operators on the cylinder is of course entirely fixed by conformal invariance,

$$\begin{align*} \frac{Z^{\rm twist}_\alpha}{Z^\alpha} = \left\langle \mathcal{T}_{(1,\alpha,\dots,2)}(0,0) \mathcal{T}_{(1,2,\dots,\alpha)}(L_A,0) \right\rangle = \frac{1}{[ \frac{\beta v}{\pi } \sinh (\frac{\pi L_A}{\beta v}) ]^{2 \Delta_\alpha}}, \nonumber \end{align*}$$

which finally leads to

$$\begin{align} \label{eq:thermal_EE} S^{({\rm EE})}_\alpha = \frac{c}{6} \left(1 + \frac{1}{\alpha} \right) \log \left(\frac{\beta v}{\pi } \sinh (\frac{\pi L_A}{\beta v}) \right) \; \mathop{\simeq}\limits_{L_A \gg \beta v } \; \left(1 + \frac{1}{\alpha} \right) \frac{\pi c}{6 \beta v} L_A. \nonumber \end{align} \tag{ 11 }$$

It is extensive: the entanglement entropy satisfies a volume law (as illustrated in figure 2). However, most of this entropy is just thermal entropy. Indeed, if we calculate the (Renyi version of the) classical thermodynamic entropy at temperature β, we find:

$$\begin{align*}S_\alpha^{\rm (class.)} =&~ \frac{1}{1-\alpha} \log \left({\rm tr} [ (\rho_\beta)^\alpha] \right) \nonumber \\ \nonumber=&~ \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} [ {\rm e}^{- \alpha \beta H} ]}{({\rm tr} [ {\rm e}^{- \beta H} ])^\alpha} \right) \nonumber \\ & \mathop{\simeq}\limits_{\Lambda \gg \beta v} \frac{1}{1-\alpha} \log \left(\frac{{\rm e}^{\frac{\pi c}{6} \frac{\Lambda}{\alpha \beta v} } }{{\rm e}^{\frac{\pi c}{6} \frac{\Lambda}{\beta v} \alpha } } \right) = \left(1+\frac{1}{\alpha} \right) \frac{\pi c }{6 \beta v} \Lambda, \nonumber \end{align*}$$

where we used equation (10) in the last line. Thus, the entropy per unit length is $\left(1+\frac{1}{\alpha} \right) \frac{\pi c }{6 \beta v}$, which precisely coincides with the prefactor of equation (11), demonstrating that the extensive part of the 'usual' entanglement entropy at finite temperature entirely comes from thermal entropy, and has nothing to do with quantum entanglement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.2.2. OSEE from CFT, and an operator area law for $\boldsymbol{\rho_\beta}$.

Now, instead of the 'usual' entanglement entropy, consider the OSEE. To calculate the latter in CFT, we use again the replica trick, in the form

$$\begin{align*} S_\alpha (\rho_\beta) = \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} [\rho_\beta^{\otimes \alpha} \cdot \mathcal{S}^A_{(1,\alpha,\alpha-1,\dots,2)} \cdot \rho_\beta^{\otimes \alpha} \cdot \mathcal{S}^A_{(1,2,\dots,\alpha)}] }{({\rm tr} [ (\rho_\beta)^2 ])^{\alpha} } \right). \nonumber \end{align*}$$

Exactly like in the calculation of the entanglement entropy which we just reviewed, the two traces inside the logarithm can be regarded as CFT partition functions on certain multi-sheeted surfaces. Again, it is useful to view the ratio of these two partition functions as a correlation function of twist operators, living on an infinitely long cylinder of circumference $2 \beta v$, as illustrated in figure 3(c). If we parametrize the cylinder by the complex coordinate $x+{\rm i} y$ with $(x, y) \in \mathbb{R} \times [0, 2\beta v]$, then the four twist operators are located at the points 0, L_A, $i \beta v$ and $L_A + {\rm i} \beta v$, and the four-point function we need to calculate is (see figure 3)

$$\begin{align*} \frac{{\rm tr} [\rho^{\otimes \alpha} \cdot \mathcal{S}^A_{\sigma^{-1}} \cdot \rho^{\otimes \alpha} \cdot \mathcal{S}_\sigma ] }{({\rm tr} [ \rho^2 ])^{\alpha} } = \left\langle \mathcal{T}_{\sigma^{-1}} (0) \mathcal{T}_\sigma ({\rm i}\beta v) \mathcal{T}_{\sigma} (L_A) \mathcal{T}_{\sigma^{-1}} ({\rm i}\beta v+L_A) \right\rangle, \nonumber \end{align*}$$

with the cyclic permutation $\sigma = (1, 2, \dots, \alpha)$. Since this is a four-point function, it is not entirely fixed by conformal invariance. In fact, the four-point function of twist operators is an object that has attracted a lot of attention in recent years in the context of the two-interval entanglement entropy[36–38], and it is now well-established that is depends on the full spectrum of the theory one is working with. Therefore, the calculation of the four-point function is usually complicated, and we want to avoid such complications here. So we stick to the information that can be most easily extracted, without relying on any specific details of the CFT. This concerns the two regimes $L_A \ll \beta v$ and $L_A \gg \beta v$, where the four-point function factorizes into a product of two two-point functions; the latter are easily calculated. This gives

$$\begin{align*} L_A\ll \beta v: &\,\,\, \left\langle \mathcal{T}_{\sigma^{-1}} (0) \mathcal{T}_\sigma ({\rm i}\beta v) \mathcal{T}_{\sigma} (L_A) \mathcal{T}_{\sigma^{-1}} ({\rm i}\beta v+L_A) \right\rangle \nonumber \\ & \quad \simeq \; \left\langle \mathcal{T}_{\sigma^{-1}} (0)\mathcal{T}_\sigma (L_A) \right\rangle \left\langle \mathcal{T}_{\sigma^{-1}} ({\rm i}\beta v+L_A) \mathcal{T}_\sigma ({\rm i}\beta v) \right\rangle = \frac{1}{[ \frac{2\beta v}{\pi } \sinh (\frac{\pi L_A}{2\beta v}) ]^{2 \Delta_\alpha}} \frac{1}{[ \frac{2\beta v}{\pi } \sinh (\frac{\pi L_A}{2\beta v}) ]^{2 \Delta_\alpha}} \nonumber \\ L_A \gg \beta v: &\,\,\, \left\langle \mathcal{T}_{\sigma^{-1}} (0) \mathcal{T}_\sigma ({\rm i}\beta v) \mathcal{T}_{\sigma} (L_A) \mathcal{T}_{\sigma^{-1}} ({\rm i}\beta v+L_A) \right\rangle \nonumber \\ & \quad \simeq \; \left\langle \mathcal{T}_{\sigma^{-1}} (0)\mathcal{T}_\sigma (i \beta v) \right\rangle \left\langle \mathcal{T}_{\sigma^{-1}} ({\rm i}\beta v+L_A) \mathcal{T}_\sigma (L_A) \right\rangle = \frac{1}{[ \frac{2\beta v}{\pi } \sin (\frac{\pi \beta v}{2\beta v}) ]^{2 \Delta_\alpha}} \frac{1}{[ \frac{2\beta v}{\pi } \sin (\frac{\pi \beta v}{2\beta v}) ]^{2 \Delta_\alpha}}, \nonumber \end{align*}$$

which leads to the following asymptotic behavior for the OSEE,

$$\begin{align} \label{eq:CFT_twice_EE} L_A \ll \beta v : \qquad S_\alpha(\rho_\beta) \simeq \frac{c}{3}\left(1+ \frac{1}{\alpha}\right) \log \left(\frac{L_A}{a_0} \right) \nonumber \end{align} \tag{ 12a }$$

$$\begin{align} \label{eq:thermal_area_law} L_A \gg \beta v : \qquad S_\alpha(\rho_\beta) \simeq \frac{c}{3}\left(1+ \frac{1}{\alpha}\right) \log \left(\frac{\beta v}{a_0 } \right), \nonumber \end{align} \tag{ 12b }$$

where we have reinstated the UV cutoff a₀ in order to get dimensionless expressions. Of course, strictly speaking, the results hold only for α integer, however we see now that the final expressions are analytic in α, and so it is very plausible that they hold for any real positive α. Such reasoning is standard in the case of the 'usual' entanglement entropy; we are simply repeating these well-known arguments here in the case of the OSEE; moreover, the analyticity in α will be supported by numerical evidence below.

Zoom In Zoom Out Reset image size

The physical meaning of equation (12a) is clear: at low temperatures, the reduced density matrix $\rho_A$ is very close to the one of the ground state $\left\vert{\psi_0}\right\rangle$, whose entanglement entropy is $S_\alpha(\left\vert{\psi_0}\right\rangle) = \frac{c}{6} (1+1/\alpha) \log (L_A/a_0)$. Thus, equation (12a) is just a particular occurrence of the relation (6). Equation (12b) is more exciting: it is the operator area law for the thermal density matrix $\rho_\beta$ (see figure 2). It shows that the OSEE is bounded at finite temperature; moreover, it increases logarithmically with inverse temperature β. As discussed in the introduction, this has an interesting practical consequence, which is that $\rho_\beta$ may be efficiently represented by an MPO with finite bond dimension. In addition, we see from equation (12b) that, when one approaches zero temperature, the bond dimension grows polynomially with inverse temperature β. We thus recover, within the framework of CFT, some of the conclusions obtained in [32]. Let us also emphasize that the question of approximability of thermal density matrices by MPOs or their higher dimensional versions has been tackled directly, with rigorous statements and proofs, in [39, 40].

2.2.3. The free fermion case.

Contrary to most CFTs, for the free Dirac fermion (c  =  1), the four-point function of twist operators is known (more generally, all correlations of twist operators are known in that case, and they are given in [15]). Thus, for completeness, we now give the complete formula for $S_\alpha (\rho_\beta)$ in the free fermion case. The four-point function on the cylinder of circumference $2\beta v$ is

$$\begin{align*} ({\rm free \; fermions})&\,\, \left\langle \mathcal{T}^\dagger_\alpha (0) \mathcal{T}_\alpha ({\rm i}\beta v) \mathcal{T}_\alpha (L_A) \mathcal{T}_\alpha^\dagger ({\rm i}\beta v+L_A) \right\rangle \nonumber \\ & \quad = \left[ \frac{\frac{2\beta v}{\pi } \sinh (\frac{\pi (L_A+{\rm i}\beta v)}{2\beta v}) \frac{2\beta v}{\pi } \sinh (\frac{\pi (L_A-{\rm i} \beta v)}{2\beta v}) }{\left(\frac{2\beta v}{\pi } \sin (\frac{\pi \beta v}{2\beta v}) \right)^2 \left(\frac{2\beta v}{\pi } \sinh (\frac{\pi L_A}{2\beta v}) \right)^2 } \right]^{2 \Delta_\alpha} = \left[ \frac{1 }{\frac{2\beta v}{\pi } \tanh (\frac{\pi L_A}{2\beta v}) } \right]^{4 \Delta_\alpha}, \nonumber \end{align*}$$

which gives the OSEE (see figure 4):

$$\begin{align} \label{eq:exact_thermal_free} ({\rm free \; fermions}) \qquad S_\alpha(\rho_\beta) = \frac{1}{3} \left(1+ \frac{1}{\alpha}\right) \log \left(\frac{2\beta v}{\pi a_0} \tanh (\frac{\pi L_A}{2\beta v}) \right). \nonumber \end{align} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

In the past years, there has been a lot of work on the problem of equilibration in closed quantum systems, and on generalized Gibbs ensembles (GGE), motivated by [41, 42]. In particular, many recent works (see, for instance, [43–53, 54]) have dealt with the problem of what charges exactly should be used as a basis set in the construction of the GGE: should they be local?, non-local?, quasi-local?, pseudo-local?, localized?, semi-local?, ultra-local? and so on (these terms are all properly defined in the previous list of references). These discussions are very advanced and concern subtle points about the thermodynamic Bethe Ansatz solution of specific models such as the XXZ chain or integrable field theories (for instance: what exactly is the minimal set of conserved charges whose expectation values are needed to uniquely specify a distribution of Bethe roots in the thermodynamic limit?). Such questions are far beyond the scope of the present paper. These recent works do inspire, however, the following more pedestrian question, which is directly related to our topic.

Assume, as usual in discussions about GGE, that the reduced density matrix of a subsystem A of size L_A (L_A small compared to the total size of the system) converges to some

$$\begin{align} \label{eq:H_GGE} \rho_{\rm GGE} = {\rm e}^{- H_{\rm GGE}} \nonumber \end{align} \tag{ 14 }$$

at large times (assume for instance that this convergence holds in L¹-norm, which then also implies convergence in L²-norm). As usual, one needs to be careful about possible revivals of the system in the definition of the convergence, for instance by averaging over time (possibly by excluding some time-window where the revival occurs); such details are well discussed in the standard treatments of GGE, and we will skip these subtle points here. We assume that $\rho_{\rm GGE}$ is properly defined, and we would like to get some idea of the locality of H_GGE not by focusing directly on the basis set for the space of conserved charges that should be used to write H_GGE, but by asking instead whether it would make sense to try to approximate $\rho_{\rm GGE} = {\rm e}^{-H_{\rm GGE}}$ by an MPO. We thus ask:

• (Q1)  
  Does the GGE density matrix $\rho_{\rm GGE}$ satisfy the operator area law?

In other words, does $\rho_{\rm GGE}$ behave morally like a thermal (Gibbs) density matrix, or is it spectacularly different from a Gibbs state, for instance with a logarithmic violation of the area law? (We note that the seemingly related question of whether $\rho_{\rm GGE}$ satisfies an area law for the mutual information has been addressed in [55]; there, a violation of the area law was found for a specific quench protocol. The question of the area law for the mutual information has also been investigated in open quantum systems, see for instance the recent [56–59].)

Since the construction of the GGE is in itself a hard task for interacting systems [60], it seems that a direct evaluation of the OSEE of $\rho_{\rm GGE}$ is out of reach at this point. However, it is of course easy to make progress in free fermion systems, which are already quite instructive for our purposes. Indeed, the free fermion case is arguably the worst when it comes to locality, in the sense that the conserved quantities that must be included in the construction of the GGE are the mode occupation themselves, which are less local than the ones that show up, say, in the XXZ chain (e.g. the 'quasi-local' ones [47]). So, we expect that, if a simple calculation in the free fermion case shows that $\rho_{\rm GGE}$ satisfies the operator area law, then it will be a good indication that the same will be true in other integrable systems. With this reasoning in mind, we now turn to a simple free-fermion calculation.

In this simple example, we will reach the conclusion that, for a global (translation invariant) quench to a critical point, as long as the initial state is the ground state of a gapped Hamiltonian, $\rho_{\rm GGE}$ satisfies the operator area law. This conclusion holds even though, perhaps surprisingly, H_GGE itself is not local² (in the sense of exponentially decaying terms). A similar observation was made previously for the mutual information in the transverse-field Ising chain: in that case, H_GGE is a sum of algebraically decaying terms [61], and yet the mutual information satisfies an area law according to the discussion in the conclusion of [55].

3.1.1. A free-fermion example.

For simplicity, we would like to treat an example of a quantum quench to the XX chain. For the initial state, a simple choice would be the Néel state as in figure 6 below, however since all occupation numbers in the Néel state are equal, one finds that $\rho_{\rm GGE} = 1$ in that case. So we can indeed conclude that $H_{\rm GGE}=0$ is local and that $\rho_{\rm GGE}=1$ satisfies the area law in that example; however it is too simplistic.

So, instead, we turn to a slightly less trivial example. We consider the following Hamiltonian for an infinitely long one-dimensional superconductor on a lattice,

$$\begin{align} H_{\delta,h} = \frac{1}{2} \int_{-\pi}^\pi \frac{{\rm d}k}{2\pi} \left(\begin{array}{@{}cc@{}} c^\dagger_k & c_{-k} \end{array} \right) \left(\begin{array}{@{}cc@{}} -\cos k - h & {\rm i} \delta \sin k \nonumber \\ -{\rm i} \delta \sin k & \cos k +h \end{array} \right) \left(\begin{array}{@{}c@{}} c_k \nonumber \\ c^\dagger_{-k} \end{array} \right), \nonumber \end{align} \tag{ 15 }$$

which is nothing but the XY chain in transverse field after a Jordan–Wigner transformation, see for instance [62, 63]. We focus on the following quench protocol: the system is prepared in the ground state of $H_{\delta, h}$, and at time t  =  0, we switch off the parameters δ and h,

$$\begin{align} H_{\delta, h} \quad {\rm at \;} t<0 \quad \longrightarrow \quad H_{0,0} \quad {\rm at \;} t>0. \nonumber \end{align} \tag{ 16 }$$

Notice that H_0,0 is the Hamiltonian of the XX chain. The ground state of $H_{\delta, h}$ is

$$\begin{align*} \left\vert{\psi_{\delta,h}}\right\rangle = \prod_{k>0} \left(\cos (\theta_k/2) + {\rm i} \sin (\theta_k/2) c^\dagger_{k} c^\dagger_{-k} \right) \left\vert{0}\right\rangle, \nonumber \end{align*}$$

where

$$\begin{align} \label{eq:theta_k} {\rm e}^{{\rm i} \theta_k} = \frac{h+\cos k + {\rm i} \delta \sin k }{\left\vert h+\cos k + {\rm i} \delta \sin k \right\vert }. \nonumber \end{align} \tag{ 17 }$$

The conserved charges of H_0,0 are the mode occupations at each wave vector k; when evaluated in the initial state, they give

$$\begin{align} n_k = \left\langle {\psi_{\delta,h}} \right\vert c^\dagger_k c_k \left\vert{\psi_{\delta,h}}\right\rangle = \sin (\theta_k/2)^2. \nonumber \end{align} \tag{ 18 }$$

The crucial property on which we are going to rely shortly is the (real-) analyticity of n_k as a function of k in the Brillouin zone. Notice that, assuming $\delta \neq 0$, n_k is always an analytic function of k, provided $\vert h\vert \neq 1$. This is simply because n_k is an analytic function of $\theta_k$, and $\theta_k$ itself is real-analytic in k as long as $h + \cos k$ is non-zero at k  =  0 and $k=\pi$, as can be seen from equation (17).

The GGE density matrix is

$$\begin{align} \rho_{\rm GGE} \propto {\rm e}^{-\sum_k \lambda_k c^\dagger_k c_k} \nonumber \end{align} \tag{ 19 }$$

where ${\rm e}^{-\lambda_k} = \frac{n_k}{1-n_k}$, such that the expectation values of the occupation modes are the same as in the initial state, ${\rm tr} [ c^\dagger_k c_k \rho_{\rm GGE} ] = \left\langle {\psi_{\delta, h}} \right\vert c^\dagger_k c_k \left\vert{\psi_{\delta, h}}\right\rangle$. Here the density matrix is the one of an infinite system; to obtain the one of a finite subsytem A, one can of course trace over all the degrees of freedom outside A.

Thus, we have $H_{\rm GGE} = -\sum_k \lambda_k c^\dagger_k c_k$, up to some irrelevant additional constant. The question of the locality of H_GGE follows from elementary results in Fourier analysis: H_GGE has exponentially decaying terms if $\lambda_k$ is real-analytic in k in the Brillouin zone. But we see that, if there is a point k where $n_k \rightarrow 0$ or $n_k \rightarrow 1$, then $\lambda_k$ is not analytic, and so H_GGE cannot be local in the sense of exponentially decaying terms. From equation (17), we see that this is what happens when $k \rightarrow 0$ and $k \rightarrow \pi$. So H_GGE cannot be local.

Let us now see why, despite the fact that H_GGE is not local (in the sense of exponentially decaying terms), $\rho_{\rm GGE}$ still satisfies the operator area law. To do that, one can use the trick that $\rho_{\rm GGE}$, which is an operator acting on a Hilbert space $\mathcal{H}$, can be viewed as a state in the space $\mathcal{H} \otimes \overline{\mathcal{H}}$ (see also figure 1 and the appendix for details on how to implement this in the free fermion case). We write the state $\left\vert{\rho_{\rm GGE} }\right\rangle$ in this larger space as

$$\begin{align} \label{eq:GGE_op_state} \rho_{\rm GGE} \propto \prod_k \left[ 1 + ({\rm e}^{-\lambda_k}-1) c^\dagger_k c_k \right] \quad \longrightarrow \quad \left\vert{\rho_{\rm GGE}}\right\rangle \propto \prod_k \left[ 1 + ({\rm e}^{-\lambda_k}-1) c^\dagger_k \tilde{c}^\dagger_k \right] \left\vert{0}\right\rangle, \nonumber \end{align} \tag{ 20 }$$

where the $\tilde{c}^\dagger_k$'s are new creation operators that anti-commute with all the c_k's, and $\left\vert{0}\right\rangle$ is the state annihilated by all the c_k's and $\tilde{c}_k$'s. The OSEE of $\rho_{\rm GGE}$ is, by definition, the 'usual' entanglement entropy of the state $\left\vert{\rho_{\rm GGE}}\right\rangle$ (see figure 1). Thus, we need to argue that the entanglement entropy of $\left\vert{\rho_{\rm GGE}}\right\rangle$ satisfies the area law. Since we are dealing with translation-invariant free fermions, this can be done by expressing the entanglement entropy in terms of a certain (block-) Toeplitz determinant, whose asymptotics depend on whether or not the symbol of the Toeplitz matrix possesses singularities [64, 65]. Since the symbol is related to n_k, it is possible to argue that the analyticity of n_k implies that the symbol itself is analytic, such that the area law then follows from the strong Szegö limit theorem. This line of arguments can presumably be made mathematically rigorous. However, for conciseness, we now take a slightly different route, which also sheds some light on the role played by the analyticity of n_k. Namely, we exhibit a local, gapped Hamiltonian acting on $\mathcal{H}\otimes \overline{\mathcal{H}}$, of which $\left\vert{\rho_{\rm GGE}}\right\rangle$ is the ground state. And we then rely on the standard fact that it is sufficient to imply the area law.

Such a Hamiltonian can be cooked up as follows. Consider the operator

$$\begin{align} \label{eq:super_ham} \sum_k \left(\begin{array}{@{}cc@{}} c^\dagger_k & \tilde{c}_k \end{array} \right) \left(\begin{array}{@{}cc@{}} \cos \alpha_k & -\sin \alpha_k \nonumber \\ -\sin \alpha_k & -\cos \alpha_k \end{array} \right) \left(\begin{array}{@{}c@{}} c_k \nonumber \\ \tilde{c}^\dagger_k \end{array} \right), \nonumber \end{align} \tag{ 21 }$$

which is a flat-band Hamiltonian, because the $2 \times 2$ matrix in the middle has eigenvalues  +1, −  1, independently of $\alpha_k$. The ground state of this operator is $\prod_k \left[ 1 + {\rm tan } (\alpha_k/2) c^\dagger_k \tilde{c}^\dagger_k \right] \left\vert{0}\right\rangle$. Comparing this to equation (20), we see that $\left\vert{\rho_{\rm GGE}}\right\rangle$ is the ground state of this flat-band hamiltonian if

$$\begin{align} {\rm tan } (\alpha_k/2) = {\rm e}^{-\lambda_k}-1 = \frac{1}{1-n_k} -2. \nonumber \end{align} \tag{ 22 }$$

This defines a function $n_k \mapsto \alpha_k (n_k)$ which is analytic in n_k. (One way of seeing that it is analytic is to write it as $\alpha_k (n_k) = 2 {\rm Im } [ \log (1-{\rm i} - (1-2i)n_k) ]$, where the branch cut of the $\log(.)$ is along the negative real axis; then observe that, for n_k real, $1-{\rm i} - (1-2i)n_k$ is never real and negative.) As a consequence, if n_k is an analytic function of k, so is $\alpha_k$. This is precisely the property we need in order to ensure that the operator (21) is local in real space (again, in the sense of exponentially decaying terms). So, in summary, provided that n_k is a real-analytic function of k in the Brillouin zone, $\left\vert{\rho_{\rm GGE}}\right\rangle$ is the ground state of a gapped, local Hamiltonian. As such, it must satisfy the area law. In contrast, in the critical case $\vert h\vert =1$, the occupation number n_k is not analytic as a function of k, therefore the operator (21) is not local in real space, leaving the possibility open that $\left\vert{\rho_{\rm GGE}}\right\rangle$ violates the area law.

So we have reached the conclusion that, in this simple example, as long as the initial XY Hamiltonian is gapped, $\rho_{\rm GGE}$ satisfies the operator area law; this is true even though H_GGE itself is usually not local (in the sense of exponentially decaying terms). The only case where $\rho_{\rm GGE}$ can possibly violate the operator area law is when the initial XY Hamiltonian is gapless, when n_k is no longer an analytic function of k.

We now focus on a quench from a translationally invariant initial state that is the ground state of some gapped hamiltonian, or global quench in the (now standard) terminology of Cardy and Calabrese [66–68]. One of the most important results they obtained about global quenches is the fact that the entanglement entropy of a subsystem A increases linearly in time, until it saturates to a value proportional to the length L_A. They interpreted this result in terms of pairs of quasi-particles being emitted from the initial state and generating a ballistic spreading of the entanglement (for a recent discussion, see also [69]). This result has dramatic consequences for numerical simulations, as it implies that any attempt at approximating the state at time t  >  0, $\left\vert{\psi(t)}\right\rangle$, with an MPS, will quickly face the exponential growth of the bond dimension. This unfortunately induces very strong limitations on MPS-based algorithms like time-dependent DMRG or TEBD.

However, in view of the discussion in the previous section, we could think of the following alternative way of simulating the dynamics after a global quench. If the system is not integrable, the reduced density matrix $\rho_A$ should converge to the thermal one, and thus should be approximable by an MPO. If the system is integrable, then $\rho_A$ is expected to converge to some $\rho_{\rm GGE}$, and we have given arguments for the approximability of this object by MPOs. So, as far as approximability is concerned, it makes no difference whether the dynamics is integrable or not: either way, the reduced density matrix at large times can be well approximated by an MPO. At very short times also, it is obviously approximable by an MPO (because the initial state itself, being the ground state of a gapped local Hamiltonian, is approximable by an MPS). So it seems that, both at short and large times, an MPO-based algorithm that would approximate the reduced density matrix $\rho_A(t)$—as opposed to an MPS approximation of the pure state $\left\vert{\psi(t)}\right\rangle$—would perform well. This raises the question³:

• (Q1-a)  
  Is it possible to overcome the exponential growth of the bond dimension after a global quench simply by trading the MPS for an MPO?

A naive sketch of an algorithm that would exploit this idea is the following (figure 5). First, start from an MPS approximation of the initial state $\left\vert{\psi_0}\right\rangle$, say with bond dimension χ, and use it to write the reduced density matrix $\rho_A$ as an MPO with bond dimension $\chi^2$. Tracing over the degrees of freedom in part B naturally gives two 'environments' E_L, E_R, that are two vectors, each of dimension $\chi^2$. Then one needs to design a protocol for updating both the bulk matrices (the V's in figure 5) and the environments E_L,R that define the MPO for $\rho_A(t)$. This could be some version of TEBD for the update of the V's; then, if one thinks of the simplest case where all the V's are equal at every time, V_j(t)  =  V(t), one could for instance use the left- and right-eigenvectors (with maximal eigenvalue) of V(t) as updates for E_L,R(t).

Zoom In Zoom Out Reset image size

More elaborate versions of such algorithms have been investigated much more thoroughly in the literature, for instance in [20].

Sadly enough, the answer to (Q1-a) is negative; the reason is the following. If one calculates the OSEE of the reduced density matrix $\rho_A$, for a bipartition of the subsystem A in the form of two consecutive intervals $A = A_1 \cup A_2$ of lengths x and L_A  −  x respectively (figure 6), then, as expected, the OSEE is small both at short and at large times. This supports the fact that both the short-time and large-time behaviors can be captured by an MPO with small bond dimension. However, the OSEE blows up linearly in the transient regime $t \lesssim \frac{L_A}{2 v}$, implying that the bond dimension that is needed to approximate $\rho_A(t)$ in that time window does blow up exponentially. This is illustrated in figure 6, where numerical results for a quench in the XX chain from a Néel initial state are displayed. In fact, the linear growth of the OSEE at times $t \leqslant \frac{x}{v}$ can be traced back to equation (6): at short times, the reduced density matrix is still very close to the one of the pure state, and one does not (yet) gain much by approximating $\rho_A (t)$ instead of $\left\vert{\psi(t)}\right\rangle$; this trick becomes efficient only at later times.

Zoom In Zoom Out Reset image size

Finally, we turn to the CFT calculation. As in previous sections, we start by rewriting $S_\alpha (\rho_A(t))$ in terms of swap operators for α integer. From equation (7), we see that $S_\alpha (\rho_A(t)) = \frac{1}{1-\alpha} \log r(t)$, where r(t) is the ratio of two expectation values in a system with $2\alpha$ replicas,

$$\begin{align} r(t) \equiv \frac{\left\langle{\psi(t)}\right\vert ^{\otimes 2 \alpha} \; \mathcal{S}^{A_1}_{(1,\alpha+2)(2,\alpha+3)\dots (\alpha-1,2\alpha) (\alpha,\alpha+1)} \mathcal{S}^{A_2}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)}\; \left\vert{\psi(t)}\right\rangle ^{\otimes 2 \alpha} }{\left\langle{\psi(t)}\right\vert ^{\otimes 2 \alpha} \mathcal{S}^A_{(1,\alpha+1)(2,\alpha+2)\dots (\alpha,2\alpha)} \left\vert{\psi(t)}\right\rangle ^{\otimes 2 \alpha} }. \nonumber \end{align} \tag{ 23 }$$

This expression is very similar to the one appearing in the calculation of the negativity in [70]; the only difference is that the permutations involved are not the same. It is then straightforward to adapt the calculation of [70] to the OSEE. We start by making a crucial assumption about the initial state $\left\vert{\psi_0}\right\rangle$: we assume that it is of the form

$$\begin{align} \label{eq:CC_initial} \left\vert{\psi(t)}\right\rangle = {\rm e}^{- \frac{\lambda}{v} H} \left\vert{b}\right\rangle, \nonumber \end{align} \tag{ 24 }$$

where $\left\vert{b}\right\rangle$ is a conformally invariant boundary state, H is the Hamiltonian of the CFT, and $\lambda >0$ is some UV cutoff, known as the extrapolation length in the field theory approach to surface critical phenomena [71]. One simple way to understand the appearance of this length is as follows (this argument appeared in [72, 73] and was later adapted by Cardy [74] to quantum quenches): since the initial state $\left\vert{\psi_0}\right\rangle$ is assumed to have short-range correlations, it must flow (under the RG) to a conformal boundary state; therefore $\left\vert{\psi_0}\right\rangle$ itself can be written as the RG fixed point $\left\vert{b}\right\rangle$ perturbed by some irrelevant operators $\phi_p$, $\left\vert{\psi_0}\right\rangle \propto {\rm e}^{- \int d x \sum_p \lambda_p \phi_p(x)} \left\vert{b}\right\rangle$. One then needs to discuss what irrelevant operators are allowed in this expression. The most generic one, which is also often the least irrelevant one, is the stress-tensor itself, which comes with a coefficient λ homogeneous to a length. But the effect of perturbing the boundary condition by the stress-tensor is precisely to shift the position of the boundary. Therefore, taking $\left\vert{\psi_0}\right\rangle$ in the form (24) is equivalent to neglecting all perturbations but the one due to the stress-tensor. Of course, in principle, one should take any possible perturbations into account (as sketched in [72] in a different context, and by Cardy in the context of quantum quenches in [74]), however this is much more difficult in practice, so we will omit all other perturbations and proceed from equation (24), as in most references dealing with global quenches in CFT, including [70].

The state at time t is then

$$\begin{align*} \left\vert{\psi(t)}\right\rangle = {\rm e}^{- {\rm i} H t} {\rm e}^{- \frac{\lambda}{v} H} \left\vert{b}\right\rangle, \nonumber \end{align*}$$

and, as usual in this type of calculations, it is useful to switch to imaginary time: we replace ivt by y, or equivalently

$$\begin{align} \left\vert{\psi(t)}\right\rangle \rightarrow {\rm e}^{- \frac{\lambda+y}{v} H} \left\vert{b}\right\rangle \qquad {\rm and} \qquad \left\langle{\psi(t)}\right\vert \rightarrow \left\langle{b}\right\vert {\rm e}^{- \frac{\lambda-y}{v} H}. \nonumber \end{align} \tag{ 25 }$$

Then, for $-\lambda <y<\lambda$, the ratio above becomes a ratio of correlation functions of three twist operators in an infinite strip (figure 7) of width $2\lambda$, parametrized by $(x, y) \in \mathbb{R} \times [- \lambda, \lambda]$,

$$\begin{align*} r(t) \rightarrow \frac{(\left\langle{b}\right\vert {\rm e}^{- \frac{\lambda-y}{v} H})^{\otimes 2 \alpha} \; \mathcal{S}^{A_1}_{(1,\alpha+2)(2,\alpha+3)\dots (\alpha-1,2\alpha) (\alpha,\alpha+1)} \mathcal{S}^{A_2}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)}\; ({\rm e}^{- \frac{\lambda+y}{v} H} \left\vert{b}\right\rangle)^{\otimes 2 \alpha} }{(\left\langle{b}\right\vert {\rm e}^{- \frac{\lambda-y}{v} H})^{\otimes 2 \alpha} \mathcal{S}^A_{(1,\alpha+1)(2,\alpha+2)\dots (\alpha,2\alpha)} ({\rm e}^{- \frac{\lambda+y}{v} H} \left\vert{b}\right\rangle)^{\otimes 2 \alpha} } = \frac{G_3}{G_2} \nonumber \\ \nonumber \end{align*}$$

with

$$\begin{align*} G_2 & \equiv \left\langle \mathcal{T}_{(1,\alpha+1)(2+\alpha+2)\dots (\alpha,2\alpha)} ({\rm i} y) \mathcal{T}_{(1,\alpha+1)(2+\alpha+2)\dots (\alpha,2\alpha)} (L_A+ {\rm i} y) \right\rangle \nonumber \\ G_3 & \equiv \left\langle \mathcal{T}_{(1,\alpha+2)(2,\alpha+3)\dots (\alpha-1,2\alpha) (\alpha,\alpha+1)}({\rm i} y) \mathcal{T}_{(1,2,3,\dots,\alpha)(\alpha+1,2\alpha,2\alpha-1,\dots,\alpha+2)}(x+{\rm i} y) \mathcal{T}_{(1,\alpha+1)(2,\alpha+2)\dots(\alpha,2\alpha)}(L_A+{\rm i} y) \right\rangle. \nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

The strip is conformally mapped to the upper half-plane by $x+{\rm i} y \mapsto \zeta = {\rm i} {\rm e}^{\frac{\pi}{2 \lambda} (x+ iy)}$; it is convenient to work with the coordinates $\zeta_1$, $\zeta_2$, $\zeta_3$ in the upper half-plane, and their cross-ratios,

$$\begin{align} \left\{\begin{array}{@{}rcl@{}} \zeta_1 &=& {\rm i} {\rm e}^{\frac{\pi}{2 \lambda} ({\rm i} y) } \nonumber \\ \zeta_2 &=& {\rm i} {\rm e}^{\frac{\pi}{2 \lambda} (L_A+ {\rm i} y) } \nonumber \\ \zeta_3 &=& {\rm i} {\rm e}^{\frac{\pi}{2 \lambda} (x+ {\rm i} y) } \end{array} \right. ; \qquad \eta_{ij} = \frac{(\zeta_i - \zeta_j)(\overline{\zeta}_i-\overline{\zeta}_j)}{(\zeta_i- \overline{\zeta}_j) (\overline{\zeta}_i- \zeta_j) }. \nonumber \end{align} \tag{ 26 }$$

Then we have

$$\begin{align*} G_2 &= \left\vert \frac{{\rm d} \zeta_1}{{\rm d} z_1} \right\vert ^{\alpha \Delta_2} \left\vert \frac{{\rm d} \zeta_2}{{\rm d} z_2} \right\vert ^{\alpha \Delta_2} \frac{1}{\vert \zeta_1-\overline{\zeta}_1\vert ^{\alpha \Delta_2} \vert \zeta_2-\overline{\zeta}_2\vert ^{\alpha \Delta_2}} \frac{1}{\eta_{12}^{\alpha \Delta_2}} \; f_2(\eta_{12}, \overline{\eta}_{12}), \nonumber \\ G_3 &= \left\vert \frac{{\rm d} \zeta_1}{{\rm d} z_1} \right\vert ^{\alpha \Delta_2} \left\vert \frac{{\rm d} \zeta_2}{{\rm d} z_2} \right\vert ^{\alpha \Delta_2} \left\vert \frac{{\rm d} \zeta_3}{{\rm d} z_3} \right\vert ^{2 \Delta_\alpha} \frac{1}{\vert \zeta_1-\overline{\zeta}_1\vert ^{\alpha \Delta_2} \vert \zeta_2-\overline{\zeta}_2\vert ^{\alpha \Delta_2} \vert \zeta_3-\overline{\zeta}_3\vert ^{2 \Delta_\alpha}} \left(\frac{\eta_{12}^{2 \Delta_\alpha - 2 \alpha \Delta_2} }{\eta_{13}^{2 \Delta_\alpha} \eta_{23}^{2 \Delta_\alpha} } \right)^{\frac{1}{2}} \; f_3(\{\eta_{ij} \}, \{\overline{\eta}_{ij} \}). \nonumber \end{align*}$$

The Jacobians in front of both expressions come from the conformal mapping $z\mapsto \zeta$, while the rest is the correlation function in the upper half-plane. At the points $\zeta_1$ and $\zeta_2$, the twist operators have scaling dimension $\alpha \Delta_2$, since they correspond to permutation with α cycles of length 2. At the point $\zeta_3$, the permutation has two cycles of length α, therefore the scaling dimension is $2 \Delta_\alpha$. There is of course not a unique way to write formulas for G₂ and G₃, since there is some arbitrariness in the choice of the functions of the anharmonic ratios f₂ and f₃. As in [70], we have chosen our expressions such that f_2,3 always go to a finite constant whenever one of the anharmonic ratios $\eta_{ij}$ approaches 0 or 1, as can be seen from the fusion rules of the twist operators. The ratio of interest is then

$$\begin{align*} \frac{G_3}{G_2} &= \left(\left\vert \frac{{\rm d} \zeta_3}{{\rm d} z_3} \right\vert \frac{1}{\vert \zeta_3-\overline{\zeta}_3 \vert } \left(\frac{\eta_{12} }{\eta_{13}\eta_{23} }\right)^{\frac{1}{2}} \right)^{2 \Delta_\alpha} \frac{f_3(\eta_{12}, \eta_{23},\eta_{13}, \overline{\eta}_{12}, \overline{\eta}_{23},\overline{\eta}_{13}) }{f_2(\eta_{12}, \overline{\eta}_{12}) } \nonumber \\ &= \left(\frac{\pi}{4 \lambda} \frac{1}{\cos \frac{\pi y}{2 \lambda} } \left(\frac{\eta_{12} }{\eta_{13}\eta_{23} }\right)^{\frac{1}{2}} \right)^{2\Delta_\alpha} \frac{f_3(\eta_{12}, \eta_{23},\eta_{13}, \overline{\eta}_{12}, \overline{\eta}_{23},\overline{\eta}_{13}) }{f_2(\eta_{12}, \overline{\eta}_{12}) }. \nonumber \end{align*}$$

The anharmonic ratios are equal to

$$\begin{align} \eta_{ij} = \frac{2 \sinh^2 \left(\frac{\pi}{4 \lambda} (x_i-x_j) \right)}{\cos \left(\frac{\pi y}{\lambda} \right) + \cosh \left(\frac{\pi}{2\lambda} (x_i-x_j) \right) } \nonumber \end{align} \tag{ 27 }$$

and when we Wick-rotate back $y\rightarrow {\rm i} vt$, they become

$$\begin{align*} \eta_{ij} =&~ \frac{2 \sinh^2 \left(\frac{\pi}{4 \lambda} (x_i-x_j) \right)}{\cosh \left(\frac{\pi v t}{\lambda} \right) + \cosh \left(\frac{\pi}{2\lambda} (x_i-x_j) \right) } \nonumber \\ &\mathop{\simeq}\limits_{x_i-x_j, vt \gg \lambda} \frac{1 }{1+ {\rm e}^{\frac{\pi }{\lambda} (vt- \frac{1}{2} \vert x_i-x_j \vert)} } \nonumber \\ & \sim \left\{\begin{array}{@{}ccl@{}} 1 & {\rm if} & \frac{1}{2} \vert x_i - x_j\vert - vt \gg \lambda >0 \nonumber \\ {\rm e}^{-\frac{\pi }{2 \lambda} (2vt- \vert x_i-x_j \vert)} &{\rm if} & vt - \frac{1}{2} \vert x_i - x_j\vert \gg \lambda >0. \end{array} \right. \nonumber \end{align*}$$

To study the regime where $\left\vert vt - \frac{1}{2} \vert x_i - x_j\vert \right\vert$ is of order λ, we would need the full dependence on the aspect ratios, which means that we would need to know the functions f₂ and f₃. In general, this is difficult. But there would be little to gain anyway, since we are most interested in results in the thermodynamic limit, where all length scales involved are much larger than λ. Under this assumption, f₂ and f₃ can be replaced by constants, and we see that

$$\begin{align} \frac{G_3}{G_2} \simeq \left\{\begin{array}{@{}ccl@{}} \left(\frac{\pi}{2 \lambda} {\rm e}^{-\frac{\pi v t}{2 \lambda}} \left(\frac{1 }{1 \times 1 }\right)^{\frac{1}{2}} \right)^{2\Delta_\alpha} \times {\rm const.} \quad &{\rm if} & vt < \frac{1}{2} {\rm min} (x, L_A -x) \nonumber \\ \left(\frac{\pi}{2 \lambda} {\rm e}^{-\frac{\pi v t}{2 \lambda}} \left(\frac{1}{{\rm e}^{-\frac{\pi }{2\lambda} (2vt- {\rm min} (x, L_A -x))} }\right)^{\frac{1}{2}} \right)^{2\Delta_\alpha} \times {\rm const.} \quad &{\rm if} & \frac{1}{2} {\rm min} (x, L_A -x) < vt < \frac{1}{2} {\rm max} (x, L_A -x) \nonumber \\ \left(\frac{\pi}{2 \lambda} {\rm e}^{-\frac{\pi v t}{2 \lambda}} \left(\frac{1}{{\rm e}^{-\frac{\pi }{2\lambda} (4vt- L_A)} }\right)^{\frac{1}{2}} \right)^{2\Delta_\alpha} \times {\rm const.} \quad &{\rm if} & \frac{1}{2} {\rm max} (x, L_A -x) < v t < L_A \nonumber \\ \left(\frac{\pi}{2 \lambda} {\rm e}^{-\frac{\pi v t}{2 \lambda}} \left(\frac{{\rm e}^{-\frac{\pi }{2\lambda} (2vt- L_A)} }{{\rm e}^{-\frac{\pi }{2\lambda} (4vt- L_A)} }\right)^{\frac{1}{2}} \right)^{2 \Delta_\alpha} \times {\rm const.} \quad &{\rm if} & L_A < v t. \end{array} \right. \nonumber \end{align} \tag{ 28 }$$

We are finally ready to write the CFT result for the OSEE, obtained by taking the logarithm of $G_3/G_2$, with the prefactor $1/(1-\alpha)$,

$$\begin{align} \label{eq:OSEE_global_quench} S_\alpha(\rho_A(t)) \simeq \frac{c}{6} \left(1+\frac{1}{\alpha}\right) \times \left\{\begin{array}{@{}ccl@{}} \frac{\pi v t}{2 \lambda} + {\rm const.} \quad &{\rm if} & vt < \frac{1}{2}{\rm min} (x, L_A -x) \nonumber \\ \frac{\pi }{4\lambda} {\rm min} (x, L_A -x) + {\rm const.} \quad &{\rm if} & \frac{1}{2} {\rm min} (x, L_A -x) < vt < \frac{1}{2} {\rm max} (x, L_A -x) \nonumber \\ \frac{\pi L_A}{4\lambda} - \frac{\pi v t}{2\lambda} + {\rm const.} \quad &{\rm if} & \frac{1}{2} {\rm max} (x, L_A -x) < v t < L_A \nonumber \\ 0+{\rm const.} \quad &{\rm if} & L_A < v t. \end{array} \right. \nonumber \end{align} \tag{ 29 }$$

The additive constant needs not be the same in the four lines, and the possible shift involved reflects the role of f₃ and f₂; however, this is always a subleading effect. A comparison between this formula and numerical results in the XX chain is displayed in figure 6 for the case x  =  L_A/2. The agreement is qualitatively good. At t  <  L_A/(2v), the numerical result is clearly compatible with the linear growth we just found. However, for t  >  L_A/(2v) the discrepancy between the numerics and the analytical prediction is larger. There are at least two possible explanations for that. One explanation could be that the initial state (24) we considered is not a good CFT representation of the Néel state for the XX chain, and that other perturbations should be included, as discussed above. Another explanation could be that the subleading effects due to the functions f₂ and f₃ are still quite large for the system sizes we probe numerically. These questions lie beyond the scope of the present paper and will be investigated more thoroughly elsewhere.

In conclusion, our CFT calculation—which focuses on the leading behavior only and could possibly be improved to grasp also subleading effects—confirms that the OSEE initially grows linearly in time, and therefore blows up in the transient regime, even though it decays at later times. This is in perfect agreement with the discussion above.

We consider once again an infinitely long system described by a CFT with a Hamiltonian H, with the corresponding evolution operator

$$\begin{align} U(t) = {\rm e}^{- {\rm i} H t}. \nonumber \end{align} \tag{ 30 }$$

We want to calculate the OSEE for the bipartition $A \cup B = (-\infty, 0] \cup [ 0, +\infty)$. At time t  =  0, U(t) is the identity operator, whose OSEE is obviously zero; at later times, U(t) is a non-trivial operator with non-zero OSEE. The question is: how fast does it increase as a function of time?

In CFT, this question is easily answered as follows. We rely again on equation (5), assuming α integer $\geqslant 2$. Here, equation (5) reads

$$\begin{align*} S_\alpha (U(t)) = \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} [({\rm e}^{{\rm i} H t})^{\otimes \alpha} \cdot \mathcal{S}_{(1,\alpha,\alpha-1,\dots,2)} \cdot ({\rm e}^{-{\rm i} H t})^{\otimes \alpha} \cdot \mathcal{S}_{(1,2,\dots,\alpha)}] }{({\rm tr} [ {\rm e}^{{\rm i} H t} {\rm e}^{-{\rm i} H t} ])^{\alpha} } \right). \nonumber \end{align*}$$

Notice that the denominator inside the logarithm is trivially equal to one; the numerator cannot be evaluated that easily though. As in the previous, it is convenient to switch to imaginary time, and to start by calculating the quantity

$$\begin{align*} \frac{{\rm tr} [({\rm e}^{-(\frac{\beta}{2}-\frac{y}{v}) H})^{\otimes \alpha} \cdot \mathcal{S}_{(1,\alpha,\alpha-1,\dots,2)} \cdot ({\rm e}^{-(\frac{\beta}{2}+\frac{y}{v}) H})^{\otimes \alpha} \cdot \mathcal{S}_{(1,2,\dots,\alpha)}] }{({\rm tr} [ {\rm e}^{- \frac{\beta}{2} H} {\rm e}^{- \frac{\beta}{2} H} ])^{\alpha} }, \nonumber \end{align*}$$

and at the end of the calculation, do the Wick rotation $y \rightarrow {\rm i} v t$ and take the limit $\beta \rightarrow 0^+$. We will keep β finite though, as it conveniently plays the role of a UV cutoff. The big advantage of looking at this imaginary-time version of the ratio is that, as long as $\vert y\vert < \beta v/2$, it is again (exactly like in section 2.2) a ratio of two partition functions of the CFT on a cylinder with α sheets, see figure 8, that boils down to a two-point function of twist operators. Indeed, if we parametrize the cylinder by the complex coordinates $x+{\rm i} y$, $(x, y) \in \mathbb{R} \times [-\beta v/2, \beta v/2]$, then the two swap operator $\mathcal{S}_\alpha$ can be replaced by a pair of twist operators at $- {\rm i} \beta v/2$ and ${\rm i} y$ respectively. There are also two operators at the left end of the cylinder (at infinity), but those two can be fused together to the identity; so they do not affect the result in the end. We are left with the two-point function

$$\begin{align*} & \frac{{\rm tr} [({\rm e}^{-(\frac{\beta}{2}-\frac{y}{v}) H})^{\otimes \alpha} \cdot \mathcal{S}_{(1,\alpha,\dots,2)} \cdot ({\rm e}^{-(\frac{\beta}{2}+\frac{y}{v}) H})^{\otimes \alpha} \cdot \mathcal{S}_{(1,2,\dots, \alpha)} ] }{({\rm tr} [ {\rm e}^{- \beta H} ])^{\alpha} } \nonumber \\ &= \left\langle \mathcal{T}_{(1,\alpha,\dots,2)} (-{\rm i} \beta v/2) \mathcal{T}_{(1,2,\dots,\alpha)} ({\rm i} y)\right\rangle \propto \frac{1}{[ \frac{\beta v }{\pi } \sin (\frac{\pi (\frac{\beta v}{2}+y)}{\beta v}) ]^{2 \Delta_\alpha}}, \nonumber \end{align*}$$

up to some non-universal constant factor that is homogeneous to a length to the power $2 \Delta_\alpha$, because the ratio of traces must be dimensionless. Then, Wick-rotating back to $y \rightarrow {\rm i} v t$, and adjusting the constant factor such that $S_\alpha (U(t))$ vanishes at t  =  0, we find

$$\begin{align*} S_\alpha (U_t) =& \frac{c}{6}\left(1+ \frac{1}{\alpha}\right) \log \left(\cosh \left(\frac{\pi t }{\beta} \right) \right) \nonumber \\ & \mathop{\simeq}\limits_{t \gg \tau } \frac{c}{6}\left(1+ \frac{1}{\alpha}\right) \frac{\pi t }{\beta}. \nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

The OSEE increases linearly in time. In other words, if one wants to approximate U(t) by an MPO, then one needs a bond dimension that blows up exponentially with time. This conclusion, of course, in not surprising, and it is intimately tied to the famous linear growth of the entanglement entropy after a global quench, first pointed out by Calabrese and Cardy [66].

In fact, there are two good reasons why the OSEE of the evolution operator must be closely related to the entanglement entropy after a global quench from an initial state with short-range correlations $\left\vert{\psi_0}\right\rangle$. First, it is because the entanglement of $U(t) \left\vert{\psi_0}\right\rangle$ obviously does not come from $\left\vert{\psi_0}\right\rangle$ itself, so it must come from U(t); in other words, since $\left\vert{\psi_0}\right\rangle$ is well approximated by an MPS with small bond dimension, the only explanation for the fact that the entanglement entropy of $U(t) \left\vert{\psi_0}\right\rangle$ blows up must be that the bond dimension of U(t) itself blows up. Second, upon using the 'operator-folding' trick shown in figure 1(c), namely replacing the operator $U(t) \in {\rm End}(\mathcal{H})$ by the state $\left\vert{U(t)}\right\rangle \in \mathcal{H} \otimes \overline{\mathcal{H}}$, on sees that U(t) itself can be viewed as resulting from a global quench,

$$\begin{align} \left\vert{U(t)}\right\rangle = {\rm e}^{-{\rm i} (H \otimes 1 - 1 \otimes H) t} \left\vert{{\rm Id}}\right\rangle. \nonumber \end{align} \tag{ 31 }$$

This implies that the linear growth of the OSEE of U(t) is valid far beyond the range of applicability of CFT: all systems with ballistic spreading of correlations are known to have an entanglement entropy that grows linearly after a global quench (see e.g. the recent discussion [69]), so they must have an evolution operator U(t) whose OSEE increases linearly in time. Again, this should not be surprising to most readers, however the point we want to emphasize here is that this is really a property of the evolution operator (and therefore of the Hamiltonian itself), which does not have anything to do with any particular quench protocol.

For the same reason, in phases where there is no ballistic spreading of correlations, namely in localized phases, the OSEE of the evolution must mimic the generic behavior of the entanglement entropy after a global quench. In free fermion systems that are (Anderson) localized, the OSEE of U(t) must be bounded, while in many-body localized phases, it must grow logarithmically [33, 77, 78]. In the (Anderson-localized) free-fermion case, this is simply because the Hamiltonian is diagonalized in a basis of orbitals that are exponentially localized. For instance, if we consider the following simple model

$$\begin{align} H = - \frac{1}{2} \sum_x [c^\dagger_x c_{x+1} + {\rm h.c.}] + \sum_x V(x) c^\dagger_x c_x, \nonumber \end{align} \tag{ 32 }$$

for some random onsite potential V(x), then it is always true that it can be put in diagonal form

$$\begin{align} H = \sum_x \varepsilon_x a^\dagger_x a_x, \nonumber \end{align} \tag{ 33 }$$

with some random energies $\varepsilon_x$, and where the new creation/annihilation modes are exponentially localized around the original ones,

$$\begin{align} u_{x,x^{\prime}} \equiv \left\langle{0}\right\vert a(x) c^\dagger(x^{\prime}) \left\vert{0}\right\rangle \sim {\rm e}^{-\vert x-x^{\prime}\vert /\xi}. \nonumber \end{align} \tag{ 34 }$$

$\xi > 0$ is the localization length, it depends on the strength of the random potential. The amplitudes $u_{x, x^{\prime}}$ satisfy the normalization condition $\sum_{x^{\prime\prime}} u_{x, x^{\prime\prime}} u^*_{x^{\prime}, x^{\prime\prime}} = \delta_{x, x^{\prime}}$. The Hilbert space $\mathcal{H}$ is the Fock space generated by the $c^\dagger$'s; it is useful to introduce a copy of this Fock space, noted $\mathcal{H}_{f^\dagger}$, generated by other fermion modes $f^\dagger$'s. The $c^\dagger$'s and $f^\dagger$'s do not act on the same space, but it is nevertheless convenient to write expressions that involve both kinds of operators; we decide, by convention, that the $c^\dagger$'s and the $f^\dagger$'s always anticommute: $\{c^\dagger_x, f_{x^{\prime}} \} =\{c^\dagger_x, f^\dagger_{x^{\prime}} \} =0$. Then we can write the following map $W: \mathcal{H} \rightarrow \mathcal{H}_{f^\dagger}$,

$$\begin{align} W = \exp\left(\sum_{x,x^{\prime}} u_{x,x^{\prime}} f^\dagger_x c_{x^{\prime}} \right), \nonumber \end{align} \tag{ 35 }$$

which takes a localized orbital $a^\dagger_x = \sum u^*_{x, x^{\prime}} c^\dagger_{x^{\prime}}$ and maps it to $f^\dagger_x$. W is local in the sense that the coefficients $u_{x, x^{\prime}}$ decay exponentially with distance; it also satisfies $W^\dagger W = 1$, so it is in fact a local unitary transformation. Because W is local, it must satisfy the operator area law: it has a bounded OSEE. The evolution operator now takes the form

$$\begin{align} ({\rm Anderson \; localized}) \qquad {\rm e}^{-{\rm i} H t} = W^\dagger \cdot \left(\bigotimes_{x} {\rm e}^{- {\rm i} t \varepsilon_x f^\dagger_x f_x } \right) \cdot W. \nonumber \end{align} \tag{ 36 }$$

The point is that the operator in the middle, which is the only piece that is time-dependent in this expression, is a product operator, therefore it cannot generate any entanglement. So, in an Anderson-localized system, the evolution operator itself satisfies the operator area law.

In many-body localized phases, the story is similar, except that the diagonalized evolution operator does not take the form of a tensor product anymore. Instead, it is a sum of commuting terms that couple the localized orbitals (the $f^\dagger$'s) at arbitrary distances, with exponentially decaying amplitudes,

$$\begin{align} ({\rm many{\rm{-}}body \; localized}) \qquad {\rm e}^{-{\rm i} H t} = W^\dagger \cdot \exp \left(- {\rm i} t \left[ \sum_x \varepsilon_x f^\dagger_x f_x + \sum_{x,x^{\prime}} v_{x,x^{\prime}} f^\dagger_x f_x f^\dagger_{x^{\prime}} f_{x^{\prime}} + \dots \right] \right) \cdot W. \nonumber \end{align} \tag{ 37 }$$

The latter are responsible for a dephasing mechanism, identified in [78], which ultimately generates a logarithmic growth of the entanglement; these arguments can be repeated here to show that the OSEE must grow logarithmically⁵.

Let us start by summarizing the observations of [2]. For simplicity, we consider only the case of a Hamiltonian of the form $H = \int \frac{{\rm d}k}{2\pi} \varepsilon(k) c^\dagger(k) c(k)$, which conserves the number of fermions. But the statements also hold for Hamiltonians with pairing terms (i.e. $c^\dagger c^\dagger$ and cc).

5.1.1. The fermion creation/annihilation operator is an exact MPO.

The first observation [2] is that

$$\begin{align*} c^\dagger_x (t) = \sum_{x^{\prime} \in \mathbb{Z}} K(x-x^{\prime},t) c^\dagger_{x^{\prime}} \nonumber \end{align*}$$

for some kernel $K(x-x^{\prime}, t)$. For the XX chain, K can be written with a Bessel function, but in general the precise form of $K(x-x^{\prime}, t)$ does not matter. What is important is to observe that this expression of $c^\dagger_x (t)$ can be interpreted as an MPO with bond dimension 2. Consequently, all finite products of fermion creation/annihilation operators, of the form $c_{x_1}^\dagger \dots c_{x_p}^\dagger c_{x^{\prime}_1} \dots c_{x^{\prime}_q}$, are exact MPOs with finite bond dimension (the bond dimension is at most 2^p+q). Of course, this implies that all such operators satisfy the operator area law, because their OSEE is obviously not larger than $(\,p+q) \log 2$.

5.1.2. The OSEE of a Jordan–Wigner string increases logarithmically in time.

The second observation [2], more interesting and intriguing, is that the (endpoint of the) Jordan–Wigner (JW) string,

$$\begin{align} {\rm JW}_x \equiv \prod_{x^{\prime} \leqslant x} (-1)^{c_{x^{\prime}}^\dagger c_{x^{\prime}}} = \prod_{x^{\prime} \leqslant x} (1 - 2 c_{x^{\prime}}^\dagger c_{x^{\prime}}), \nonumber \end{align} \tag{ 38 }$$

when represented in Heisenberg picture, namely

$$\begin{align} {\rm JW}_x(t) = \prod_{x^{\prime} \leqslant x} [ 1 - 2 c_{x^{\prime}}^\dagger (t) c_{x^{\prime}}(t) ], \nonumber \end{align} \tag{ 39 }$$

has an OSEE that grows logarithmically at large time t. The prefactor of the logarithm of t is always a rational number (for $\alpha=1$); Prosen and Pižorn conjectured a formula for this prefactor in [34]; we shall derive this prefactor for the XX chain below.

5.1.3. All local operators have an OSEE that increases at most logarithmically.

We now turn to the explicit calculation of the OSEE of a JW string in the XX chain. Without loss of generality, the endpoint of the string can be placed at the origin, so we are interested in the OSEE of ${\rm JW}_0 (t)$. We proceed as follows (see also the appendix for more details about how to calculate the OSEE in free fermion systems). We start by turning the operator into a state

$$\begin{align} {\rm JW}_0(t) \longrightarrow \left\vert{{\rm JW}_0 (t)}\right\rangle \equiv \prod_{x \leqslant 0} [ 1 - 2 c_x^\dagger(t) \tilde{c}^\dagger_x(t) ] \left\vert{0}\right\rangle, \nonumber \end{align} \tag{ 40 }$$

namely we replace c(t) by some new creation operator $\tilde{c}^\dagger (t)$, which anticommutes with all the c's. $\left\vert{0}\right\rangle$ is the state that is annihilated by all the c's and all the $\tilde{c}$'s. As sketched in figure 1(c), the main interest of rewriting operators as states in an enlarged Hilbert space (here a fermion Fock space) resides in the fact that the OSEE becomes the usual entanglement entropy; we will make use of this shortly.

It is convenient to think of this state as coming from the time-evolution

$$\begin{align} \left\vert{{\rm JW}_0(t)}\right\rangle = {\rm e}^{-{\rm i} (H\otimes \tilde{1} - 1\otimes \tilde{H}) t} \left\vert{{\rm JW}_0(0)}\right\rangle, \nonumber \end{align} \tag{ 41 }$$

where $H\otimes \tilde{1} - 1\otimes \tilde{H}$ acts as

$$\begin{align} \begin{array}{@{}lll@{}} \left[H\otimes \tilde{1} - 1\otimes \tilde{H}, c^\dagger(k) \right] = \varepsilon(k) c^\dagger(k) & \qquad \quad& \left[H\otimes \tilde{1} - 1\otimes \tilde{H}, \tilde{c}^\dagger(k) \right] = - \varepsilon(k) \tilde{c}^\dagger(k) \nonumber \\ \nonumber \\ \left[H\otimes \tilde{1} - 1\otimes \tilde{H}, c(k) \right] = -\varepsilon(k) c(k) & \qquad \quad& \left[H\otimes \tilde{1} - 1\otimes \tilde{H}, \tilde{c}(k) \right] = \varepsilon(k) \tilde{c}(k). \end{array} \nonumber \end{align} \tag{ 42 }$$

Next, we perform the following Bogoliubov transformation,

$$\begin{align} \label{eq:bogoliubov} \left(\begin{array}{@{}c@{}} b^\dagger_x \nonumber \\ d_x \end{array} \right) \equiv \left(\begin{array}{@{}cc@{}} \frac{2}{\sqrt{5}} & - \frac{1}{\sqrt{5}} \nonumber \\ \frac{1}{\sqrt{5}} & \frac{2}{\sqrt{5}} \end{array} \right) \left(\begin{array}{@{}c@{}} c^\dagger_x \nonumber \\ \tilde{c}_x \end{array} \right). \nonumber \end{align} \tag{ 43 }$$

This transformation does not do anything to the Hamiltonian $H\otimes \tilde{1} - 1\otimes \tilde{H}$, which already was (and still is) in diagonal form. However, it simplifies the initial state:

$$\begin{align} \left\vert{{\rm JW}_0(0)}\right\rangle \propto \prod_{x \leqslant 0} d^\dagger_x b^\dagger_x \left\vert{0}\right\rangle. \nonumber \end{align} \tag{ 44 }$$

Notice that this is nothing but a product state, of the form $\left\vert{{\rm DWIS}_b }\right\rangle \otimes \left\vert{{\rm DWIS}_d }\right\rangle =$ $\left(\prod_{x \leqslant 0} d^\dagger_x \left\vert{0_b}\right\rangle \right) \otimes \left(\prod_{x \leqslant 0} b^\dagger_x \left\vert{0_b}\right\rangle \right)$, where $\left\vert{0_{b, d}}\right\rangle$ is the vacuum for the b- and d-modes respectively. The two independent initial states for the b's and the d's are domain-wall initial states (DWIS), in the terminology of [80, 82]. At this point, the b- and d-modes are completely decoupled, and we are left with two independent time-evolutions,

$$\begin{align} \left\vert{{\rm JW}_0(t)}\right\rangle = {\rm e}^{-{\rm i} H_b t} \left\vert{{\rm DWIS}_b }\right\rangle \otimes {\rm e}^{-{\rm i} H_d t} \left\vert{{\rm DWIS}_d }\right\rangle \nonumber \end{align} \tag{ 45 }$$

where $H_b = \int \frac{{\rm d}k}{2\pi} \varepsilon (k) b^\dagger(k) b(k)$ and $H_d = -\int \frac{{\rm d}k}{2\pi} \varepsilon (k) d^\dagger(k) d(k)$. For the XX chain, the dispersion relation is $\varepsilon(k) = -\cos k$.

The entanglement entropy of the XX chain evolving from a DWIS was calculated recently numerically in [83], and from CFT in [81], relying on previous results [80]. The result is the following. For a bipartition $A\cup B$ with $A=(-\infty, x]$ and $B=(x, +\infty)$, the entanglement entropy of ${\rm e}^{- {\rm i} H t} \left\vert{{\rm DWIS}}\right\rangle$ remains constant while $v t < \vert x\vert$, and then behaves as $S_\alpha ({\rm e}^{- {\rm i} H t} \left\vert{{\rm DWIS}}\right\rangle) = \frac{1}{12} (1+\frac{1}{\alpha}) \log \left[ t (1-(x/vt){\hspace{0pt}}^2){\hspace{0pt}}^{3/2} \right]$ when $vt > \vert x\vert$. We can import this result and apply it directly to the present situation. Since, here, we have two copies of the XX chain (one for the b's and one for the d's), the OSEE of the JW string is exactly twice the entanglement entropy of the latter,

$$\begin{align} S_\alpha ({\rm JW}_0(t)) = \frac{1}{6} \left(1+\frac{1}{\alpha}\right) \log \left[ t \left(1-(x/t)^2\right)^{\frac{3}{2}} \right]. \nonumber \end{align} \tag{ 46 }$$

We note that the recent paper [84] studied a closely related problem about the growth of the (usual) entanglement entropy in the Ising chain.

5.2.1. Breakdown of conformal invariance.

As we just emphasized, the free fermion case strongly suggests that it is not possible to get the generic behavior of the OSEE of operators in Heisenberg picture as easily as, say, the generic linear growth of the OSEE of the evolution operator (section 4). The traditional CFT framework [68], which we used in sections 2–4, apparently breaks down. In this last section, we briefly discuss what happens if we try to attack this problem within this standard setup, and argue that it cannot be correct. Indeed, simple arguments lead to the conclusion that the OSEE should either be bounded or increase linearly in time; both scenarios are clearly ruled out by the numerical results in figure 9.

5.3.1. OSEE of operator in Heisenberg picture as an out-of-time-ordered correlator.

We assume once again that we are dealing with an infinitely long one-dimensional quantum critical system, described by a CFT with gapless excitations that propagate at velocity v. We consider a local operator ϕ that acts at point x  =  0, and we want to calculate the OSEE of $\phi(t) = {\rm e}^{{\rm i} H t} \phi {\rm e}^{-{\rm i} Ht}$ for the bipartition $A \cup B$, $A = (-\infty, x)$, $B= [ x, +\infty)$. We use equation (5) for integer α,

$$\begin{align*} S_\alpha (\phi(t)) = \frac{1}{1-\alpha} \log \left(\frac{{\rm tr} [ \phi^\dagger(t)^{\otimes \alpha} \cdot \mathcal{S}^A_{(1,\alpha,\alpha-1,\dots,2)} \cdot \phi (t)^{\otimes \alpha} \cdot \mathcal{S}^A_{(1,2,3,\dots,\alpha)} ] }{{\rm tr} [ \phi^\dagger(t)^{\otimes \alpha} \cdot \phi (t)^{\otimes \alpha}] } \right). \nonumber \end{align*}$$

The numerator of the expression in the logarithm involves a correlator with operators ϕ and $\phi^\dagger$ acting at time t, while the two swap operators act at time t  =  0; we see, however, that they appear in the order $\Phi^\dagger(t) \mathcal{S}^A \Phi(t) \mathcal{S}^A$—with $\Phi^\dagger(t) \equiv \phi(t){\hspace{0pt}}^\alpha$—, as opposed to the time-ordered version which would look like $\Phi^\dagger(t) \Phi(t) \mathcal{S}^A \mathcal{S}^A$. This type of correlator, dubbed out-of-time-ordered correlator or OTOC, has attracted a lot of attention recently in relation with quantum chaos and black holes and the AdS/CFT correspondence, following [85–88].

As in sections 2–4, one can rewrite the correlator appearing in equation (5) in terms of twist operators, replacing the swap $\mathcal{S}^A$ by a pair of twist operators, $\mathcal{S}^A \rightarrow \mathcal{T}_{(1, 2, \dots, \alpha)} \mathcal{T}_{(1, \alpha, \dots, 2)}$, with the two $\mathcal{T}$'s located at the two ends of the interval $A = (-\infty, x)$. It is convenient to think of the normalization of the twist operators as not being fixed at this point, so that $\mathcal{S}^A \propto \mathcal{T}(-\infty) \mathcal{T}(x)$, with some proportionality constant left undetermined. The two twist operators at $-\infty$, associated with $\mathcal{S}^A_{(1, 2, \dots, \alpha)}$ and $\mathcal{S}^A_{(1, \alpha, \dots, 2)}$ respectively, can be fused together to give the identity. So, in CFT, what we need to calculate is a ratio of the form

$$\begin{align} \label{eq:OTOC_twists} \frac{{\rm tr} [ \phi^\dagger(t)^{\otimes \alpha} \cdot \mathcal{T}_{(1,\alpha,\alpha-1,\dots,2)} \cdot \phi (t)^{\otimes \alpha} \cdot \mathcal{T}_{(1,2,3,\dots,\alpha)} ] }{{\rm tr} [ \phi^\dagger(t)^{\otimes \alpha} \cdot \phi (t)^{\otimes \alpha}] \; {\rm tr} [ \mathcal{T}_{(1,\alpha,\alpha-1,\dots,2)} \cdot \mathcal{T}_{(1,2,3,\dots,\alpha)} ]} \nonumber \end{align} \tag{ 47 }$$

where the two twist operators act at position x and at time 0. The new correlator of $\mathcal{T} \mathcal{T}$ in the denominator has been included in order to have a globally well-normalized expression, in the following sense. There are various operators that are located at the same position in equation (47), so there are several divergences that must be regulated in some way. But, whatever the regularization scheme, the ratio should always be one if ϕ is the identity operator; this is why the $\mathcal{T} \mathcal{T}$ factor in the denominator must be there.

5.3.2. Tentative regularization scheme for the out-of-time-ordered correlator, and expected outcome.

In order to make sense of the expression (47) in CFT, we use a regularization scheme that is similar to the one in [89, 90]. Once again, we switch to imaginary time; the order of insertion of the operators is taken into account by ensuring that they are inserted from right to left in increasing order of imaginary time. Thus, introducing some small parameter $\beta >0$, we see that one way of getting the numerator of equation (47) is to first evaluate the expression

$$\begin{align*} {\rm tr} [ ({\rm e}^{-(\frac{\beta}{4}-\frac{y}{v}) H } \phi^\dagger {\rm e}^{- (\frac{\beta}{4}+\frac{y}{v}) H })^{\otimes \alpha} \cdot \mathcal{T}_{(1,\alpha,\alpha-1,\dots,2)} \cdot ({\rm e}^{-(\frac{\beta}{4}-\frac{y}{v}) H } \phi {\rm e}^{-(\frac{\beta}{4}+\frac{y}{v}) H })^{\otimes \alpha} \cdot \mathcal{T}_{(1,2,3,\dots,\alpha)} ] \nonumber \end{align*}$$

with $-\frac{\beta v}{4} <y < \frac{\beta v}{4}$, then do the analytic continuation $y \rightarrow {\rm i} v t$, and then finally take β to zero. This expression (and similar ones for the denominator) is identical to an expectation value at finite inverse temperature β, involving operators at different imaginary times. So, for notational simplicity, we just write $\left\langle. \right\rangle$ for such expectation values, and use the complex coordinate $z = x+{\rm i} y$ to indicate the position of the different operators. Then the imaginary-time, regularized, version of the ratio (47) reads

$$\begin{align} \label{eq:OTOC_im} \frac{\left\langle \phi^{\dagger \otimes \alpha}(z_1, \overline{z}_1) \mathcal{T}_{(1,\alpha,\alpha-1,\dots,2)} (z_2, \overline{z}_2) \phi^{\otimes \alpha} (z_3, \overline{z}_3) \mathcal{T}_{(1,2,3,\dots,\alpha)}(z_4, \overline{z}_4) \right\rangle }{\left\langle \phi^{\dagger \otimes \alpha}(z_1, \overline{z}_1) \phi^{\otimes \alpha} (z_3, \overline{z}_3) \right\rangle \left\langle \mathcal{T}_{(1,\alpha,\alpha-1,\dots,2)} (z_2, \overline{z}_2) \mathcal{T}_{(1,2,3,\dots,\alpha)}(z_4, \overline{z}_4) \right\rangle }, \nonumber \end{align} \tag{ 48 }$$

with

$$\begin{align} \label{eq:list_z} z_1 = {\rm i} y+ {\rm i} \frac{\beta v}{4}, \qquad z_2 = x, \qquad z_3 = {\rm i} y-{\rm i} \frac{\beta v}{4}, \qquad z_4 = x - {\rm i} \frac{\beta v}{2}. \nonumber \end{align} \tag{ 49 }$$

It is convenient to map the infinite cylinder of circumference $\beta v$, where the coordinates z_j live, onto the plane, with the coordinates

$$\begin{align} \label{eq:zeta_z} \zeta_j = {\rm e}^{\frac{2\pi}{\beta v} z_j }. \nonumber \end{align} \tag{ 50 }$$

Now, what we need to analyze is the behavior of the four-point function $\frac{\left\langle \Phi_1^\dagger \mathcal{T}_2^\dagger \Phi_3 \mathcal{T}_4 \right\rangle }{\left\langle \Phi_1^\dagger \Phi_3\right\rangle \left\langle \mathcal{T}_2^\dagger \mathcal{T}_4 \right\rangle }$ in the plane (we use the notations $\Phi = \phi^{\otimes \alpha}$ and $\Phi_i = \Phi(z_i)$ for better readability), in the regime that corresponds to $y \rightarrow {\rm i} v t$, $t \rightarrow \infty$. Carrying out this exercise explicitly is a little bit technical, so we defer it to the appendix. For our purposes, it is however sufficient to make a few elementary remarks, which already show that the present framework cannot lead to sensible results.

As any other four-point correlator in the plane, the one that is of interest to us now must be expressible as a sum of conformal blocks,

$$\begin{align} \label{eq:block_G} \frac{\left\langle \Phi_1^\dagger \mathcal{T}_2^\dagger \Phi_3 \mathcal{T}_4 \right\rangle }{\left\langle \Phi_1^\dagger \Phi_3\right\rangle \left\langle \mathcal{T}_2^\dagger \mathcal{T}_4 \right\rangle } = \sum_p \mathcal{G}_p (\eta) \overline{\mathcal{G}}_p (\overline{\eta}), \qquad \quad \eta = \frac{\zeta_{12} \zeta_{34}}{\zeta_{13} \zeta_{24} }, \; \overline{\eta} = \frac{\overline{\zeta}_{12} \overline{\zeta}_{34}}{\overline{\zeta}_{13} \overline{\zeta}_{24} }, \nonumber \end{align} \tag{ 51 }$$

where $\zeta_{ij} \equiv \zeta_i-\zeta_j$. The sum runs over primary fields of the replicated CFT with α replicas, indexed by p. When $t \rightarrow \infty$, we see that $\eta \sim {\rm e}^{\frac{2 \pi t}{\beta}}$. In reasonable CFTs, such as for instance rational CFTs, the blocks behave algebraically in η (cases where this algebraic behavior is violated correspond to so-called logarithmic CFTs—for an introduction, see the special volume of J. Phys. A [91]—, where the behavior is logarithmic instead; however such cases appear in non-unitary theories which cannot arise as the low-energy description of one-dimensional quantum systems evolving unitarily, and are therefore of no interest to us here). In particular, if h_p is the conformal dimension of the primary operator in the intermediate channel, then the block behaves as $\mathcal{G}_p(\eta) \sim \eta^{-h_p}$. Thus, there are only two possibilities:

• the leading non-zero contribution comes from the identity block with h_p  =  0, such that $\frac{\left\langle \Phi^\dagger_1 \mathcal{T}^\dagger_2 \Phi_3 \mathcal{T}_4 \right\rangle }{\left\langle \Phi_1^\dagger \Phi_3\right\rangle \left\langle \mathcal{T}^\dagger_2 \mathcal{T}_4 \right\rangle }$ goes to a non-zero constant when $\eta \rightarrow \infty$. Then the OSEE at large time is just proportional to the logarithm of that constant. In other words, the OSEE of $\phi(t)$ is bounded (at least for α integer and $\geqslant 2$)
• for some reason, the contribution from the identity block vanishes, and it is some other primary operator (with h_p  >  0) that gives the leading contribution $\frac{\left\langle \Phi^\dagger_1 \mathcal{T}^\dagger_2 \Phi_3 \mathcal{T}_4 \right\rangle }{\left\langle \Phi_1^\dagger \Phi_3\right\rangle \left\langle \mathcal{T}^\dagger_2 \mathcal{T}_4 \right\rangle } \sim \eta^{- h_p} \overline{\eta}^{- \overline{h}_p}$. Since $\eta \sim {\rm e}^{\frac{2\pi t}{\beta}}$, we see that the OSEE of $\phi(t)$ increases linearly in time.

There is no room for other possibilities, such as e.g. logarithmic or power-law growths of the OSEE, based on this analysis and on this regularization scheme. Hence, this appears to be in contradiction with our free fermion calculation above, and with the numerical results diplayed in figure 9.

To conclude this section: it is unclear whether or not it is possible to use CFT methods similar to the ones of sections 2–4 to tackle the OSEE of local operators in Heisenberg picture for generic interacting systems. Certainly, the free fermion case cannot be tackled by such methods, because of the appearance of a domain-wall initial state that breaks conformal invariance. But whether or not a similar breakdown of conformal invariance appears also in generic interacting systems is unclear. We leave this as another open question, which deserves a more thorough investigation; it is a problem that is interesting both for practical numerical purposes, and from a low-dimensional field theory perspective.

Let α be a positive integer. Tensoring α replicas of a CFT—say in the plane $\mathbb{C}$—with central charge c, one gets a new CFT with central charge $\alpha c$. It is well-known that one can consider twist operators (dating back⁶ to [94, 95]), such that, turning around the operator, one goes from replica j to replica j  +  1 (mod α): Let us call $\mathcal{T}(w, \overline{w})$ this twist operator. Any operator acting in a single copy $\phi_j(z) = 1 \otimes 1\otimes \phi_j (z) \otimes 1 \otimes\dots \otimes 1$ must have a monodromy around $\mathcal{T}(w, \overline{w})$,

$$\begin{align} \label{eq:twist} \phi_j(w + (z-w){\rm e}^{{\rm i} \theta}) \mathcal{T}(w,\overline{w}) \rightarrow \phi_{j+1}(z) \mathcal{T}(w,\overline{w}) \nonumber \end{align} \tag{ B.1 }$$

when θ goes from zero to $2\pi$. This may not be sufficient to characterize $\mathcal{T}(w, \overline{w})$ entirely; there may be many operators that do that in the $\alpha^{\rm th}$-replicated theory, possibly with arbitrary high dimensions. We thus impose that $\mathcal{T}$ is, among the many operators that satisfy (B.1), the one with the smallest possible scaling dimension, which we write as $\Delta_\alpha = h_\alpha + \overline{h}_\alpha$. Notice that $\mathcal{T}$ treats right-moving and left-moving modes on the same footing, and so it must have spin zero: $h_\alpha = \overline{h}_\alpha$.

We then proceed to show that $\mathcal{T}(w, \overline{w})$ is a primary operator. To see that, notice that (B.1) holds, in particular, for the stress-tensor on the $j^{\rm th}$ sheet, T_j(z), such that the total stress-tensor of the full theory, $T(z) = \sum_{j=1}^\alpha T_j (z)$, has no monodromy. It follows that its OPE with $\mathcal{T}(w)$ is of the form $T(z) \mathcal{T}(w) \simeq \sum_{p\in \mathbb{Z}} (z-w){\hspace{0pt}}^{-2-p} \mathcal{T}_p(w)$, for some operators $\mathcal{T}_p$ with conformal dimensions $\Delta_\alpha + p$. But, by definition, there cannot be such operators with a scaling dimension smaller than $\Delta_\alpha$. So all the terms with p  <  0 actually vanish. The OPE must then be of the form

$$\begin{align} T(z) \mathcal{T}(w) \simeq \frac{h_\alpha}{(z-w)^2} \mathcal{T}(w) + \frac{1}{z-w} \partial \mathcal{T}(w) + {\rm regular\,\, terms}, \nonumber \end{align} \tag{ B.2 }$$

as claimed. The conformal weight $h_\alpha$ is easily evaluated. One way of calculating it is to look at the full theory, with central charge $\alpha c$, on a long cylinder of circumference $\ell$, and of length $\Lambda \gg \ell$: the partition function behaves as

$$\begin{align} Z \sim {\rm e}^{\Lambda \frac{\pi \alpha c}{6 \ell}}. \nonumber \end{align} \tag{ B.3 }$$

Now insert the twist operator at both ends of the cylinder, at $\pm\! \Lambda/2$. This results in a new partition function

$$\begin{align} Z_{\rm twist} &= \left\langle \mathcal{T}(-\Lambda/2) \mathcal{T}(+\Lambda/2) \right\rangle \times Z \nonumber \\ & \sim {\rm e}^{- \Lambda \frac{2\pi}{\ell} \left(\Delta_\alpha -\frac{\alpha c}{12} \right) }. \nonumber \end{align} \tag{ B.4 }$$

On the other hand, as one turns around the cylinder α times, it is clear that one goes through every replica exactly once, and so Z_twist is nothing but the partition function of a CFT with central charge c, on a cylinder of circumference $\alpha \ell$. Therefore,

$$\begin{align} Z_{\rm twist} \sim {\rm e}^{\Lambda \frac{\pi c}{6 \alpha\ell}}, \nonumber \end{align} \tag{ B.5 }$$

and, as a consequence,

$$\begin{align} \Delta_\alpha = \frac{c}{12} \left(\alpha - \frac{1}{\alpha} \right). \nonumber \end{align} \tag{ B.6 }$$

This argument for fixing the dimension of twist operators seems to be very well-known, it appeared for instance⁷ in [96].

It is very natural to generalize this to arbitrary permutations in $\sigma \in S_\alpha$, namely to define the twist operator $\mathcal{T}_{\sigma}$ as the operator with smallest possible scaling dimension such that

$$\begin{align} \label{eq:sigma} \phi_j(w + (z-w){\rm e}^{{\rm i} \theta}) \mathcal{T}(w,\overline{w}) \rightarrow \phi_{\sigma(\,j)}(z) \mathcal{T}(w,\overline{w}) \nonumber \end{align} \tag{ B.7 }$$

when θ goes from zero to $2\pi$. It is clear, from the same argument as above, that this operator is primary. If σ has cycles of lengths $\alpha_1 +\dots +\alpha_c = \alpha$, then the dimension of $\mathcal{T}_\sigma$ is obviously

$$\begin{align} \nonumber\Delta_\sigma &= \Delta_{\alpha_1} + \dots + \Delta_{\alpha_c} \nonumber \\ &= \frac{c}{12} \left(\alpha - \frac{1}{\alpha_1} - \dots -\frac{1}{\alpha_c} \right). \nonumber \end{align} \tag{ B.8 }$$

Notice that $\Delta_\sigma$ is maximal for circular permutations. The OPE of these arbitrary twist operators must mimic the group law of $S_\alpha$, which is a non-abelian group. So, to define the fusion of two twist operators, one need to a little bit careful. One way to do this is to fix some curve $\Gamma \subset \mathbb{C}$, and to assume that the twist operators are all lying on this curve; like beads on a string. Then the twist operators are naturally ordered: they appear in a given order as one moves along $\Gamma$. This allows to write the non-abelian OPE satisfied by the twist operators,

$$\begin{align} \mathcal{T}_{\sigma_1}(z) \mathcal{T}_{\sigma_2}(w) \simeq \; \left\vert z-w\right\vert ^{-\Delta_{\sigma_1} - \Delta_{\sigma_2} + \Delta_{\sigma_1 \circ \sigma_2} } \mathcal{T}_{\sigma_1 \circ \sigma_2}(w) + \dots \nonumber \end{align} \tag{ B.9 }$$

where z and w lie on $\Gamma$. In general, calculating correlation functions of more than three twist operators is a difficult task; early attempts include [96].

For simplicity, let us assume that ϕ is a primary operator with conformal spin zero, and that $\phi = \phi^\dagger$. Our starting point is the four-point function of ϕ in the plane $\mathbb{C}$, expressed as a sum over conformal blocks

$$\begin{align*} \left\langle \phi_1 \phi_2 \phi_3 \phi_4 \right\rangle _{\mathbb{C}} = \frac{1}{\vert w_{13} \vert ^{2 \Delta} \vert w_{24} \vert ^{2 \Delta}} \; \sum_p \mathcal{F}_p (\omega ) \overline{\mathcal{F}}_p \left(\overline{\omega} \right), \qquad \quad \omega = \frac{w_{12} w_{34}}{w_{13} w_{24} }. \nonumber \end{align*}$$

Here we are thinking of a diagonal theory, such as one of the minimal models of the A_n series. As usual, there is some arbitrariness in the choice of the form of the prefactor $1/ (\vert z_{13} \vert ^{2 \Delta} \vert z_{24} \vert ^{2 \Delta})$ and of the blocks themselves; here we have made our choice such that the conformal blocks transform conveniently under crossing, or 'F-move',

$$\begin{align*} \mathcal{F}_p (1-\omega) = \sum_q F_{p q} \mathcal{F}_q (\omega) \nonumber \end{align*}$$

with a matrix F that has constant entries (they do not depend on ω). Pictorially, this is of course the famous relation

As $\omega \rightarrow 0$, the conformal blocks behave as $\mathcal{F}_p (\omega) \sim \frac{c_p}{\omega^{2 h_\phi - h_p} }$ where $h_p \geqslant 0$ is the conformal dimension of the operator in the intermediate channel, and c_p is a constant equal to a product of OPE coefficients. If we label the identity block by p  =  1, then c₁  =  1 in the standard convention for normalization of the OPE coefficient $C^1_{\phi \phi} =1$ (or, equivalently, the standard normalization of the two-point function $\left\langle \phi \phi \right\rangle$). Then we have

$$\begin{align} \label{eq:limit0} \omega^{\Delta_\phi} \mathcal{F}_p (\omega) \quad \mathop{\longrightarrow}\limits_{\omega \rightarrow 0} \quad \left\{\begin{array}{@{}ll@{}} 1 & {\rm if} \; p=1 \nonumber \\ 0 & {\rm otherwise}.\end{array} \right. \nonumber \end{align} \tag{ C.2 }$$

Next, we use the conformal map $\zeta \mapsto w(\zeta) = \sqrt{\zeta}$ to write the two-point correlation function of the doubled operator $\Phi \equiv \phi \otimes \phi$ living on the two-sheeted plane with twist operators $\mathcal{T} \equiv \mathcal{T}_{(1, 2)}$ inserted at 0 and $\infty$. With $w_1 = \sqrt{\zeta_1}$, $w_2 = \sqrt{\zeta_3}$, $w_3 = -\sqrt{\zeta_1}$, $w_4 = -\sqrt{\zeta_3}$, we get

$$\begin{align*} \nonumber \left\langle \mathcal{T} (\infty) \Phi (\zeta_1) \Phi (\zeta_3) \mathcal{T} (0) \right\rangle = \left\vert \frac{1}{16 \zeta_1 \zeta_3} \right\vert ^{2\Delta_\phi} \sum_p \mathcal{F}_p \left(\frac{2 - \frac{\sqrt{\zeta_1}}{\sqrt{\zeta_3}} - \frac{\sqrt{\zeta_3}}{\sqrt{\zeta_1}} }{4} \right) \overline{\mathcal{F}}_p \left(\frac{2 - \frac{\sqrt{\bar{\zeta}_1}}{\sqrt{\bar{\zeta}_3}} - \frac{\sqrt{\bar{\zeta}_3}}{\sqrt{\bar{\zeta}_1}} }{4} \right). \nonumber \end{align*}$$

With a special conformal transformations $\zeta \mapsto \frac{a \zeta+b}{c \zeta+d}$, we can get the same correlator with twists inserted at arbitrary positions $\zeta_2$ and $\zeta_4$. The result can be put in the form

$$\begin{align*} \left\langle \Phi_1 \mathcal{T}_2 \Phi_3 \mathcal{T}_4 \right\rangle = \frac{1}{\left\vert \zeta_{13} \right\vert ^{4 \Delta_\phi} \left\vert \zeta_{24} \right\vert ^{2 \Delta_2} } \sum_p \mathcal{G}_p (\eta ) \overline{\mathcal{G}}_p \left(\overline{\eta}\right), \qquad \eta = \frac{\zeta_{12} \zeta_{34}}{\zeta_{13} \zeta_{24}}, \nonumber \end{align*}$$

where $\Delta_2 = \frac{c}{8}$ is the scaling dimension of the twist operator $\mathcal{T}=\mathcal{T}_{(1, 2)}$, and where we define the new conformal blocks in terms of the $\mathcal{F}$'s,

$$\begin{align} \label{eq:defG} \mathcal{G}_p (\eta) \equiv \left(\frac{1}{16 \eta(1-\eta)} \right)^{\Delta_\phi} \mathcal{F}_p \left(\frac{2-\sqrt{\frac{\eta}{\eta-1}} - \sqrt{\frac{\eta-1}{\eta}} }{4} \right). \nonumber \end{align} \tag{ C.3 }$$

Notice that, as a consequence of equation (C.2), we have

$$\begin{align} \label{eq:limit1} \mathcal{G}_p (\eta) \quad \mathop{\longrightarrow}\limits_{\eta \rightarrow \infty} \quad \left\{\begin{array}{@{}ll@{}} 1 & {\rm if} \; p=1 \nonumber \\ 0 & {\rm otherwise}.\end{array} \right. \nonumber \end{align} \tag{ C.4 }$$

The quotient by the two-point functions is

$$\begin{align*} \frac{\left\langle \Phi_1 \mathcal{T}_2 \Phi_3 \mathcal{T}_4 \right\rangle }{\left\langle \Phi_1 \Phi_3 \right\rangle \left\langle \mathcal{T}_2 \mathcal{T}_4 \right\rangle } = \sum_p \mathcal{G}_p (\eta ) \overline{\mathcal{G}}_p \left(\overline{\eta}\right), \nonumber \end{align*}$$

which is of the form (51) discussed in the main text. $\mathcal{G}_p$ has a branch cut between 0 and 1: indeed, as η turns around 1, intersecting the interval (0, 1) once, we see that we pick the other root in the argument of $\mathcal{F}$ in (C.3). So, if we let $\mathcal{G}^{\prime}_p (\eta)$ be the analytic continuation of $\mathcal{G}_p (\eta)$ when η crosses the branch cut, then

$$\begin{align*} \nonumber\mathcal{G}^{\prime}_p (\eta) &= \left(\frac{1}{16 \eta (1-\eta)} \right)^{\Delta_\phi} \mathcal{F}_p \left(\frac{2+\sqrt{\frac{\eta}{\eta-1}} + \sqrt{\frac{\eta-1}{\eta}} }{4} \right) \nonumber \\ \nonumber &= \left(\frac{1}{16 \eta (1-\eta)} \right)^{\Delta_\phi} \sum_q F_{p q} \mathcal{F}_q \left(\frac{2-\sqrt{\frac{\eta}{\eta-1}} - \sqrt{\frac{\eta-1}{\eta}} }{4} \right) \nonumber \\ &= \sum_q F_{pq} \mathcal{G}_q (\eta). \nonumber \end{align*}$$

Finally, we can come back to our problem: we set $\zeta_1 = {\rm e}^{\frac{2\pi}{\beta v} ({\rm i} y + {\rm i} \frac{\beta v}{4}) }$, $\zeta_2 = {\rm e}^{\frac{2\pi}{\beta v} x}$, $\zeta_3 = {\rm e}^{\frac{2\pi}{\beta v} ({\rm i} y- {\rm i}\frac{\beta v}{4}) }$ and $\zeta_4 = {\rm e}^{\frac{2\pi}{\beta v} (x - {\rm i} \frac{\beta v}{2})}$, see equations (49)–(50); then the cross-ratios become

$$\begin{align*} \eta = \frac{1}{2} \left({\rm e}^{{\rm i} \frac{\pi}{4}} {\rm e}^{\frac{\pi}{\beta v} (x-{\rm i} y) } + {\rm e}^{-{\rm i} \frac{\pi}{4}} {\rm e}^{-\frac{\pi}{\beta v} (x-{\rm i} y) } \right)^2 \qquad \overline{\eta} = \frac{1}{2} \left({\rm e}^{-{\rm i} \frac{\pi}{4}} {\rm e}^{\frac{\pi}{\beta v} (x+{\rm i} y) } + {\rm e}^{{\rm i} \frac{\pi}{4}} {\rm e}^{-\frac{\pi}{\beta v} (x+{\rm i} y) } \right)^2. \nonumber \end{align*}$$

Taking the analytic continuation $y \rightarrow {\rm i} v t$, this gives

$$\begin{align*} \eta = \frac{1}{2} \left({\rm e}^{{\rm i} \frac{\pi}{4}} {\rm e}^{\frac{\pi}{\beta v} (x+ vt) } + {\rm e}^{-{\rm i} \frac{\pi}{4}} {\rm e}^{-\frac{\pi}{\beta v} (x+ v t) } \right)^2 \qquad \overline{\eta} = \frac{1}{2} \left({\rm e}^{-{\rm i} \frac{\pi}{4}} {\rm e}^{\frac{\pi}{\beta v} (x- v t) } + {\rm e}^{{\rm i} \frac{\pi}{4}} {\rm e}^{-\frac{\pi}{\beta v} (x- v t) } \right)^2. \nonumber \end{align*}$$

Notice that $\overline{\eta}$ is not the complex conjugate of η after this operation. From now on, we focus on the case x  <  0, and we assume $\beta v \ll \vert x\vert$. Then as t goes from zero to infinity, we see that $\overline{\eta}$ always stays close to $\infty$. On the contrary, η starts from $\infty$ at t  =  0, then moves close to 1, turns around 1, crossing the segment [0,1] once, and then goes back to $\infty$ as $t\rightarrow +\infty$. In other words, while $\overline{\mathcal{G}}_p (\overline{\eta})$ quickly converges to $\overline{\mathcal{G}}_p (\infty)$, the behavior of $\mathcal{G}_p (\eta)$ is slightly more subtle, because the branch cut between 0 and 1 is crossed once when t goes from 0 to $+\infty$. In consequence, we have $\mathcal{G}_p (\eta) \rightarrow \mathcal{G}^{\prime}_p (\infty) = \sum_q F_{pq}\mathcal{G}_q (\infty)$ when $t \rightarrow +\infty$. Putting everything together, we find

$$\begin{align} \label{eq:OTOC_F11} \frac{\left\langle \Phi_1 \mathcal{T}_2 \Phi_3 \mathcal{T}_4 \right\rangle }{\left\langle \Phi_1 \Phi_3 \right\rangle \left\langle \mathcal{T}_2 \mathcal{T}_4 \right\rangle } \mathop{\longrightarrow}\limits_{t\rightarrow + \infty} \sum_{p,q} \overline{\mathcal{G}}_p (\infty) F_{p q} \mathcal{G}_q (\infty) = F_{11}, \nonumber \end{align} \tag{ C.5 }$$

where we used equation (C.4). Equation (C.5) is the main result of our attempt at deriving the OSEE of $\phi(t)$: we see that, according to the discussion in section 5.3.2, the OSEE should either increase linearly or remain bounded, depending whether F₁₁ is zero or non-zero. At least, this is the result that we obtain within the setup of section 5.3. As already pointed out in the main text, this is in contradiction both with the free fermion case, where an alternative analytical calculation is available (section 5.2) and leads to a logarithmic growth of the OSEE, and with the numerics of figure 9.

For concreteness, let us treat the example of the σ-field in the Ising model explicitly. σ is a spinless primary operator of scaling dimension $\Delta_\sigma=1/8$. In the plane $\mathbb{C}$, we have

$$\begin{align*} \left\langle \sigma_1 \sigma_2 \sigma_3 \sigma_4 \right\rangle = \frac{1}{\vert w_{13}\vert ^{\frac{1}{4}} } \frac{1}{\vert w_{24}\vert ^{\frac{1}{4}} } \frac{\left\vert 1+\sqrt{1-\omega} \right\vert + \left\vert 1-\sqrt{1-\omega} \right\vert }{2 \vert \omega (1-\omega) \vert ^{\frac{1}{4}}}, \qquad \quad \omega = \frac{w_{12} w_{34}}{w_{13} w_{24}}. \nonumber \end{align*}$$

The two conformal blocks are

$$\begin{align*} \mathcal{F }_1 (\omega) = \frac{\sqrt{1+ \sqrt{1-\omega} } }{\sqrt{2} \omega^{\frac{1}{8}} (1-\omega)^{\frac{1}{8}}} \qquad \mathcal{F}_2 (\omega) = \frac{\sqrt{1- \sqrt{1-\omega} } }{\sqrt{2} \omega^{\frac{1}{8}} (1-\omega)^{\frac{1}{8}}}, \nonumber \end{align*}$$

with the F-matrix

$$\begin{align*} F = \frac{1}{\sqrt{2}} \left(\begin{array}{@{}cc@{}} 1 & 1 \nonumber \\ 1 & -1 \end{array} \right). \nonumber \end{align*}$$

Notice that equation (C.2) is satisfied. Thus, we find that the OTOC is, at large time,

$$\begin{align*} \frac{\left\langle \sigma^{\otimes 2} \mathcal{T} \sigma^{\otimes 2} \mathcal{T} \right\rangle }{\left\langle \sigma^{\otimes 2} \sigma^{\otimes 2} \right\rangle \left\langle \mathcal{T} \mathcal{T} \right\rangle } \mathop{\longrightarrow}\limits_{t\rightarrow + \infty} \frac{1}{\sqrt{2}} \neq 0. \nonumber \end{align*}$$

This result predicts that the OSEE of $\sigma(t)$ should be bounded, which is in contradiction with the discussion of sections 5.1 and 5.2: after the Jordan–Wigner transformation, the σ-field in the transverse Ising chain is an operator that is attached to a JW string, so its OSEE is known to grow logarithmically.