Universal scrambling in gapless quantum spin chains

TL;DR

This paper investigates how scale invariance in gapless quantum spin chains influences information scrambling, showing that universal behavior of the out-of-time-ordered correlator (OTOC) is linked to the conformal field theory description with central charge c=1.

Contribution

It demonstrates that in the gapless phase of the XXZ model, the OTOC decays universally, independent of anisotropy, and connects this to conformal field theory predictions.

Findings

Faster OTOC decay in gapless vs gapped regions.

OTOC behavior is universal across anisotropy in the gapless phase.

Estimated central charge c is close to 1, confirming conformal invariance.

Abstract

Information scrambling, characterized by the out-of-time-ordered correlator (OTOC), has attracted much attention, as it sheds new light on chaotic dynamics in quantum many-body systems. The scale invariance, which appears near the quantum critical region in condensed matter physics, is considered to be important for the fast decay of the OTOC. In this paper, we focus on the one-dimensional spin-1/2 XXZ model, which exhibits quantum criticality in a certain parameter region, and investigate the relationship between scrambling and the scale invariance. We quantify scrambling by the averaged OTOC over the Pauli operator basis, which is related to the operator space entanglement entropy (OSEE). Using the infinite time-evolving block decimation (iTEBD) method, we numerically calculate time dependence of the OSEE in the early time region in the thermodynamic limit. We show that the averaged OTOC decays faster in the gapless region than in the gapped region. In the gapless region, the averaged OTOC behaves in the same manner regardless of the anisotropy parameter. This result is consistent with the fact that the low energy excitations of the gapless region belong to the same universality class as the Tomonaga-Luttinger liquid with the central charge c = 1. Furthermore, we estimate c by fitting the numerical data of the OSEE with an analytical result of the two-dimensional conformal field theory, and confirmed that c is close to unity. Thus, our numerical results suggest that the scale invariance is crucial for the universal behavior of the OTOC.

Tables

Table 1. Table 1: Fitting results for the OSEE shown in Fig.  5 with a fitting function S A OSEE ​ ( 2 ) = c 4 ​ ln ⁡ L eff + a superscript subscript 𝑆 𝐴 OSEE 2 𝑐 4 subscript 𝐿 eff 𝑎 {S_{A}^{\mathrm{OSEE\>(2)}}}=\frac{c}{4}\ln L_{\mathrm{eff}}+a .

Fitting range	$c$	$a$
$[2.5,4]$	0.957 $\pm$0.008	0.571 $\pm$ 0.007
$[3,4]$	0.946 $\pm$0.016	0.581 $\pm$ 0.014
$[3.5,4]$	0.96 $\pm$0.06	0.569 $\pm$ 0.05

Table 2. Table 2: Fitting results for the OSEE with a fitting function S A OSEE ​ ( 2 ) = c 4 ​ ln ⁡ L eff + a superscript subscript 𝑆 𝐴 OSEE 2 𝑐 4 subscript 𝐿 eff 𝑎 {S_{A}^{\mathrm{OSEE\>(2)}}}=\frac{c}{4}\ln L_{\mathrm{eff}}+a with χ = 256 𝜒 256 \chi=256 .

Fitting range	$c$	$a$
$[2.5,4]$	0.939 $\pm$0.016	0.584 $\pm$ 0.013
$[3,4]$	0.90 $\pm$0.04	0.62 $\pm$ 0.03
$[3.5,4]$	0.89 $\pm$0.10	0.62 $\pm$ 0.10

Table 3. Table 3: Fitting results for the OSEE with a fitting function S A OSEE ​ ( 2 ) = c 4 ​ ln ⁡ L eff + a superscript subscript 𝑆 𝐴 OSEE 2 𝑐 4 subscript 𝐿 eff 𝑎 {S_{A}^{\mathrm{OSEE\>(2)}}}=\frac{c}{4}\ln L_{\mathrm{eff}}+a and χ = 384 𝜒 384 \chi=384 .

FItting range	$c$	$a$
$[2.5,4]$	0.969 $\pm$0.014	0.5663 $\pm$ 0.011
$[3,4]$	0.94 $\pm$0.03	0.59 $\pm$ 0.02
$[3.5,4]$	0.90 $\pm$0.08	0.63 $\pm$ 0.07

Equations

$$F(t):=⟨\hat{W}(t)\hat{V}(0)\hat{W}(t)\hat{V}(0){⟩}_{\beta },$$

$$⟨\hat{W}(t)\hat{V}\hat{W}(t)\hat{V}{⟩}_{\beta }\to \text{tr}\text{quantity}[\hat{y}\hat{W}(t)\hat{y}\hat{V}\hat{y}\hat{W}(t)\hat{y}\hat{V}],$$

$${\overline{F}}_{AB}(t)$$

$${\overline{F}}_{\mathrm{AB}}(t)\le {\overline{F}}_{{\mathrm{A}}^{\prime }{\mathrm{B}}^{\prime }}(t)$$

$$|\hat{W}⟩:=\hat{W}\otimes \hat{I}|\Phi ⟩,$$

$$S_{\mathrm{A}}^{\mathrm{OSEE}}\text{quantity}[\hat{W}]$$

$$\hat{\rho }$$

$$S_{\mathrm{A}}^{\mathrm{OSEE}(n)}\text{quantity}[\hat{W}]:=\frac{1}{1-n}\log (\frac{\text{tr}\left[{\left({\hat{W}}^{†}\right)}^{\otimes n}{\hat{S}}_{\mathrm{A}}^{(1,n,n-1,\cdots ,2)}{\hat{W}}^{\otimes n}{\hat{S}}_{\mathrm{A}}^{(1,2,3,\cdots ,n)}\right]}{\text{tr}\left[{\hat{W}}^{†}\hat{W}\right]}),$$

$${\overline{F}}_{\mathrm{AB}}(t)$$

$$\frac{Z(\beta /2{)}^2}{2^{L}Z(\beta )}{\mathrm{e}}^{-S_{\mathrm{A}}^{\mathrm{OSEE}(2)}}\le {\overline{F}}_{{\mathrm{A}}^{\prime }{\mathrm{B}}^{\prime }}(t)$$

$$S_{\mathrm{A}}^{\mathrm{OSEE}(\mathrm{n})}\text{quantity}[\hat{U}(\beta ,t)]=\frac{c}{6}\text{quantity}(1+\frac{1}{n})\ln (\cosh (\frac{2\pi t}{\beta }))+S_0,$$

$${\overline{F}}_{\mathrm{AB}}(t)\propto \exp \text{quantity}(-S_{\mathrm{A}}^{\mathrm{OSEE}(2)}\text{quantity}[\hat{U}(\beta ,t)])\propto \frac{1}{\text{quantity}(\cosh (\frac{2\pi t}{\beta }){)}^{\frac{c}{4}}}.$$

$$\frac{⟨\hat{W}(t)\hat{V}\hat{W}(t)\hat{V}{⟩}_{\beta }}{⟨\hat{W}\hat{W}{⟩}_{\beta }⟨\hat{V}\hat{V}{⟩}_{\beta }}\sim {\mathrm{e}}^{-4\pi h_{v}t},$$

$${\hat{H}}_{\mathrm{XXZ}}:=J\sum _{i=-\infty }^{\infty }\text{quantity}[{\hat{\sigma }}_i^{\mathrm{X}}{\hat{\sigma }}_{i+1}^{\mathrm{X}}+{\hat{\sigma }}_i^{\mathrm{Y}}{\hat{\sigma }}_{i+1}^{\mathrm{Y}}+\Delta \text{quantity}({\hat{\sigma }}_i^{\mathrm{Z}}{\hat{\sigma }}_{i+1}^{\mathrm{Z}}-1)],$$

$${\mathrm{e}}^{h(\hat{A}+\hat{B})}={\mathrm{e}}^{h\hat{A}/2}{\mathrm{e}}^{h\hat{B}}{\mathrm{e}}^{h\hat{A}/2}+\text{order}{h}^3,$$

$$|{\mathrm{e}}^{-\beta \hat{H}}⟩\simeq \frac{1}{\sqrt{2}}\text{quantity}(|{\psi }_{\mathrm{GS}}^0⟩\otimes |{\psi }_{\mathrm{GS}}^0⟩+|{\psi }_{\mathrm{GS}}^1⟩\otimes |{\psi }_{\mathrm{GS}}^1⟩).$$

$$\begin{array}{rl}& S_A^{\mathrm{OSEE}(n)}\text{quantity}[{\mathrm{e}}^{-\frac{\beta }{4}\hat{H}}]\\ & \simeq \frac{1}{2}{S}_A^{\mathrm{OSEE}(n)}\text{quantity}[|{\psi }_{\mathrm{GS}}^0⟩⟨{\psi }_{\mathrm{GS}}^0|]+\frac{1}{2}{S}_A^{\mathrm{OSEE}(n)}\text{quantity}[|{\psi }_{\mathrm{GS}}^1⟩⟨{\psi }_{\mathrm{GS}}^1|]+\ln 2\\ & =S_A^{\mathrm{OSEE}(n)}\text{quantity}[|{\psi }_{\mathrm{GS}}^0⟩⟨{\psi }_{\mathrm{GS}}^0|]+\ln 2,\end{array}$$

$$⟨\mathcal{T}(u)\overline{\mathcal{T}}(v){⟩}_{\mathrm{cyl}}\propto L_{\mathrm{eff}}^{-{\Delta }_n},$$

$$\text{innerproduct}{\hat{O}}_{\mathrm{A}}{\hat{O}}_{\mathrm{B}}=\text{tr}\left[{\hat{O}}_{\mathrm{A}}^{†}{\hat{O}}_{\mathrm{B}}\right].$$

$$|{\hat{O}}_{\mathrm{A}}{\hat{O}}_{\mathrm{B}}⟩={\hat{O}}_{\mathrm{A}}\otimes \hat{I}|{\hat{O}}_{\mathrm{B}}⟩=\hat{I}\otimes {\hat{O}}_{\mathrm{B}}^T|{\hat{O}}_{\mathrm{A}}⟩.$$

$$\begin{array}{rl}F(t)& =\text{tr}[{\hat{O}}_B(\beta ,t){\hat{O}}_A{\hat{O}}_B(\beta ,t){\hat{O}}_A]\\ & =\text{innerproduct}{\hat{O}}_A{\hat{O}}_B(\beta ,t){\hat{O}}_B(\beta ,t){\hat{O}}_A\\ & =⟨{\hat{O}}_A|{\hat{O}}_B(\beta ,t)\otimes {\hat{O}}_B(\beta ,t{)}^{\ast }|{\hat{O}}_A⟩.\end{array}$$

$$\begin{array}{rl}F(t)& =\text{tr}[{\hat{O}}_A{\hat{O}}_B(\beta ,t){\hat{O}}_A{\hat{O}}_B(\beta ,t)]\\ & =\text{innerproduct}{\hat{O}}_B(\beta ,t){\hat{O}}_A{\hat{O}}_A{\hat{O}}_B(\beta ,t)\\ & =⟨{\hat{O}}_B(\beta ,t)|{\hat{O}}_A\otimes {\hat{O}}_A^{\ast }|{\hat{O}}_B(\beta ,t)⟩.\end{array}$$

$${\hat{O}}_{\mathrm{A}}={\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }}{\hat{O}}_{{\mathrm{A}}^{\prime }},$$

$$[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{{\mathrm{A}}^{\prime }}]=0.$$

$$\text{tr}[{\hat{O}}_{\mathrm{B}}(\beta ,t){\hat{O}}_{\mathrm{A}}{\hat{O}}_{\mathrm{B}}(\beta ,t){\hat{O}}_{\mathrm{A}}]$$

$$=\text{tr}[{\hat{O}}_{\mathrm{B}}(\beta ,t){\hat{O}}_{{\mathrm{A}}^{\prime }}{\hat{O}}_{\mathrm{B}}(\beta ,t){\hat{O}}_{{\mathrm{A}}^{\prime }}]$$

$$+\frac{1}{2}\text{tr}[{\hat{O}}_{{\mathrm{A}}^{\prime }}[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]{\hat{O}}_{{\mathrm{A}}^{\prime }}[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]].$$

$$\begin{array}{rl}& F_{\mathrm{AB}}(t)\\ =& \text{tr}[{\hat{O}}_{\mathrm{B}}(\beta ,t){\hat{O}}_{{\mathrm{A}}^{\prime }}{\hat{O}}_{\mathrm{B}}(\beta ,t){\hat{O}}_{{\mathrm{A}}^{\prime }}]\\ & +\frac{1}{2}\text{tr}[{\hat{O}}_{{\mathrm{A}}^{\prime }}[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]{\hat{O}}_{{\mathrm{A}}^{\prime }}[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]]\\ =& F_{{\mathrm{A}}^{\prime }\mathrm{B}}(t)-\frac{1}{2}\text{tr}[{\hat{O}}_{\mathrm{A}}^{\prime }[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t){]}^{†}{\hat{O}}_{\mathrm{A}}^{\prime }[O_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]]\\ =& F_{{\mathrm{A}}^{\prime }\mathrm{B}}(t)-\frac{1}{2}\text{innerproduct}[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]{\hat{O}}_{{\mathrm{A}}^{\prime }}{\hat{O}}_{{\mathrm{A}}^{\prime }}[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]\\ =& F_{{\mathrm{A}}^{\prime }\mathrm{B}}(t)-\frac{1}{2}⟨[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]|{\hat{O}}_{{\mathrm{A}}^{\prime }}\otimes {\hat{O}}_{{\mathrm{A}}^{\prime }}^{\ast }|[{\hat{O}}_{\mathrm{A}\setminus {\mathrm{A}}^{\prime }},{\hat{O}}_{\mathrm{B}}(\beta ,t)]⟩.\end{array}$$

$$\sum _{{\hat{O}}_{{\mathrm{A}}^{\prime }}\in {\mathcal{P}}_{{\mathrm{A}}^{\prime }}}{\hat{O}}_{{\mathrm{A}}^{\prime }}\otimes {\hat{O}}_{{\mathrm{A}}^{\prime }}^{\ast }=\frac{1}{2^{a^{\prime }}}|{\hat{I}}_{{\mathrm{A}}^{\prime }}⟩⟨{\hat{I}}_{{\mathrm{A}}^{\prime }}|,$$

$$\sum _{{\hat{O}}_{{\mathrm{A}}^{\prime }}\in {\mathcal{P}}_{{\mathrm{A}}^{\prime }}}{F}_{\mathrm{AB}}(t)\le \sum _{{\hat{O}}_{{\mathrm{A}}^{\prime }}\in {\mathcal{P}}_{{\mathrm{A}}^{\prime }}}{F}_{{\mathrm{A}}^{\prime }\mathrm{B}}(t).$$

$${\overline{F}}_{\mathrm{AB}}(t)\le {\overline{F}}_{{\mathrm{A}}^{\prime }\mathrm{B}}(t).$$

$$R_{\text{quantity}(\alpha {\alpha }^{\prime }),\text{quantity}(\beta {\beta }^{\prime })}$$

Full text

††preprint: APS/123-QED

I Introduction

While classical chaos is characterized by the Lyapunov exponent Gallavotti (1999), the definition of chaos in quantum systems is non-trivial, and many studies have been conducted to understand various aspects of quantum chaos Stöckmann (1999); Haake (2001). In addition to the conventional indicators of quantum chaos such as the level statistics Stöckmann (1999); Haake (2001) and the Loschmidt echo Gorin et al. (2006), scrambling recently attracts attention because it sheds new light on chaotic behaviors in quantum many-body systems Hayden and Preskill (2007); Sekino and Susskind (2008). Scrambling represents the delocalization process of locally encoded quantum information. This is characterized by the out-of-time-ordered correlator (OTOC), which is defined as

[TABLE]

where $\hat{V}$ and $\hat{W}$ are local Hermitian operators, and $\langle\cdots\rangle_{\beta}$ represents the canonical ensemble average at the inverse temperature $\beta$. We note that the ordering of the times in the OTOC is different from that for the usual correlation functions. As the OTOC is closely related to the square of the commutator $C(t):=-\langle\quantity[\hat{W}(t),\hat{V}]^{2}\rangle$, it quantifies how the operators $W(t)$ and $V(0)$ become non-commutable as the spatial support of $W(t)$ increases in time evolution. In quantum systems, the squared commutator is expected to increase exponentially in the early time regime, if the system scrambles fast enough, which leads to the exponential behavior of the OTOC of the form $F(t)\approx a-b\mathrm{e}^{\lambda_{L}t}$ Maldacena et al. (2016a); Kitaev (2015). Here, $\lambda_{\mathrm{L}}$ is referred to as the quantum Lyapunov exponent. In terms of the quantum information theory, the averaged OTOC over the operator basis is equivalent to tripartite mutual information Hosur et al. (2016) and represents quantum pseudorandomness Roberts and Yoshida (2017).

As an important feature of the OTOC, its decay rate has an upper bound that is referred to as the chaotic bound or the Maldacena-Shenker-Stanford (MSS) bound Maldacena et al. (2016a). This bound is achieved by the Sachdev-Ye-Kitaev (SYK) model Kitaev (2015); Sachdev and Ye (1993); Sachdev (2015); Maldacena and Stanford (2016); Fu et al. (2016); Cotler et al. (2017), which exhibits the partial scale invariance SL(2,$\mathbf{R}$). Moreover, the two-dimensional conformal field theory (2D CFT) Francesco (1999) in the large central charge limit also achieves the upper bound Roberts and Stanford (2015). We emphasize that these two models exhibits the scale invariance. We note that the role of scale invariance has also been studied with the Tomonaga-Luttinger liquid (TLL) Dóra and Moessner (2017) and the Bose-Hubbard model Shen et al. (2017); Bohrdt et al. (2017).

In this paper, we focus on the relationship between the scale invariance and scrambling with a simple quantum spin chain. Specifically, we numerically investigate the one-dimensional quantum spin-$1/2$ anti-ferromagnetic XXZ chain, which exhibits a quantum phase transition between the gapless phase and the gapped phase. The low-energy excitations in the gapless phase is characterized by the TLL and those in the gapped phase are described by the Ising-like excitations Giamarchi (2003). For this model, we calculate the averaged OTOC that is equivalent to the operator space entanglement entropy (OSEE) Zanardi (2001); Prosen and Pizorn (2007); Znidaric et al. (2008); Pizorn and Prosen (2009); Xu and Swingle .

We show that there is a qualitative difference between the behavior of OTOC in the gapless region and that of gapped region: the decay of OTOC in the gapless region is faster than that of the gapped one. This result supports the expectation that the scale invariance of CFT in a quantum system leads to fast scrambling Maldacena et al. (2016a); Kitaev (2015). Furthermore, our numerical study shows that the decaying behavior of OTOC in the gapless region is always similar and does not depend on the anisotropy parameter. We therefore suggest that it is universal in the gapless regime but only at low temperatures. Here we argue that such a universal behavior comes from the fact that the low-energy excitations of the gapless XXZ model are described by the TLL and belong to the same universality class of the 2D CFT with $c=1$ Korepin (1993); Giamarchi (2003). We thus expect that the universal decaying behavior of OTOC in the gapless regime should hold for various quantum systems including non-integrable ones if their low-energy spectra are described by the 2D CFT with $c=1$. In fact, the TLL or the 2D CFT with $c=1$ is almost ubiquitous in the low-energy excitations not only for integrable gapless quantum spin chains such as the spin-1/2 XXZ chain but also for non-integrable ones.

Our numerical method is the infinite time-evolving block decimation (iTEBD) method Vidal (2003, 2007); Orús and Vidal (2008), which enables us to directly address the thermodynamic limit. This is a benefit of the iTEBD compared to other methods such as numerical exact diagonalization. We estimate the central charge $c$ by fitting the numerical results with an analytical result of the 2D CFT Dubail (2017). As a consequence, we confirmed that $c$ is close to unity, which is consistent with the fact that the XXZ model in the gapless region belongs to the universality class of $c=1$.

This paper is organized as follows. In Sec. II, we explain the averaged OTOC, the OSEE, and the relationship between them. In Sec. III, we show our main numerical results of the average OTOC in the gapped and the gapless regions in the thermodynamic limit. In Sec. IV, we summarize this paper and make some remarks. In Appendix A, we prove an inequality concerning the averaged OTOC. In Appendix B, we introduce the canonical form of the infinite matrix product operator (iMPO) Verstraete and Cirac (2004); Zwolak (2004); Barthel (2017). In Appendix C, we discuss the details of numerical errors.

II averaged OTOC and OSEE

In this section, we define the averaged OTOC and discuss its relationship to the operator space entanglement entropy (OSEE).

We define the OTOC at finite temperature by Eq. (1). In the quantum field theory (QFT), the OTOC requires appropriate regularizations and the following regularization is commonly adopted Maldacena and Stanford (2016); Roberts and Stanford (2015):

[TABLE]

where $\hat{y}^{4}:=\mathrm{e}^{-\beta\hat{H}}/Z(\beta)$ and $Z(\beta):=\tr\quantity[\mathrm{e}^{-\beta\hat{H}}]$. We adopt this regularized OTOC, which provides a platform to compare our study and the previous studies of QFT.

For simplicity, we assume that the total system is on a lattice and each local site is described by a two-level system (i.e., a qubit). Let S be the set of all sites and $L:=|{\rm S}|$ be the number of them. The Hilbert space of the total system is $\mathcal{H}=\quantity(\mathbb{C}^{2})^{\otimes L}$.

We consider subsystems $\mathrm{A},\mathrm{B}\subset\mathrm{S}$. The number of sites in A and B are denoted by $a$ and $b$, respectively. These subsystems are defined as $\mathrm{A}:=\{i_{1},i_{2},\cdots,i_{a}\}$ and $\mathrm{B}:=\{j_{1},j_{2},\cdots,j_{b}\}$. An operator belonging to the Pauli operator basis of A (B) is respectively written as $\hat{O}_{\mathrm{A}}$ ($\hat{O}_{\mathrm{B}}$). The operator is explicitly written by a tensor product $\hat{O}_{\mathrm{A}}=\hat{o}_{i_{1}}\otimes\hat{o}_{i_{2}}\otimes\cdots\otimes\hat{o}_{i_{a}}$ and $\hat{O}_{\mathrm{B}}=\hat{o}_{j_{1}}\otimes\hat{o}_{j_{2}}\otimes\cdots\otimes\hat{o}_{j_{b}}$, where $\hat{o}_{i}\in\quantity{\hat{I}_{i},\hat{\sigma}_{i}^{\mathrm{X}},\hat{\sigma}_{i}^{\mathrm{Y}},\hat{\sigma}_{i}^{\mathrm{Z}}}$ is a Pauli matrix at site $i\in{\rm S}$. We denote the sets of the operator basis of A and B as $\mathcal{P}_{\mathrm{A}}$ and $\mathcal{P}_{\mathrm{B}}$, respectively.

We now define the averaged OTOC as follows Hosur et al. (2016); Roberts and Yoshida (2017):

[TABLE]

where we averaged the OTOC over the operator bases $\mathcal{P}_{\mathrm{A}}$ and $\mathcal{P}_{\mathrm{B}}$. For the special case of $\beta=0$, it has been shown that the decay of the averaged OTOC is related to the negativity of tripartite mutual information Hosur et al. (2016).

We introduce a basic property of the averaged OTOC. For A and B, we consider subsystems $\mathrm{A}^{\prime}\subseteq\mathrm{A}$ and $\mathrm{B}^{\prime}\subseteq\mathrm{B}$, respectively. Then, the averaged OTOC satisfies

[TABLE]

(see the Appendix A for the proof). This inequality holds even when A and B have an intersection.

In the same setup, we define the OSEE. For the Hilbert space of the total system $\mathcal{H}$, we consider a copy of $\mathcal{H}$ and denote it as $\mathcal{\tilde{H}}$. We regard a linear operator $\hat{W}$ acting on $\mathcal{H}$ as a state on the extended Hilbert space $\mathcal{H}\otimes\mathcal{\tilde{H}}$ as

[TABLE]

where $\ket{\Phi}\coloneqq\prod_{i=1}^{L}(\ket{0}_{i}\ket{0}_{\tilde{i}}+\ket{1}_{i}\ket{1}_{\tilde{i}})$ is the EPR pair between the original system and the copied system and $\hat{I}$ is the identity operator. Here, the index $i$ represents the original system and $\tilde{i}$ represents its copy. We note that this mapping is essentially the same as the channel-state duality Choi (1975). When the total system is split into A and B, i.e., $\mathrm{S}=\mathrm{A}\cup\mathrm{B},\mathrm{A}\cap\mathrm{B}=\emptyset$, the OSEE of the operator $\hat{W}$ in the region A is defined as

[TABLE]

where

[TABLE]

with $\tilde{\mathrm{A}}$ being the copy of $\mathrm{A}$ and $\tilde{\mathrm{B}}$ being the copy of $\mathrm{B}$. Furthermore, the Rényi-$n$ OSEE of $\hat{W}$ is defined as

[TABLE]

where $\hat{S}_{\mathrm{A}}^{(1,2,3,\cdots,n)}$ is a permutation operator acting on $\mathcal{H}^{\otimes n}$, which cyclically permutes the replicas of the subsystem A. Explicitly, the action of $\hat{S}_{\mathrm{A}}^{(1,2,3,\cdots,n)}$ is given by $\hat{S}_{\mathrm{A}}^{(1,2,3,\cdots,n)}\ket{i_{1}}_{\mathrm{A}}\otimes\ket{j_{1}}_{\mathrm{B}}\otimes\ket{i_{2}}_{\mathrm{A}}\otimes\ket{j_{2}}_{\mathrm{B}}\otimes\cdots\otimes\ket{i_{n}}_{\mathrm{A}}\otimes\ket{j_{n}}_{\mathrm{B}}=\ket{i_{2}}_{\mathrm{A}}\otimes\ket{j_{1}}_{\mathrm{B}}\otimes\ket{i_{3}}_{\mathrm{A}}\otimes\ket{j_{2}}_{\mathrm{B}}\cdots\otimes\ket{i_{n}}_{\mathrm{A}}\otimes\ket{j_{n-1}}\otimes\ket{i_{1}}_{\mathrm{A}}\otimes\ket{j_{n}}_{\mathrm{B}}$, where $\ket{i}_{\mathrm{A}}$ belongs to a basis of the subsystem A and $\ket{j}_{\mathrm{B}}$ belongs to a basis of the subsystem B. We note that $\hat{S}_{\mathrm{A}}^{(1,n,n-1,\cdots,2)}$ is the inverse operator of $\hat{S}_{\mathrm{A}}^{(1,2,3,\cdots,n)}$.

The averaged OTOC has the following relation to the OSEE Hosur et al. (2016). When $\mathrm{S}=\mathrm{A}\cup\mathrm{B}$ and $\mathrm{A}\cap\mathrm{B}=\emptyset$, the averaged OTOC is written as

[TABLE]

where $\hat{U}(\beta,t)\coloneqq\mathrm{e}^{-\left(\frac{\beta}{4}+\mathrm{i}t\right)\hat{H}}$. From Eq. (4),

[TABLE]

holds for $\mathrm{A}^{\prime}\subseteq\mathrm{A}$ and $\mathrm{B}^{\prime}\subseteq\mathrm{B}$. Thus, the averaged OTOC $\overline{F}_{\mathrm{A}^{\prime}\mathrm{B}^{\prime}}(t)$ is bounded from below with the Rényi-2 OSEE. Although the OTOC with local observables is investigated in many previous studies, we focus on Eq. (9) to consider a lower bound of the averaged OTOC.

The OSEE is also investigated in QFT. The OSEE of $\hat{U}(\beta,t)$ in quantum critical systems has been calculated analytically by the 2D CFT Dubail (2017). In QFT, the total system is not on a lattice and subsystems are represented by intervals. With the bipartition $\mathrm{A}=\left(-\infty,0\right]$ and $\mathrm{B}=\left(0,\infty\right)$, the OSEE of $\hat{U}(\beta,t)$ is given by

[TABLE]

where $S_{0}$ is an unimportant constant. This formula is derived with the replica trick as the case of the Calabrese-Cardy formula Calabrese and Cardy (2004). We make a remark on the case of another partition $\mathrm{A}=\left[0,L_{\mathrm{A}}\right]$ and $\mathrm{B}=\mathrm{S}\setminus\mathrm{A}$. In the limit of $L_{\mathrm{A}}\gg 1$, the OSEE ${S_{\mathrm{A}}^{\mathrm{OSEE}\>(n)}}\quantity[\hat{U}(\beta,t)]$ becomes double of Eq. (11) because there are two boundaries.

We obtain an expression of the averaged OTOC in terms of the OSEE by substituting Eq. (11) into Eq. (9):

[TABLE]

In the same manner, for the partition with $\mathrm{A}=\left[0,L_{\mathrm{A}}\right]$, the averaged OTOC is obtained as $\overline{F}_{\mathrm{AB}}(t)\propto\quantity(\cosh(\frac{2\pi t}{\beta}))^{-c/2}$.

At early times $t/\beta\ll 1$, the above expression of the averaged OTOC is different from $F(t)\approx a-b\mathrm{e}^{\lambda_{L}t}$ that is expected from the chaotic bound Maldacena et al. (2016a); Kitaev (2015). This discrepancy comes from the fact that we did not take the semi-classical limit (e.g., the large-$c$ limit).

On the other hand, at late times $t/\beta\gg 1$, Eq. (12) reduces to the exponential decay $\overline{F}_{AB}(t)\sim\mathrm{e}^{-\frac{c\pi t}{\beta}}$. We note that from the 2D CFT with the large-$c$ limit Roberts and Stanford (2015), the OTOC with the primary fields $\hat{V}$ and $\hat{W}$ also decays exponentially at $t/\beta\gg 1$:

[TABLE]

where $h_{v}$ is the conformal dimension of $V$.

The exponent in Eq. (13) is different from that of the averaged OTOC, because the conformal dimension of the twist operator for $n=2$ is $c/16$. This is not a contradiction and our results are consistent with Eq. (13). This is because Eq. (13) is obtained in the parameter region $1\ll h_{v}\ll h_{w}$ Roberts and Stanford (2015) with $h_{w}$ being the conformal dimension of $W$, which is not satisfied in our setup. Besides, the positions of the primary fields in the spacetime in Ref. Roberts and Stanford (2015) are different from those in our setup.

III Numerical results by iTEBD

In this section, we present our main numerical results. Specifically, by using the iTEBD, we calculate the OSEE of $\hat{U}(\beta,t)=\mathrm{e}^{-\left(\frac{\beta}{4}+\mathrm{i}t\right)\hat{H}}$ and discuss the qualitative difference between the gapless and gapped regions.

III.1 Setup

In this subsection, we describe our setup for the numerical calculation and the protocol of the iTEBD. The Hamiltonian of the one-dimensional quantum spin-$1/2$ anti-ferromagnetic XXZ chain is defined as

[TABLE]

where $\hat{\sigma}_{i}^{\mathrm{X}},\hat{\sigma}_{i}^{\mathrm{Y}}$ and $\hat{\sigma}_{i}^{\mathrm{Z}}$ are the Pauli operators at site $i$, $\Delta$ is the anisotropy parameter, and $J$ $(>0)$ is the strength of the exchange interaction. For simplicity, we set $J=1$. For $0\leq\Delta\leq 1$, the low energy excitations from the ground state are gapless and this model exhibits quantum criticality. From the 2D CFT, it is known that the gapless XXZ model belongs to the universality class with the conformal charge $c=1$ Giamarchi (2003); Korepin (1993). We divide the total system into two semi-infinite subsystems $\mathrm{A}=\quantity{-\infty,\cdots,-1,0}$ and $\mathrm{B}=\quantity{1,2,\cdots\infty}$ and calculate the Rényi-2 OSEE for this bipartition.

We explain the protocol to approximately calculate the imaginary-time evolution operator by the iTEBD. We perform the Trotter decomposition Trotter (1959) of $\mathrm{e}^{-\frac{\beta}{4}\hat{H}}$ and obtain the products of $\mathrm{e}^{-\Delta\beta\hat{H}}$, where $\Delta\beta:=\beta/(4N_{\mathrm{Trotter}})$ and $N_{\mathrm{Trotter}}$ is the number of decomposition. In this study, we adopt the Trotter decomposition of the second order

[TABLE]

where we set $h=-\Delta\beta$. We split $\hat{H}=\hat{A}+\hat{B}$ into $\hat{A}$ and $\hat{B}$. Then, $\hat{A}$ ($\hat{B}$) corresponds to the interactions defined in the odd (even) bonds, respectively.

Then, we approximate the products of the infinitesimal time evolution operators in the form of the iMPO by repeatedly using the singular value decomposition (SVD). In general, the MPO provides a representation of a many-body operator, which decomposes it into a one-dimensional sequence of locally defined tensors. With the bond dimension $\chi$, we control the entanglement in the operator space. By assuming the translational invariance, the MPO can be applied to an infinite chain and the extended MPO for that case is referred to as the iMPO. Figure 1(a) shows the approximation protocol of $\mathrm{e}^{-\frac{\beta}{4}\hat{H}}$ by the iTEBD and the iMPO. By construction of the iTEBD protocol with the Trotter decomposition, we finally obtain a translational invariant iMPO composed of two tensors $\quantity{A_{1},A_{2}}$ with the two-site unit cell.

After performing the iTEBD, we transform the iMPO into the canonical form Orús and Vidal (2008) to calculate the OSEE. As illustrated in Fig. 1(b), the canonical form for the two-site translational invariant case is composed of four tensors $\quantity{\lambda_{\mathrm{A}_{1}},\Gamma_{\mathrm{A}_{1}},\lambda_{\mathrm{A}_{2}},\Gamma_{\mathrm{A}_{2}}}$, where $\Gamma_{\mathrm{A}_{1}}$ and $\Gamma_{\mathrm{A}_{2}}$ are complex-valued tensors with four indices, and $\lambda_{\mathrm{A}_{1}}$ and $\lambda_{\mathrm{A}_{2}}$ are diagonal matrices with non-negative elements. The definition of the canonical form is explained in Appendix B. Using the canonical form, we can calculate the Rényi-2 OSEE for the bipartition $\mathrm{A}=\quantity{-\infty,\cdots,-1,0}$ and $\mathrm{B}=\quantity{1,2,\cdots\infty}$.

We can also approximate the complex-time evolution operator $\hat{U}(\beta,t)$ in the same manner. In the case of the complex-time evolution, we first approximate $\hat{U}(\beta,0)$ in the form of the iMPO. By taking it as the initial state, we execute the iTEBD for $\hat{U}(0,t)$.

We note that the MPO formalism itself does not work very efficient for the direct calculation of the OTOC Hemery  Pollmann  and Luitz . In our study, we calculated OSEE by using Eq. (12) which relates the OSEE and the averaged OTOC.

III.2 Imaginary-time evolution

We now show our numerical results. Figure 2 shows the imaginary-time dependence of the Rényi-2 OSEE ${S_{A}^{\mathrm{OSEE\>(2)}}}$ for $\hat{U}(\beta,0)=\mathrm{e}^{-\frac{\beta}{4}\hat{H}}$. The anisotropy parameters are taken as $\Delta=0.25,0.5,0.75,1,1.6,1.8,2,3$. We note that a similar calculation of imaginary time evolution has been performed for a finite chain Barthel (2017), while our results are for an infinite chain in the thermodynamic limit.

In the high-temperature region with small $\beta$, the OSEE increases monotonically as $\beta$ increases regardless of the existence of the gap. However, in the low-temperature region with large $\beta$, the difference between the gapless and gapped regions becomes visible. While the OSEE continues to increase monotonically in the gapless region, the OSEE in the gapped region with $1<\Delta$ shows non-monotonic behavior. This result implies that the OSEE captures the difference between the excitation structure near the ground state, which is accessible if the temperature is sufficiently low.

We discuss the reason why the OSEE in the gapped region exhibits the non-monotonic behavior. The XXZ model in the gapped region has the two-fold degeneracy in the ground states, which exhibits the $\mathrm{Z}_{2}$ symmetry. Let $E_{0}$ be the eigenenergy of the ground states and $\ket{\psi^{0}_{\mathrm{GS}}}$ and $\ket{\psi^{1}_{\mathrm{GS}}}$ be the corresponding eigenvectors. In the gapped region and at sufficiently low temperature $\beta\gg 1$, we obtain $\mathrm{e}^{-\beta\hat{H}}\simeq\mathrm{e}^{-\beta E_{0}}\quantity(\ket{\psi^{0}_{\mathrm{GS}}}\bra{\psi^{0}_{\mathrm{GS}}}+\ket{\psi^{1}_{\mathrm{GS}}}\bra{\psi^{1}_{\mathrm{GS}}})$. Thus, the corresponding state vector in the extended Hilbert space is approximated as

[TABLE]

In our numerical calculation, however, the iMPO obtained by the iTEBD is two-site translation invariant, while the operator $\mathrm{e}^{-\frac{\beta}{4}\hat{H}}$ is one-site translation invariant. This implies that the one-site translation invariance is not guaranteed and can be broken in our results. In such a case, the ground-state degeneracy is removed by the artificial symmetry breaking, resulting in an artificial non-monotonic behavior in numerical calculations Jiang and Balents .

Due to the removal of the degeneracy, $\ket{\mathrm{e}^{-\beta\hat{H}}}$ is described by only one of $\ket{\psi^{0}_{\mathrm{GS}}}$ or $\ket{\psi^{1}_{\mathrm{GS}}}$. This is the reason why the OSEE in the gapped region shows the non-monotonic $\beta$-dependence in Fig. 2. As the temperature decreases, the removal of the degeneracy becomes more visible, and the OSEE decays more rapidly. As a special case, if the anisotropy parameter is very large $\Delta\gg 1$, the degenerate ground states are represented as $\ket{\psi^{0}_{\mathrm{GS}}}\simeq\ket{\uparrow\downarrow\uparrow\downarrow\cdots}$ and $\ket{\psi^{1}_{\mathrm{GS}}}\simeq\ket{\downarrow\uparrow\downarrow\uparrow\cdots}$. In this case, the OSEE of $\mathrm{e}^{-\beta\hat{H}}$ at $\beta\gg 1$ is given by

[TABLE]

where we used the fact that the entanglement entropies of the two ground states are the same to obtain the third line. When the degeneracy is removed, the OSEE decreases by $\ln 2$ in the above equation, because the OSEE reduces to that of a single ground state.

Figure 3 shows the inverse-temperature dependence of the OSEE with $\Delta=3$, which is already shown in Fig. 2. The open black circles depict ${S_{A}^{\mathrm{OSEE\>(2)}}}+\ln 2$ at low temperature $(7.2\leq\beta)$ and the black dashed line shows the average of them. We see that if Eq. (17) holds approximately, the OSEE at low temperature saturates near the black dashed line.

III.3 Real-time evolution

We next consider the real-time evolution of the OSEE. In the following, we adopt sufficiently high temperatures where the above-mentioned artificial symmetry breaking does not occur.

Figure 4 (a) shows the real-time dependence of the OSEE for $\mathrm{e}^{-\left(\frac{\beta}{4}+\mathrm{i}t\right)\hat{H}}$ at $\beta=4$ and with $\Delta=0.25,0.5,0.75,1,1.6,1.8,2$. The OSEE in the gapped region increases less than that in the gapless region. This implies that the corresponding averaged OTOC decays faster in the gapless region.

In the gapless region ($\Delta=0.25,0.5,0.75,1$), the averaged OTOC exhibits the same behavior for the different values of $\Delta$. As we will argue in the next subsection, this is because the XXZ model in the gapless region belongs to the same universality class with the central charge $c=1$. In the gapped region, the energy gap $E_{\mathrm{G}}$ is written as an increasing function of $\Delta$: $E_{\mathrm{G}}\simeq 4\pi J\exp\quantity{-\frac{\pi^{2}}{2\sqrt{2\quantity(\Delta-1)}}}$ Baxter (1982); Takahashi (1999). From Fig. 4 (a), the increase of the OSEE is smaller if $E_{\mathrm{G}}$ is larger. This implies that the decay of the averaged OTOC becomes slower, as the parameter goes apart from the gapless region and the energy gap $E_{\mathrm{G}}$ becomes larger.

III.4 Role of the scale invariance

The numerical result in previous subsections implies that the system with the scale invariance scrambles fast. We will further investigate the role of scale invariance.

In the context of the 2D CFT, the Rényi-$n$ OSEE is represented by a correlation function of the twist operators on a cylinder. Let $L_{\mathrm{eff}}=\frac{\beta v_{s}}{2\pi}\cosh\frac{2\pi}{\beta}t$ be the effective length. The two-point function of the twist operator is given by

[TABLE]

where $\Delta_{n}=\frac{c}{12}(n-\frac{1}{n})$. The OSEE with the bipartition $\mathrm{A}=\left(-\infty,0\right]$ and $\mathrm{B}=\left(0,-\infty\right)$ is written as ${S_{A}^{\mathrm{OSEE}\>(n)}}\quantity[\mathrm{e}^{-\frac{\beta}{4}\hat{H}}]=\frac{1}{1-n}\ln\langle\mathcal{T}(u)\overline{\mathcal{T}}(v)\rangle_{\mathrm{cyl}}$, and the Rényi-$n$ OSEE is given by ${S_{A}^{\mathrm{OSEE\>(2)}}}=\frac{c}{2}\ln L_{\mathrm{eff}}+a$, where $a$ is a constant. Thus, we expect that the numerical results (even with different temperatures) of the OSEE are characterized by a single function that depends only on $L_{\mathrm{eff}}$ if the scale invariance emerges.

Figure 5 shows the effective-length dependence of the OSEE for the complex-time evolution operator $\mathrm{e}^{-\left(\frac{\beta}{4}+\mathrm{i}t\right)\hat{H}}$ in the gapless region ($\Delta=0.5$). The horizontal axis is $\ln L_{\mathrm{eff}}=\ln\quantity(\frac{\beta v_{s}}{2\pi}\cosh(\frac{2\pi}{\beta}t))$. We plot the Rényi-2 OSEE at $\beta=2,3,4,5,6$. For each parameter, we calculate the effective length dependence by fixing the inverse temperature and by changing the real time. As shown in Fig. 5, the data points almost collapse onto a single curve, which implies that the OSEE depends only on the effective length $L_{\mathrm{eff}}$. In the right side of the figure, however, we find that there are some deviations. We consider that these deviations are caused by the fact that the bond dimension $\chi$ used in our numerical calculation is not large enough. In Appendix C, we evaluate the numerical errors from the Trotter decomposition and the bond dimension, where the numerical errors are sufficiently small in the range of $\ln L_{\mathrm{eff}}\leq 4$.

We fit the numerical data points for $\beta=2,3,4,5,6$ against the fitting function ${S_{A}^{\mathrm{OSEE\>(2)}}}=\frac{c}{4}\ln L_{\mathrm{eff}}+a$, which is from the 2D CFT with fitting parameters $c$ and $a$. The fitting range is taken as $\ln L_{\mathrm{eff}}\leq 4$. We expect that the universal behavior of the correlation function appears when $L_{\mathrm{eff}}$ is sufficiently large.

We perform the fitting by changing the lower end of the fitting range within $\ln L_{\mathrm{eff}}\leq 4$, i.e., for intervals $[2.5,4]$, $[3,4]$, and $[3.5,4]$. The results are summarized in Table 1. The obtained values of $c$ are close to unity, which is consistent with the fact that the XXZ model with $\Delta=0.5$ at low temperature is a quantum critical system described by the 2D CFT with the central charge $c=1$. Thus, we conclude that the numerical results in the gapless region show universal behavior caused by the scale invariance.

IV Summary and discussion

We have investigated the relationship between the scale invariance at quantum criticality and the OTOC for the one-dimensional anti-ferromagnetic quantum spin-$1/2$ XXZ chain. Specifically, we have numerically calculated the time evolution of the OSEE for $\hat{U}(\beta,t)=\mathrm{e}^{-\left(\frac{\beta}{4}+\mathrm{i}t\right)\hat{H}}$ both in the gapless region and in the gapped region.

We have shown that the OSEE exhibits qualitatively different behaviors in the gapless and the gapped regions as demonstrated in Fig. 4. In the gapped region, the amount of increase of the OSEE becomes smaller as the system goes away from the gapless region, which implies that the decay of the averaged OTOC becomes slower. While these numerical results are for bipartition of the chain (i.e., $\mathrm{A}=\quantity{-\infty,\cdots,-1,0}$ and $\mathrm{B}=\quantity{1,2,\cdots\infty}$), we expect a similar behavior for general $\mathrm{A}^{\prime}$, $\mathrm{B}^{\prime}$ because of inequality Eq. (4). Furthermore, as shown in Fig. 5, the OSEE for $\hat{U}(\beta,t)$ in the gapless region at low temperature exhibits the similar behavior as that derived by the 2D CFT analytically, which is consistent with the fact that the 2D CFT is an effective theory of the XXZ model.

In this study, we have performed numerical simulations in the range of $\Delta\geq 0$. We argue that the averaged OTOC in $-1<\Delta<0$ behaves similarly to that in $0\leq\Delta\leq 1$, because the low-energy physics is effectively described by the 2D CFT of $c=1$ in all of the gapless regime $-1<\Delta\leq 1$ Korepin (1993); Giamarchi (2003). In $\Delta\leq-1$, on the other hand, it is known that the low-energy excitation above the ground state is gapped, where magnon bound states appear Fukuhara et al. (2013). Since the propagation velocity of free magnons is faster than that of magnon bound states, the decay of the OTOC in the early time regime with $\Delta<-1$ is dominated by the contribution of free magnons and thus similar to that in the gapped region $\Delta>1$.

On the basis of our numerical results, we have concluded that the scale invariance (or quantum criticality) is important for the fast decay of the OTOC. On the other hand, we did not observe that the OTOC in the XXZ model exhibits the exponential decay and saturates the MSS bound Maldacena et al. (2016a), in contrast to the SYK model Maldacena and Stanford (2016) and the 2D CFT Roberts and Stanford (2015). This is consistent with previous results Maldacena et al. (2016a), because the saturation of the MSS bound would come from the semiclassical limit such as the large central charge limit $c\gg 1$ of the 2D CFT, while this limit cannot be taken in our setup of the XXZ spin chain. Investigating whether the MSS bound can be saturated in a spin model with local interactions is a future issue. The quantum simulator can experimentally address this problem with ultracold atoms Gross and Bloch (2017) and superconducting qubits Wendin (2017).

Acknowledgements.

We are grateful to Shunsuke Furukawa, Youhei Yamaji, and Masatoshi Imada for valuable discussions. E.I. and T.S. are supported by JSPS KAKENHI Grant Number JP16H02211. E.I. is also supported by JSPS KAKENHI Grant Number JP15K20944.

Appendix A Inequality for OTOC

In this Appendix, we prove inequality (4). We first discuss the representation of the OTOC in the operator space Xu and Swingle . By using the notation of Eq. (5), we write the inner product between operators $\hat{O}_{\mathrm{A}}$ and $\hat{O}_{\mathrm{B}}$ as

[TABLE]

Also, the product $\hat{O}_{\mathrm{A}}\hat{O}_{\mathrm{B}}$ can be written as the following state:

[TABLE]

With the above notations, we represent the regularized OTOC. We now assume that $\hat{O}_{A}$ and $\hat{O}_{B}$ are Hermitian. We define the Heisenberg representation of $\hat{O}_{B}$ for the complex-time evolution as $\hat{O}_{B}(\beta,t)\coloneqq\mathrm{e}^{-\quantity(\frac{\beta}{4}-\mathrm{i}t)\hat{H}}\hat{O}_{B}\mathrm{e}^{-\quantity(\frac{\beta}{4}+\mathrm{i}t)\hat{H}}/\sqrt{Z}$. The OTOC is then written as

[TABLE]

Using the cyclic property of the trace, we also obtain

[TABLE]

We now prove inequality (4). A subsystem of A is denoted as $\mathrm{A}^{\prime}$. A basis operator $\hat{O}_{\mathrm{A}}$ can be decomposed into a basis operator of the subsystem $\mathrm{A}^{\prime}$ and that of $\mathrm{A}\setminus\mathrm{A}^{\prime}$ as

[TABLE]

With this decomposition, we obtain Xu and Swingle

[TABLE]

We then calculate the OTOC for the subsystems A and B as

[TABLE]

Finally, we take the average over the basis operator $\hat{O}_{\mathrm{A}^{\prime}}$ in the subsystem $\mathrm{A}^{\prime}$ while $\hat{O}_{\mathrm{A}\setminus\mathrm{A}^{\prime}}$ being fixed, and obtain

[TABLE]

where $\hat{I}_{\mathrm{A}^{\prime}}$ is the identity operator for the subsystem $\mathrm{A}^{\prime}$. By taking this average, the second term in the most right-hand side of Eq. (24) can be regarded as the expectation value of a reduced density operator $\rho_{\mathrm{A}^{\prime}}=\frac{1}{2^{a^{\prime}}}\ket{\hat{I}_{\mathrm{A}^{\prime}}}\bra{\hat{I}_{A^{\prime}}}$ with a state vector $\ket{[\hat{O}_{\mathrm{A}\setminus\mathrm{A}^{\prime}},\hat{O}_{\mathrm{B}}(\beta,t)]}$. Since the reduced density operator is positive semidefinite, this term should be non-negative. Thus, we obtain

[TABLE]

By taking the average over $\hat{O}_{\mathrm{B}}(\beta,t)$ and $\hat{O}_{\mathrm{A}\setminus\mathrm{A}^{\prime}}$, we obtain the following inequality for the averaged OTOC:

[TABLE]

By the decomposition of the subsystem B, we can perform the same procedure and finally obtain Eq. (4).

Appendix B Canonical form of iMPO

In parallel to the canonical form of iMPS defined in Ref. Orús and Vidal (2008), we define the canonical form of the iMPO in the following. For simplicity, we consider a one-dimensional lattice system composed of $d$-level local systems. We can define the iMPO with the single-site translational invariance symmetry by two tensors $\quantity{\Gamma,\lambda}$, where $\Gamma$ is a complex-valued tensor with four indices and $\lambda$ is a real-valued tensor with a single index. The elements of these tensors are represented by $\Gamma^{ij}_{\alpha\beta}$, and $\lambda$ is represented as a diagonal matrix with non-negative elements $\lambda_{\alpha}$. While the indices $\alpha$ and $\beta$ represent virtual degrees of freedom (virtual bond indices) and run within $1,\cdots,\chi$, the indices $i$ and $j$ represent physical degrees of freedom (physical bond indices) and run within $1,\cdots,d$.

In order to define the canonial form, we introduce the following tensors $R$ and $L$:

[TABLE]

The condition for the iMPO $\quantity{\Gamma,\lambda}$ to be the canonical form is written as

[TABLE]

With the canonical form, the Rényi-$n$ OSEE for the bipartition $\mathrm{A}=\quantity{-\infty,\cdots,-1,0}$ and $\mathrm{B}=\quantity{1,2,\cdots\infty}$ can be calculated as

[TABLE]

Appendix C Numerical error

In this Appendix, we discuss numerical errors in our calculation. We show that the errors from the Trotter decomposition and the bond dimensions are sufficiently small.

C.1 Imaginary-time evolution

Figure 6 shows the dependence on the width in the Trotter decomposition of the imaginary-time evolution of the OSEE in the gapless region $\Delta=0.5$ and the gapped region $\Delta=3$. The width in the Trotter decomposition is taken as $\Delta\beta=0.04,0.02,0.01,0.005$ and the bond dimension is fixed as $\chi=512$.

In the gapless region, the numerical errors from the Trotter decomposition are smaller than the point size. We estimate the numerical errors $\Delta S^{\mathrm{OSEE}}$ quantitatively from the difference between the maximum and minimum values for $\Delta\beta=0.04,0.02,0.01,0.005$. We find that $\Delta S^{\mathrm{OSEE}}\sim 0.008$ at $\beta=8$ and $\Delta S^{\mathrm{OSEE}}\sim 0.01$ at $\beta=15.2$.

In the gapped region, while the numerical errors from the Trotter decomposition are smaller than the point size in the high-temperature region with small $\beta$, they are not as small in the low-temperature region with large $\beta$. We consider that this error in the low-temperature region comes from the artificial symmetry breaking in the ground state, which is discussed in Sec. III.2. In the high-temperature region, we estimate the error of the OSEE in the same manner as the gapless region and we obtain $\Delta S^{\mathrm{OSEE}}\sim 0.001$ at $\beta=4$.

Figure 7 shows the $\chi$-dependence of the imaginary-time evolution of the OSEE in the gapless region $\Delta=0.5$ and the gapped region $\Delta=3$. The bond dimension is taken as $\chi=256,384,512$ and the width in the Trotter decomposition is fixed to $\Delta\beta=0.01$.

In the gapless region, the numerical errors from the bond dimensions are smaller than the point size. We estimate the numerical errors $\Delta S^{\mathrm{OSEE}}$ in the same manner as the case of the Trotter decomposition. We find that $\Delta S^{\mathrm{OSEE}}\sim 0.009$ at $\beta=8$ and $\Delta S^{\mathrm{OSEE}}\sim 0.006$ at $\beta=15.2$.

In the gapped region, while the numerical errors from the bond dimensions are smaller than the point size in the high-temperature region, they are not small in the low-temperature region. This behavior is due to the artificial symmetry breaking in the ground state. We estimate the error of the OSEE and obtain $\Delta S^{\mathrm{OSEE}}\sim 0.006$ at $\beta=4$.

C.2 Real-time evolution

When we calculate the real-time evolution at finite temperature, we first calculate $\hat{U}(\beta,0)=\mathrm{e}^{-\frac{\beta}{4}\hat{H}}$ by the imaginary-time evolution and then calculate the real-time evolution to obtain $\hat{U}(\beta,t)=\mathrm{e}^{-\quantity(\frac{\beta}{4}+\mathrm{i}t)\hat{H}}$. The width for the imaginary-time evolution is fixed as $\Delta\beta=0.005$ and the inverse temperature is $\beta=4$. Figure 8 shows the real-time evolution of the OSEE with the width in the Trotter decomposition $\Delta t=0.02,0.01,0.005,0.0025$. We fix the bond dimension as $\chi=512$.

In the gapless region, the numerical errors from the Trotter decomposition are smaller than the point size. The numerical errors of the OSEE $\Delta S^{\mathrm{OSEE}}$ are estimated as $\Delta S^{\mathrm{OSEE}}\sim 0.006$ at $t/\beta=0.2$ and $\Delta S^{\mathrm{OSEE}}\sim 0.006$ at $t/\beta=0.4$.

In the gapped region, though the numerical errors from the Trotter decomposition are not smaller than the point size, we have the result that the averaged OTOC decays faster in the gapless region than in the gapped region. The numerical errors are estimated as $\Delta S^{\mathrm{OSEE}}\sim 0.004$ at $t/\beta=0.2$ and $\Delta S^{\mathrm{OSEE}}\sim 0.01$ at $t/\beta=0.4$. When we consider $S^{\mathrm{OSEE}\>(n)}-S_{0}$ where $S_{0}$ is the OSEE at $t=0$, the difference between the OSEEs in the gapless and the gapped regions is approximately $0.08$ at $t/\beta=0.2$ and $0.34$ at $t/\beta=0.4$. The numerical errors of the OSEE from the Trotter decomposition are sufficiently smaller than the changes of $S^{\mathrm{OSEE}\>(n)}-S_{0}$ for our discussion.

Figure 9 shows the real-time evolution of the OSEE with $\chi=256,384,512$. We fix the width in the Trotter decomposition as $\Delta t=0.005$.

In the gapless region, the numerical errors from the bond dimensions are smaller than the point size. We estimate the numerical errors $\Delta S^{\mathrm{OSEE}}$ and find that $\Delta S^{\mathrm{OSEE}}\sim 0.006$ at $t/\beta=0.2$ and $\Delta S^{\mathrm{OSEE}}\sim 0.006$ at $t/\beta=0.4$.

In the gapped region, though the numerical errors from the bond dimensions are not smaller than the point size, the result does not change because the errors are sufficiently smaller than the change of the OSEE. The numerical errors are estimated as $\Delta S^{\mathrm{OSEE}}\sim 0.008$ at $t/\beta=0.2$ and $\Delta S^{\mathrm{OSEE}}\sim 0.03$ at $t/\beta=0.4$. The difference between $S^{\mathrm{OSEE}\>(n)}-S_{0}$ in the gapless and the gapped regions is approximately $0.08$ at $t/\beta=0.2$ and $0.34$ at $t/\beta=0.4$. Thus, the numerical errors of the OSEE from the bond dimensions do not affect our conclusion.

C.3 Role of scale invariance

We calculate the real-time evolution of the Rényi-2 OSEE with $\chi=256,384,512$, $\Delta\beta=0.005$ and $\Delta t=0.005$. Figure 10 shows the dependence of the OSEE on the effective length $L_{\mathrm{eff}}=\frac{\beta v_{s}}{2\pi}\cosh\frac{2\pi}{\beta}t$. Each data series corresponds to the same $\chi$ and contains data with different temperatures $\beta=2,3,4,5,6$. When the effective length is sufficiently large as $\ln L_{\mathrm{eff}}\sim 4$, the numerical errors from the bond dimensions become larger than the point size.

Figures 11 and 12 show $\ln L_{\mathrm{eff}}=\ln\quantity(\frac{\beta v_{s}}{2\pi}\cosh(\frac{2\pi}{\beta}t))$ dependence of the Rényi-2 OSEE of the complex-time evolution operator $\mathrm{e}^{-\left(\frac{\beta}{4}+\mathrm{i}t\right)\hat{H}}$ with $\chi=256,384$. We obtain the data by calculating real-time evolution with fixing $\beta$ for each $\beta=2,3,4,5,6$. For each bond dimension, we fit the numerical data of all $\beta$’s with a fitting function ${S_{A}^{\mathrm{OSEE\>(2)}}}=\frac{c}{4}\ln L_{\mathrm{eff}}+a$, which comes from the 2D CFT. The fitting variables are $c$ and $a$. We adopt three fitting ranges $[2.5,4]$, $[3,4]$, and $[3.5,4]$, where the numerical errors from the bond dimensions are sufficiently small. Table 2 and 3 show the results of the fitting for $\chi=256$ and $384$, respectively. The obtained value of $c$ is close to unity, which is consistent with the fact that the XXZ model with $\Delta=0.5$ at low temperature is a quantum critical system effectively described by the 2D CFT.