Temperature driven quenches in the Ising model: appearance of negative Rényi mutual information
TTemperature driven quenches in the Ising model:appearance of negative R´enyi mutual information
M´arton Kormos
BME-MTA Statistical Field Theory Research Group, Institute of Physics, BudapestUniversity of Technology and Economics, H-1111 Budapest, HungaryE-mail: [email protected]
Zolt´an Zimbor´as
Dahlem Center for Complex Quantum Systems, Freie Universit¨at Berlin, 14195Berlin, GermanyWigner Research Centre for Physics, Hungarian Academy of Sciences, P.O. Box 49,H-1525 Budapest, HungaryE-mail: [email protected]
Abstract.
We study the dynamics of the transverse field Ising chain after a localquench in which two independently thermalised chains are joined together and are leftto evolve unitarily. In the emerging non-equilibrium steady state the R´enyi mutualinformation with different indices are calculated between two adjacent segments of thechain, and are found to scale logarithmically in the subsystem size. Surprisingly, forR´enyi indices α > α > α < . Interestingly, even for α >
1. Introduction
In the past decade, the study of correlations between subsystems in lattice modelsand in field theories has allowed for a deeper understanding of the physics of many-bodysystems, in particular in relation to quantum criticality [1, 2, 3, 4, 5], equilibration [6, 7],and topological order [8, 9]. In pure states (e.g. in ground states), the correlation betweentwo complementary subsystems is entirely quantum mechanical and can be measured bythe entanglement entropy. The ground state entanglement entropy for one dimensionalgapped local Hamiltonians was proved to obey an area law [10, 11], while for critical a r X i v : . [ c ond - m a t . s t a t - m ec h ] J un emperature driven quenches in the Ising model H A ⊗ H B , the correlation between subsystems A and B can be characterised by themutual information (MI), I ( A : B ) = S ( ρ A ) + S ( ρ B ) − S ( ρ AB ) , (1)where ρ A and ρ B are the reduced density matrices of the subsystems A and B , and S denotes the von Neumann (or entanglement) entropy, S ( ρ ) = − Tr ρ log ρ . (2)The MI has many nice properties. It is positive due to the subadditivity of thevon Neumann entropy, and is zero if and only if ρ AB = ρ A ⊗ ρ B , i.e. when thestate is uncorrelated. An operational interpretation of the MI is that it measures(asymptotically) the minimal amount of noise needed to erase the correlation in thestate by turning it into a product state [20]. This is related to the fact that the MI isequal to the relative entropy (or quantum Kullback–Leibler divergence) between ρ AB and ρ A ⊗ ρ B . The relative entropy, defined as D ( ρ || σ ) = Tr ρ (log ρ − log σ ) , (3)is a measure of distinguishability between two quantum states ρ and σ . There has beenan increasing activity on the relative entropy in field theory, see [21] and referencestherein. The mentioned relation to the MI can be shown by the following standardderivation D ( ρ AB || ρ A ⊗ ρ B ) = Tr H A ⊗H B ρ AB (log ρ AB − log ρ A ⊗ ρ B )= − S ( ρ AB ) − Tr H A ⊗H B ( ρ AB log ρ A ⊗ B + ρ AB log 1l A ⊗ ρ B )= − S ( ρ AB ) − Tr H A ρ A log ρ A − Tr H B ρ B log ρ B = I ( A : B ) . Alternatively, the mutual information can also be charaterised as the following minimum I ( A : B ) = min σ B D ( ρ AB || ρ A ⊗ σ B ) , (4)where σ B is any density matrix on the Hilbert space H B of subsystem B .It was shown that in finite temperature Gibbs states of local Hamiltonians a strictarea law holds for the mutual information [22, 23]. Until now, the only examples of aviolation of the area law in the MI outside the zero temperature regime was found toappear in non-equilibrium steady states (NESS) of spin chains [24]. These states can bewritten as Gibbs states of infinite-range Hamiltonians, and thus the theorems of Refs.[22, 23] which use the local structure of the interactions do not apply. emperature driven quenches in the Ising model α is defined as S ( α ) ( ρ ) = 11 − α log Tr( ρ α ) . (5)Note that in the α → α indices (and with α >
1) appear as natural quantities in conformal fieldtheory through the replica approach [14, 28, 29, 30]. These quantities are also easier tocalculate than the von Neumann entropy in various analytical and numerical settings[31, 32, 33, 34], and one can use them to detect criticality [1] and topological order [35].Moreover, the experimental determination of R´enyi entropies also seems more feasible[36, 37, 38, 39, 40]. Following this line of studies, as a natural generalisation of Eq. (1),also the R´enyi mutual information with index α was introduced as I ( α ) ( A : B ) = S ( α ) ( ρ A ) + S ( α ) ( ρ B ) − S ( α ) ( ρ AB ) . (6)The R´enyi mutual information has been shown to exhibit universal scaling behaviourin ground, excited and thermal states [41, 42, 43, 44, 45]. Furthermore, also for post-measurement states [46], in holographic settings [47], and non-equilibrium scenarios [48]certain universal features show up.However, when discussing the extensive use of R´enyi MI in many-body physics,it should be mentioned that, unlike the standard MI, this quantity has in general nooperational meaning and may even be negative. Thus, in quantum information theorya different R´enyi generalisation of MI is used. First, the relative entropy is generalisedby defining the α -R´enyi divergences [49, 50], D ( α )1 ( ρ || σ )= 1 α − (cid:0) ρ α σ − α (cid:1) , D ( α )2 ( ρ || σ )= 1 α − (cid:16) σ − α α ρσ − α α (cid:17) α , (7)which are used, e.g., in state discrimination theory [51]. Building on these divergences,by generalising Eq. (4), a regularised α -R´enyi mutual information can be introduced as I ( α ) j ( A : B ) = min σ B D ( α ) j ( ρ AB || ρ A ⊗ σ B ) (8)for j = 1 ,
2. These quantities are not only positive by definition, but have also othernice properties including operational interpretations [52, 53, 54, 55].Let us return to the R´enyi MI defined by Eq. (6). Despite the mentioned generalproblems with this quantity, for particular families of states it was found to be useful.For bosonic Gaussian states it was proved that the 2-R´enyi mutual information ispositive, or equivalently that the 2-R´enyi entropy satisfies the subadditivity condition[56]. This allowed for the introduction of new types of correlation measures, e.g., steeringquantifiers, which had no counterparts among quantities based on the conventional von emperature driven quenches in the Ising model I ( α ) isnegative for α >
2. Transverse field Ising model: temperature driven quench andnon-equilibrium steady state
The system we study in this work consists of two half-infinite chains thermalisedat different temperatures and brought to contact at time zero. This setup belongsto a more general scheme that is called in the literature “cut and glue quench” or“partitioning approach” which also includes the case of different chemical potentials(or magnetisation) on each side [59]. The non-equilibrium steady state was constructedfor the two-temperature case in the XX and XY spin chains in [60, 61, 62]. In [63, 64]the spatial profile of the magnetisation density and current was determined in the XXspin chain, while the Ising spin chain was studied in [65, 66]. These systems can bemapped to free spinless fermions unlike the integrable XXZ spin chain investigated in[67, 68, 69]. Continuum theories were investigated as well, including the free bosonic [70]and fermionic [71, 72] systems as well as integrable quantum field theories [73, 74, 75].There is a growing body of results in conformal field theories where the energy density,the full distribution of the current and fluctuation relations in the NESS were obtained emperature driven quenches in the Ising model
The Hamiltonian of the transverse field Ising spin chain of length N is H = − N − (cid:88) j =1 σ xj σ xj +1 − N (cid:88) j =1 hσ zj , (9)where σ αj are the Pauli matrices and we consider open boundary conditions. Bythe Jordan–Wigner transformation, c j = (cid:81) j − k =1 ( − σ zk ) σ − j , c † j = (cid:81) j − k =1 ( − σ zk ) σ + j , theHamiltonian is mapped on that of free spinless fermions: H = − N − (cid:88) j =1 (cid:104) c † j c j +1 + c † j +1 c j + c † j c † j +1 + c j +1 c j (cid:105) − h N (cid:88) j =1 (cid:18) c † j c j − (cid:19) , (10)where { c j , c † k } = δ j,k , { c j , c k } = { c † j , c † k } = 0 . It is useful to introduce the Majoranafermion operators a j − = c j + c † j , a j = i ( c j − c † j ) . (11)The Hamiltonian is a bilinear form which can be diagonalised by a linear transformationleading to the fermionic mode operators η k = 12 N (cid:88) j =1 [ φ k ( j ) a j − − iψ k ( j ) a j ] , η † k = 12 N (cid:88) j =1 [ φ k ( j ) a j − + iψ k ( j ) a j ] , (12) a j − = (cid:88) k φ k ( j )( η † k + η k ) , a j = − i (cid:88) k ψ k ( j )( η † k − η k ) (13)with the functions φ k ( j ) = A k sin( kj − θ k ) , (14a) ψ k ( j ) = − A k sin( kj ) , (14b)where 0 < θ k < π is the Bogoliubov angle satisfyingtan θ k = sin kh + cos k (15) emperature driven quenches in the Ising model A − k = (cid:80) Nj =1 sin ( kj ) is the normalisation. The modes satisfy { η k , η † k (cid:48) } = δ k,k (cid:48) , { η k , η k (cid:48) } = 0 , and in terms of them the Hamiltonian reads H = (cid:88) k ε k η † k η k + const. (16)with the dispersion relation ε k = √ h cos k + h . (17)In finite volume, the “momentum” k can only take quantised values according to thecondition k ( N + 1) − θ k = nπ , n ∈ Z . (18) Our initial state corresponds to two independent, disjoint chains of length N thermalised at different temperatures T L and T R , so the initial density matrix is ρ = ρ L ( T L ) ⊗ ρ R ( T R ) . (19)At time t = 0 the two halves are joined and let evolve by the Hamiltonian H of thechain of length 2 N. In other words, we turn on the coupling between site 0 and site 1, H = H L + H R − σ x σ x = (cid:88) k ε k γ † k γ k + const. , (20)where H L/R are the Hamiltonian of the left and right chain of length N, respectively,and the γ k are the mode operators that diagonalise H. Let us denote the mode functions of H R by φ q and ψ q , then the mode functionsof H L are φ L q ( j ) = ψ R q (1 − j ) and ψ L q ( j ) = φ R q (1 − j ) . As the first site of the full chainof length 2 N has index − N + 1 , the eigenfunctions of H are ϕ k ( j ) = φ k ( j + N ) and χ k ( j ) = ψ k ( j + N ) , where the momenta { k n } are quantised with 2 N instead of N in Eq.(18). Note that the functional form of ε k and θ k are the same in all cases (we do notquench the Ising interaction or the transverse field). The Majorana operators on the fullchain are related to the modes by a j − = (cid:88) k ϕ k ( j )( γ † k + γ k ) , a j = − i (cid:88) k χ k ( j )( γ † k − γ k ) . (21)In order to compute the time evolution of correlation functions of spin operators,we first need to compute the building blocks given by the Majorana correlations (cid:104) a n ( t ) a m ( t ) (cid:105) . The time evolved operators in the Heisenberg picture are a n ( t ) = (cid:88) j (cid:104) a j | a n ( t ) (cid:105) a j , (22) emperature driven quenches in the Ising model (cid:104) a j − | a n − ( t ) (cid:105) = (cid:88) k ϕ k ( j ) ϕ k ( n ) cos( ε k t ) , (23a) (cid:104) a j | a n ( t ) (cid:105) = (cid:88) k χ k ( j ) χ k ( n ) cos( ε k t ) , (23b) (cid:104) a j − | a n ( t ) (cid:105) = −(cid:104) a n | a j − ( t ) (cid:105) = (cid:88) k ϕ k ( j ) χ k ( n ) sin( ε k t ) . (23c)In the infinite volume limit, N → ∞ , the sum over k turns into an integral. Droppinghighly oscillating terms in the integrands we obtain (cid:104) a j − | a n − ( t ) (cid:105) = (cid:104) a j | a n ( t ) (cid:105) = (cid:90) π − π d k π ˜ ϕ ∗ k ( j ) ˜ ϕ k ( n ) cos( ε k t ) , (24a) (cid:104) a j − | a n ( t ) (cid:105) = −(cid:104) a n | a j − ( t ) (cid:105) = (cid:90) π − π d k π ˜ ϕ ∗ k ( j ) ˜ χ k ( n ) sin( ε k t ) , (24b)where the infinite volume mode functions are˜ ϕ k ( j ) = e − ikj + iθ k , ˜ χ k ( j ) = − e − ikj . (25)The time dependent Majorana two-point functions can be written as (cid:104) a n ( t ) a m ( t ) (cid:105) = (cid:88) j,l (cid:104) a j | a n ( t ) (cid:105)(cid:104) a l | a m ( t (cid:48) ) (cid:105) (cid:104) a j a l (cid:105) . (26)Clearly, the initial correlations will be non-zero only if j, l ≥ j, l ≤ , so thecorrelation function (26) splits into two parts corresponding to the contributions of theleft and right half chains. We now rewrite the Majorana operators in terms of the modeoperators η q diagonalising the left and right half chains and use for each half chain (cid:104) η † q η † q (cid:48) (cid:105) = (cid:104) η q η q (cid:48) (cid:105) = 0 and (cid:104) η † q η q (cid:48) (cid:105) = δ q,q (cid:48) f q , where f q = 11 + e ε q /T L/R (27)is the thermal Fermi–Dirac distribution function. Exploiting the completeness of themode functions, we arrive at (cid:104) a j − a l − (cid:105) = (cid:104) a j a l (cid:105) = δ j,l , j, l ≥ j, l ≤ , (28a) (cid:104) a j − a l (cid:105) = −(cid:104) a l a j − (cid:105) = − i (cid:88) q φ R q ( j ) ψ R q ( l )(1 − f R q ) , j, l ≥ , (28b) (cid:104) a j − a l (cid:105) = −(cid:104) a l a j − (cid:105) = − i (cid:88) q φ L q ( j ) ψ L q ( l )(1 − f L q ) , j, l ≤ . (28c)In the N → ∞ limit these become integral expressions. emperature driven quenches in the Ising model δ n,m in Eq. (26), resulting in (cid:104) a n ( t ) a m ( t ) (cid:105) = δ n,m + ∞ (cid:88) j,l = −∞ (cid:104) (cid:104) a j − | a n ( t ) (cid:105)(cid:104) a l | a m ( t ) (cid:105)−(cid:104) a j − | a m ( t ) (cid:105)(cid:104) a l | a n ( t ) (cid:105) (cid:105) (cid:104) a j − a l (cid:105) , (29)where both the coefficients and the initial correlations are written in the infinite N limitin the form of integrals and any explicit dependence on N disappeared. For numericalsimulations, however, we use the finite N expressions involving finite sums. The non-equilibrium steady state corresponds to the limit t → ∞ in Eq. (29) with n, m fixed. In this limit, the gradients of all observables tend to zero resulting in atranslationally invariant state. Expression (29) can be greatly simplified in this limit[83], which leads to the asymptotic correlations first derived in [61], (cid:104) a n − ( t ) a m ( t ) (cid:105) NESS = i (cid:90) π − π d k π (cid:0) − f R k − f L k (cid:1) e iθ k e ik ( m − n ) , (30a) (cid:104) a n − ( t ) a m − ( t ) (cid:105) NESS = (cid:104) a n ( t ) a m ( t ) (cid:105) NESS = δ n,m + (cid:90) π − π d k π ( f R k − f L k ) e ik ( m − n ) sgn( k ) . (30b)Thus, the NESS is a fermionic Gaussian state defined by the covariance matrixcorresponding to the correlation functions (30). The Gaussianity of the NESS impliesthat it can be interpreted as a Gibbs state of an effective quadratic Hamiltonian, ρ NESS = 1 Z exp( − β H eff ) , (31)where Z = Tr[exp( − β H eff )] and we chose β = (cid:16) T L + T R (cid:17) . Since ρ NESS must commutewith the original Hamiltonian generating the dynamics, the effective Hamiltonian canonly be a sum of conserved charges of the transverse field Ising model H eff = ∞ (cid:88) n =0 µ + n I + n + µ − n I − n , (32)where the charges are given as [84] I + n = i (cid:88) j a j ( a j +2 n +1 + a j − n +1 ) − h a j ( a j +2 n − + a j − n − ) , (33) I − n − = − i (cid:88) j ( a j a j +2 n + a j − a j +2 n − ) . (34)A state that is the exponential of a linear combination of conserved chargescorresponding to a given Hamiltonian is usually called a Generalized Gibbs Ensemble emperature driven quenches in the Ising model h = 1 critical case, one can immediately determine the coefficients µ ± n from the expectation values (30), µ + n = δ n, , µ − n = 16 πβ (cid:18) T L − T R (cid:19) n + 14 n + 8 n + 3 , (35)which implies that H eff decays algebraically, which also remains true in the h (cid:54) = 1 case.In summary, the NESS of the transverse field Ising model can be written in a GGE-likeform but with a long ranged effective Hamiltonian.
3. R´enyi mutual information in the NESS
In this section we show that in the NESS the R´enyi mutual information of twoadjacent intervals of length L has logarithmic dependence on L and present analyticalresults for the prefactor of the logarithm for indices α = 1 , , , , m . The analyticexpressions are compared with results obtained by numerical evaluation of the mutualinformation. Furthermore, we also study R´enyi MI of higher index and study itsdependence on the index and on the parameters of the system. From the results of Section 2, in particular from Eqs. (30), one can immediatelydetermine the Majorana covariance matrix Γ x,y = i (cid:104) [ a x , a y ] (cid:105) of the NESS. Due totranslational invariance, the covariance matrix is a block-Toeplitz matrix, and thus canbe expressed as the Fourier transform of a 2 × symbol ofthe block-Toeplitz matrix, (cid:32) Γ n − , m − Γ n − , m Γ n, m − Γ n, m (cid:33) = (cid:90) π − π d k π e ik ( m − n ) Λ( k ) . (36)For the transverse Ising NESS, the symbol is given asΛ( k ) = (cid:32) i ( f R k − f L k ) sgn( k ) ( f L k + f R k − e iθ k − ( f L k + f R k − e − iθ k i ( f R k − f L k ) sgn( k ) (cid:33) , (37)where θ k and f R/L k are defined in Eqs. (15) and (27), (17), respectively. As the NESS is aGaussian state, the von Neumann and R´enyi entropies of a subsystem of L consecutivespins can be calculated from the eigenvalues ± iλ ( L ) j of the 2 L × L reduced covariancematrix Γ L through the formula S ( α ) L = L (cid:88) j =1 s ( α ) (cid:16) λ ( L ) j (cid:17) , (38)where s ( α ) ( λ ) = 11 − α log (cid:20)(cid:18) λ (cid:19) α + (cid:18) − λ (cid:19) α (cid:21) when α (cid:54) = 1 , (39) s (1) ( λ ) = − (cid:18) λ (cid:19) log (cid:18) λ (cid:19) − (cid:18) − λ (cid:19) log (cid:18) − λ (cid:19) . (40) emperature driven quenches in the Ising model S ( α ) L = L (cid:88) j =1 s ( α ) ( λ ( L ) j ) = 12 πi (cid:73) C d λ s ( α ) ( λ ) L (cid:88) j =1 λ − λ ( L ) k = 12 πi (cid:73) C d λ s ( α ) ( λ ) d ln D L ( λ )2 d λ , (41)where D L ( λ ) = det( λ − i Γ L ) and the contour C in the complex plane is encirclingthe real interval [ − , L can be calculated as I ( α ) L = 2 S ( α ) L − S ( α )2 L usingthe above contour integration. A similar calculation was done for the NESS of the XXchain in Ref. [24]. The big difference between the two cases is that instead of a simpleToeplitz matrix the covariance matrix is of block-Toeplitz type. The rather lengthycalculation is delegated to Appendix A and we just state the results here. The R´enyimutual information can be shown to have a logarithmic asymptotics in the subsystemsize, I ( α ) L = σ ( α ) log L + const. , (42)and for α = 1 , , , , m the prefactor of the logarithmic term can be explicitly obtained.Introducing the notations a = 1 e − (1+ h ) /T R + 1 , b = 1 e − (1+ h ) /T L + 1 , (43a) a = 1 e − (1 − h ) /T R + 1 , b = 1 e − (1 − h ) /T L + 1 , (43b)and defining η ( w ) = (cid:40) πi log( w ) when arg( w ) ∈ [0 , π ) , − πi log( w ) when arg( w ) ∈ [ − π, , (44) emperature driven quenches in the Ising model σ (1) = 12 π (cid:88) i =1 (cid:20) a i Li (cid:18) a i − b i a i (cid:19) +(1 − a i )Li (cid:18) b i − a i − a i (cid:19) + b i Li (cid:18) b i − a i b i (cid:19) +(1 − b i )Li (cid:18) a i − b i − b i (cid:19)(cid:35) , (45a) σ (2) = − − π Re (cid:34) (cid:88) j =1 log (cid:18) − a j + 1 + i b j − − i (cid:19) + η (cid:18) a j − − i b j − − i (cid:19)(cid:35) , (45b) σ (3) = − − π Re (cid:34) (cid:88) j =1 log (cid:32) − a j + 1 + i/ √ b j − − i/ √ (cid:33) + η (cid:32) a j − − i/ √ b j − − i/ √ (cid:33)(cid:35) , (45c) σ (4) = − − π Re (cid:34) (cid:88) j =1 log (cid:18) − a j +1+ i tan π b j − − i tan π (cid:19) + η (cid:18) a j − − i tan π b j − − i tan π (cid:19) + log (cid:18) − a j +1+ i tan π b j − − i tan π (cid:19) + η (cid:18) a j − − i tan π b j − − i tan π (cid:19)(cid:35) , (45d) σ (2 m ) = − m − m − − π (2 m − m − (cid:88) k =1 2 (cid:88) j =1 Re (cid:34) log (cid:32) − a j +1+ i tan (2 k − π m +1 b j − − i tan (2 k − π m +1 (cid:33) + η (cid:32) a j − − i tan (2 k − π m +1 b j − − i tan (2 k − π m +1 (cid:33)(cid:35) , (45e)where Li ( w ) denotes the dilogarithm function. Due to the structure of Eq. (45e), wecan also get the analytical form of the R´enyi prefactor in the α → ∞ limit: σ ( ∞ ) = − − π (cid:88) j =1 (cid:90) π − π d θ (cid:20) log (cid:18) − a j +1+ i tan θ b j − − i tan θ (cid:19) + η (cid:18) a j − − i tan θ b j − − i tan θ (cid:19)(cid:21) . (45f)Let us make a few comments about these results. Equality of the initialtemperatures, T L = T R , implies a j = b j and the above formulas give σ ( α ) = 0 forall α. Moreover, the expressions (45) are invariant under the transformation( h, T L , T R ) −→ ( h − , h/T L , h/T R ) . (46)This is a simple manifestation of the Kramers–Wannier duality, which in the Ising casecan be regarded as a “half-shift” transformation, i.e. a x → a x +1 (however, we emphasisethat this transformation is non-local and does not strictly leave the R´enyi MI invariant,but the change can be only manifest in the subleading terms). Finally, let us note that for α > emperature driven quenches in the Ising model L I H Α L (a) h = 0 . , T L = 0 . , T R = 5 L I H Α L (b) h = 0 . , T L = 0 . , T R = 5 Figure 1.
R´enyi mutual information I ( α ) as a function of the subsystem size L for(from top to bottom) α = 1 , , , h = 0 . T L = 0 . , T R = 5 on(a) linear and (b) logarithmic scale. Numerical results are plotted in (a) dots and (b)continuous lines while the dashed lines show the function σ ( α ) log L + const., wherethe analytic results for σ ( α ) given in Eqs. (45a)-(45d) and the constant is adjusted byhand. L I H Α L (a) h = 0 . , T L = 0 . , T R = 1 . - - - - - - - - L I H Α L (b) h = 10 , T L = 1 , T R = 15 Figure 2.
R´enyi mutual information I ( α ) as a function of the subsystem size L in the NESS at (a) h = 0 . , T L = 0 . , T R = 1 . α = 4 , , , , , h = 10 , T L = 1 , T R = 15 for (from bottom to top) α = 8 , , , , , . Continuous lines (a) and dots (b) are numerical results andthe dashed lines show the analytic result (45) with an adjusted additive constant.
We check our analytic expressions in Eq. (45) by comparing them to numericalresults obtained by computing the sums (38) over the eigenvalues of the covariancematrix (36). This is shown in Figs. 1 and 2 where dots and continuous lines are numericalresults and the dashed lines show the analytic expression, I ( α ) = σ ( α ) log L + const. , withthe constant shift adjusted by hand. We find excellent agreement in all cases.A surprising feature of the exact prefactors (45) is that for α ≥ emperature driven quenches in the Ising model I H Α L (a) h = 0 . , T L = 0 . , T R = 0 . ‰‰‰‰‰ ‰ ‰ ‰ ‰ ‰‰‰ ‰ ç ç ç ç ç ç ç ç ç
10 100 1000 - - - - - Α Σ H Α L (b) h = 0 . , T L = 0 . , T R = 1 . Figure 3. (a) R´enyi MI I ( α ) as a function of L for α = 1 . , , . σ ( α ) on the R´enyi index α. The analytic result (45e) is shown inempty circles, numerical fits are shown in crosses, and the dashed line indicates σ ( ∞ ) in (45f). - T R Σ H Α L (a) h = 0 . , T L = 0 . - - - h Σ H Α L (b) T L = 0 . , T R = 4 Figure 4.
Prefactors σ ( α ) with α = 2 , , , , ∞ (from top to bottom) (a) as a functionof the right temperature T R for h = 0 . T L = 0 .
1; (b) as a function of h for T L = 0 . ,T R = 4 on log-linear scale. In (a) the dashed line represents σ ( ∞ ) at T R = ∞ . monotonically decreases with the size L of the subsystems for large enough L. As aconsequence, the R´enyi MI can be negative and arbitrarily large in absolute value. Thisis shown in Fig. 2 both for a ferromagnetic and a paramagnetic case. The two cases aredual to each other under the transformation (46), so the leading logarithmic contributionis the same but the constant shift and the subleading terms in general are different. For h > , the scaling form I ( α ) = σ ( α ) log L + const. is reached at smaller values of L andthe corrections to it are much smaller than for h < . In the latter case we also findan even-odd oscillating behaviour in L which however decays as L is increased. Fig. 2bshows an example where the R´enyi MI is not only decreasing but it is also negative.We have shown that σ (3) can take negative values, but it would be interesting todetermine for which values of α can σ ( α ) be negative. For non-integer R´enyi index α we could only study this question numerically; and based on these investigations, we emperature driven quenches in the Ising model - - T L Σ H Α L L I H L H L L Figure 5.
Low temperature behaviour of the R´enyi mutual information at the criticalpoint, h = 1. (a) Temperature dependence of the prefactor σ ( α ) for α = 1 , , , T R = 2 T L . (b) R´enyi mutual information I (4) ( L ) as a function ofthe subsystem size L in a thermal Gibbs state of T = 0 . T = 0 . T L = 0 . , T R = 0 . conjecture that for α > h, T L , T R such that σ ( α ) < , while for α < α = 1 . , , . h = 0 . , T L = 0 . , T R = 0 . . For α = 1 . I ( α ) isincreasing while for α = 2 . I ( α ) ( L ) converges to a limiting function as α → ∞ . Note that theapproach is not monotonic, for example, in Fig. 2b, I (8) ( L ) < I (16) ( L ) < · · · < I ( ∞ ) ( L )
4. Time evolution of the R´enyi mutual information
So far we have focused on the R´enyi MI in the non-equilibrium steady state. Anotherand much more complicated aspect is the dynamics leading to the NESS. For the Isingmodel, the evolution of the von Neumann and R´enyi entropies have already been studiedfor global [87, 88] and local quenches [89, 90, 91]. We continue this line of investigationsby asking how the R´enyi MI evolves in time after joining the two halves of the systemand how it reaches its stationary value.We do not attempt any analytic derivation here but resort to numericalinvestigations using time-evolution equations (23) and (26). Some representative resultsin the paramagnetic phase are shown in Fig. 6 for interval length L = 40 . Note thatthere are cases when, quite oddly, the R´enyi MI decreases after joining the two halfchains. We find that similarly to the XX model [24], after an initial transient the MIevolves logarithmically in time up to t ≈ L : I ( α ) ( t ) ≈ ˜ σ ( α ) log t + const. (48)In our normalisation the maximal quasiparticle velocity is v max = max(d ε ( k ) / d k ) = 1 for h > , so t = L is the time necessary for the fastest quasiparticles to fly through and leavethe interval. For t > L there is a decay to the NESS value which is indicated by the dottedhorizontal lines in Fig. 6. It is computed using the formula I ( α ) L = σ ( α ) log L + const.with the analytic prefactors σ ( α ) and the constant adjusted by hand as in Fig. 2b.Based on our numerical findings we conjecture that the prefactor of log t is equalto the prefactor of the log L term in the NESS, that is,˜ σ ( α ) = σ ( α ) (conjecture) . (49)Proving this equality is beyond the scope of our paper, but our numerical results stronglysupport this conjecture. In Fig. 6 we plot in dashed line the function (48) using (49)and our analytic results for σ ( α ) . Similarly to the NESS fits, the constant is adjusted by emperature driven quenches in the Ising model t I H L (a) α = 1 t I H L (b) α = 2 - - - - - - t I H L (c) α = 4 - - - - - - - - t I H ¥ L (d) α = ∞ Figure 6.
Time evolution of the R´enyi MI for α = 1 , , , ∞ at h = 10 , T L = 1 , T R = 15and L = 40 . Numerical results are shown in solid line, the dashed lines correspondto I ( α ) = σ ( α ) log( t ) + const. with the constants adjusted by hand. The horizontaldotted lines are the NESS result I ( α ) = σ ( α ) log( L ) + const. where the constants weredetermined in Fig. 2b. hand. The agreement between the conjectured expression and the numerical results isexcellent.In Fig. 7 we present similar results in the ferromagnetic phase h < L = 60 . Herethe amplitude of the oscillations are stronger than in the paramagnetic phase makingthe analysis more difficult. The agreement with Eqs. (48) and (49) is still satisfactory.Interestingly, for h < h < v max = h, so the expected time where the logarithmic behaviour breaks down is t = L/h.
In Fig. 7 we find, however, that for large R´enyi index α Eq. (48) ceases to hold alreadyfor t ≈
500 instead of t = 60 / . . One is tempted to speculate that the higherindex R´enyi entropies may be more sensitive to the finite size of the quasiparticles.However, the physical reason behind this behaviour is presently unclear and deservesfurther study.Decreasing T R while keeping h and T L fixed the logarithmic behaviour (48) graduallydisappears for h < (cid:28) t (cid:28) L/h but this scaling regime is pushed towards larger emperature driven quenches in the Ising model t I H L (a) α = 1 t I H L (b) α = 2 t I H L (c) α = 32 t I H ¥ L (d) α = ∞ Figure 7.
Time evolution of the R´enyi MI for α = 1 , , , ∞ at h = 0 . , T L =0 . , T R = 10 and L = 60 . Numerical results are shown in solid line, the dashed linecorresponds to I ( α ) = σ ( α ) log( t ) + const. with the constants adjusted by hand. Thehorizontal dotted lines are the NESS result I ( α ) = σ ( α ) log( L ) + const. where theconstants were determined similarly to Fig. 2a. times, which requires larger intervals that are more difficult to study numerically.
5. Discussion and outlook
In this work we studied the von Neumann and the R´enyi mutual informationbetween two touching intervals of length L at the edges of two half infinite quantumIsing spin chains thermalised at different temperatures and subsequently glued together.Asymptotically a non-equilibrium steady state (NESS) is formed around the junction.We showed that in the NESS all the different types of MI depend logarithmically onthe length of the intervals in the leading order, I ( α ) = σ ( α ) log( L ) + const. We derivedclosed form exact analytic expressions (45) for the prefactor σ ( α ) for α = 2 m for all m = 0 , , , . . . as well as for α = 3. We found that the dependence on α is notmonotonic (c.f. Fig. 3b). Taking m to infinity allowed us to study the α → ∞ limitwhere we found a simple analytic expression (45f) for σ ( ∞ ) . We compared our analyticresults to numerical calculations in finite systems.Our most interesting finding is that the R´enyi MI can assume negative values. We emperature driven quenches in the Ising model α = 2 is a threshold value, that is I ( α ) can be negative for α > α ≤ − R´enyi MI is always positive as was found for bosonic Gaussian states [56].We also studied the non-equilibrium time evolution of the R´enyi MI after joining thetwo chains. It was shown that in certain cases the R´enyi MI can decrease in the courseof the non-equilibrium time evolution. Based on our numerical results, we conjecturethat there is a time domain after the initial transient and before saturation takes placewhere the MI evolves logarithmically in time, I ( α ) = ˜ σ ( α ) log( t ) + const. Moreover, weconjecture that the prefactor of the logarithmic term coincides with the prefactor ofthe logarithmic dependence on the subsystem size in the NESS, σ ( α ) = ˜ σ ( α ) . The samebehaviour was found in [24] for the XX spin chain. The proof of this statement is left forfuture work. Furthermore, we also made the observation that the logarithmic evolutionof the MI stops earlier for higher R´enyi indices than what a naive quasiparticle picturewould suggest.It is natural to expect that the negativity of the R´enyi MI can also be observedin other systems and physical situations. The most probable candidates are current-carrying non-equilibrium steady states in other settings. It would also be interesting tostudy spin chains that cannot be mapped to free fermions, e.g., the integrable XXZ spinchain.Finally, it would be worthwhile to investigate in this context the alternative R´enyiMI defined in terms of the R´enyi divergences by Eq. (8). Studying these quantities wouldnot only have natural consequences in the quantum information task of discriminatingbetween many-body states (see, e.g. [92]), but, presumably, would also yield a new toolfor understanding many-body correlations.
Acknowledgment
We would like to thank P. Calabrese, V. Eisler, M. Mezei, and G. Szirmai foruseful discussions and valueble inputs. M.K. acknowledges funding from a “Pr´emium”Postdoctoral Grant of the Hungarian Academy of Sciences and was partially supportedby NKFIH grant no. K 119204. Z.Z. was supported by the DFG (CRC183, EI 519/9-1,EI 519/7-1) and the ERC (TAQ). Z.Z. would also like to thank the Simons Center forGeometry and Physics for hospitality, where some of the work has been carried out. emperature driven quenches in the Ising model Appendix A: Analytical calculation of the R´enyi mutual informationasymptotics
As stated in Section 3, one can calculate the R´enyi MI through the contourintegral representation (41) and by the use of a generalisation of the Fisher–Hartwigconjecture. Let us first recapitulate a simple version of the original conjecture. Considerfor increasing L Toeplitz matrices T L of dimension L × L defined by the symbol ϕ ( k ),i.e. ( T L ) nm = (cid:90) π − π d k π e ik ( m − n ) ϕ ( k ) . (50)Provided that φ ( k ) has the following factorisation form φ ( k ) = ψ ( k ) R (cid:89) r =1 v γ r , q r ( k ) , (51)where ψ ( q ) is a continuously differentiable function and v γ r , q r describe jumps at positions k = q r in the following form v γ r , q r ( k ) = exp[ − iγ r ( π − k + q r )] , q r < k < π + q r , (52)then the L → ∞ asymptotics of the determinant isdet( T L ) = ( F [ ψ ]) L (cid:32) R (cid:89) r =1 L − γ r (cid:33) E [ ψ, { γ r } , { q r } ] , (53)where F [ ψ ] = exp (cid:16) π (cid:82) π ln ψ ( k )d k (cid:17) , and the E term does not depend on L .For translation invariant Gaussian states, the covariance matrix is usually not asimple Toeplitz matrix, rather a block-Toeplitz matrix (composed of 2 × L = (cid:32) − (cid:33) ⊗ T L (54)with T L a Toeplitz matrix, the asymptotics of det Γ L can be derived from the Fisher–Hartwig conjecture. For such models, using the contour integral representation, itcan be shown that the F [ ψ ] L factor will give rise to the linear term in the entropyasymptotics (which is zero for pure Gaussian states), while a logarithmic subleading termis induced by the (cid:81) Rr =1 L − γ r factor, and the final E [ ψ, { γ r } , { q r } ] factor only provides aconstant term [85, 86, 93, 94, 95]. Let us also mention that when calculating the mutualinformation, the linear term drops out and the logarithmic term provides the leadingorder in the asymptotics [24].The covariance matrix of the Ising NESS cannot be factorised in the form ofEq. (54), thus one has to use generalisations of the Fisher–Hartwig conjecture forobtaining the mutual information asymptotics. The Fisher–Hartwig method presented emperature driven quenches in the Ising model q ) that is continuously differentiable apart from a finitenumber of points k = q r ( r = 1 . . . R ) where it has jumps satisfying the conditionlim (cid:15) → [Λ( q r − (cid:15) ) , Λ( q r + (cid:15) )] = 0. Such 2 × k ) = U † ( k ) (cid:32) Ψ( k ) R (cid:89) r =1 V r ( k ) (cid:33) U ( k ) , (55)where Ψ( k ) is continuously differentiable diagonal 2 × U ( k ) is acontinuously differentiable function of unitary matrices that diagonalise Λ( k ), and thejump matrices V r ( k ) are of the form V r ( k ) = (cid:32) exp[ − iγ r ( π − k + q r )] 00 exp[ − iδ r ( π − k + q r )] (cid:33) , q r < k < π + q r . (56)According to the generalised Fisher–Hartwig conjecture, the determinant can again befactorised in a similar form than that of Eq. (53), i.e. det(Γ L ) = F L ( (cid:81) Rr =1 L − γ r − δ r ) E ,where F and E do not depend on L . In the block-Toeplitz case the explicit form of F and E is not known. However, since the linear term of the entropy (corresponding to F L ) drops out from the mutual information asymptotics, we are able to calculate theleading order correction of I ( α ) L .Let us turn attention to the particular case of the Ising NESS. As discussed earlier,in order to calculate the von Neumann and R´enyi entropies using the contour integral(41), we have to consider the Toeplitz matrix corresponding to the symbol λ − i Λ( k ).There are two jumps in this symbol, at k = 0 and at k = π/
2. The diagonal elementsof the jump matrices V and V are the following: γ ( λ )= 12 πi log (cid:18) λ − e − ( h +1) /T R +1) − + 1 λ − e − ( h +1) /T L +1) − + 1 (cid:19) , (57a) δ ( λ )= 12 πi log (cid:18) λ − e − ( h − /T R +1) − + 1 λ − e − ( h − /T L +1) − + 1 (cid:19) , (57b) γ ( λ )= 12 πi log (cid:18) λ − e ( h +1) /T R + 1) − + 1 λ − e ( h +1) /T L + 1) − + 1 (cid:19) , (57c) δ ( λ )= 12 πi log (cid:18) λ − e ( h − /T R + 1) − + 1 λ − e ( h − /T L + 1) − + 1 (cid:19) . (57d)Taking the logarithm of D L ( λ ) = det( λ − i Γ L ) , log D L ( λ ) = L log F ( λ ) − ( γ ( λ ) + δ ( λ ) + γ ( λ ) + δ ( λ )) log L + log E ( λ ) . (58)We will drop the log E ( λ ) term, as it only gives an L -independent value and we calculate emperature driven quenches in the Ising model I ( α ) L up to O (1) in L . Taking the derivative of log D L , we obtain:d log D L ( λ )d λ = d log( F ( λ ))d λ L + (cid:20) b − a ) πi (cid:18) γ ( λ )(2 a − − λ )(2 b − − λ ) + γ ( λ )(1 − a − λ )(1 − b − λ ) (cid:19) +2( b − a ) πi (cid:18) δ ( λ )(2 a − − λ )(2 b − − λ ) + δ ( λ )(1 − a − λ )(1 − b − λ ) (cid:19)(cid:21) log L , (59)where a , a , b , b are defined in Eq. (43).When calculating the mutual information the term proportional to L drops out,thus, using Eq. (41), we obtain that I ( α ) L = a − b π (cid:18)(cid:73) C d λ s ( α ) ( λ ) γ ( λ )(2 a − − λ )(2 b − − λ ) + (cid:73) C d λ s ( α ) ( λ ) γ ( λ )(1 − a − λ )(1 − b − λ ) (cid:19) log L + a − b π (cid:18)(cid:73) D d λ s ( α ) ( λ ) δ ( λ )(2 a − − λ )( b − − λ ) + (cid:73) D d λ s ( α ) ( λ ) δ ( λ )(1 − a − λ )(1 − b − λ ) (cid:19) log L , (60)where, due to the position of the divergences and cuts of the integration kernel, thecontour C encircling the interval [ − ,
1] could be broken up to four smaller contours C , C , D , and D which encircle the branch cuts of the four denominators. For example,contour C encircles the interval ‡ [2 a − , b − . Another simplification occurs byobserving the symmetry of the problem under the exchange of variables λ → − λ : onehas γ (1 − λ ) = − γ ( λ ) and δ (1 − λ ) = − δ ( λ ); here the negative sign cancels out withthe reversal of the directions C → −C and C → −C of the contours upon reflection.Hence two pairs of the four contributions in Eq. (60) are equal, which yields I ( α ) L = (cid:20) a − b π (cid:73) C d λ s ( α ) ( λ ) γ ( λ )(2 a − − λ )( b − λ ) + a − b π (cid:73) D d λ s ( α ) ( λ ) δ ( λ )(2 a − − λ )(2 b − − λ ) (cid:21) log L . (61)The cuts of the functions γ and δ are along the intervals (2 a − , b −
1) and(2 a − , b − γ ( x + i0 ± ) = 12 πi (cid:20) log 2 a − − x b − − x ∓ i ( π − + ) (cid:21) = γ ( x ) ∓ (cid:0) − + (cid:1) , x ∈ (2 a − , b − , (62)and similarly, δ ( x + i0 ± ) = δ ( x ) ∓ (cid:0) − + (cid:1) , x ∈ (2 a − , b − . (63)Using Eqs. (62) and (63), one can further decompose the contour integral along C and D in (61) to integrations of the jump on the intervals (2 a − (cid:15), b − − (cid:15) ) and ‡ Here and below we assume b > a but the calculation is analogous in all other cases. emperature driven quenches in the Ising model Figure 8.
The integration contour for the integrals containing a and b in Eq. (64). (2 a − (cid:15), b − − (cid:15) ) and along circular contours around the points 2 a j − b j − j = 1 , I ( α ) L = lim (cid:15) →∞ (cid:20) a − b π (cid:90) b − − (cid:15) a − (cid:15) d λ s ( α ) ( λ )(2 a − − λ )(2 b − − λ ) + a − b π (cid:90) b − − (cid:15) a − (cid:15) d λ s ( α ) ( λ )(2 a − − λ )(2 b − − λ )+ a − b π (cid:32)(cid:73) B (cid:15),a d λ s ( α ) ( λ ) γ ( λ )(2 a − − λ )(2 b − − λ ) + (cid:73) B (cid:15),b d λ s ( α ) ( λ ) γ ( λ )(2 a − − λ )(2 b − − λ ) (cid:33) + a − b π (cid:32)(cid:73) B (cid:15),a d λ s ( α ) ( λ ) δ ( λ )(2 a − − λ )(2 b − − λ ) + (cid:73) B (cid:15),b d λ s ( α ) ( λ ) δ ( λ )(2 a − − λ )(2 b − − λ ) (cid:33)(cid:35) log L , (64)where B (cid:15),v denotes a circular contour of radius (cid:15) with the point 2 v − v = a , after substituting λ = 2 a − (cid:15)e iθ ,one can evaluate this principal value integral aslim (cid:15) → (cid:73) B (cid:15),a d λ πi s ( α ) ( λ ) log( λ − a + 1) − log( λ − b + 1)( λ − a − λ − b −
1) =lim (cid:15) → (cid:90) π − π d θ π s ( α ) (2 a −
1) log(2 b − a ) − log( (cid:15) ) − iθ b − a ) =lim (cid:15) → s ( α ) (2 a − b − a ) log (cid:18) b − a ) (cid:15) (cid:19) . (65)The other principal value integrals can be calculated analogously, and we obtain theexpressionslim (cid:15) → a − b π (cid:32)(cid:73) B (cid:15),a d λ s ( α ) ( λ ) γ ( λ )(2 a − − λ )(2 b − − λ ) + (cid:73) B (cid:15),b d λ s ( α ) ( λ ) γ ( λ )(2 a − − λ )(2 b − − λ ) (cid:33) =lim (cid:15) → s ( α ) (2 a −
1) + s ( α ) (2 b − π log (cid:18) (cid:15) b − a ) (cid:19) , (66)lim (cid:15) → a − b π (cid:32)(cid:73) B (cid:15),a d λ s ( α ) ( λ ) δ ( λ )(2 a − − λ )(2 b − − λ ) + (cid:73) B (cid:15),b d λ s ( α ) ( λ ) δ ( λ )(2 a − − λ )(2 b − − λ ) (cid:33) = (67)lim (cid:15) → s ( α ) (2 a −
1) + s ( α ) (2 b − π log (cid:18) (cid:15) b − a ) (cid:19) . (68) emperature driven quenches in the Ising model I (1) L , we can evaluatethe line integrals in Eq. (64) by usinglim (cid:15) → (cid:90) b − − (cid:15) a − (cid:15) d λ s (1) ( λ )(2 a − − λ )(2 b − − λ ) = lim (cid:15) → (cid:90) b − − (cid:15) a − (cid:15) d λ − λ log λ − − λ log − λ (2 a − − λ )(2 b − − λ ) =12( a − b ) (cid:20) a Li (cid:18) a − ba (cid:19) + (1 − a )Li (cid:18) b − a − a (cid:19) + b Li (cid:18) b − ab (cid:19) + (1 − b )Li (cid:18) a − b − b (cid:19) (cid:21) + lim (cid:15) → s (1) (2 a −
1) + s (1) (2 b − b − a ) log (cid:18) (cid:15) b − a ) (cid:19) , (69)obtaining the final formula Eq. (45a).For the R´enyi entropy with integer α > (cid:15) → (cid:90) b − − (cid:15) a − (cid:15) d λ (cid:20) log( λ − z )(2 a − − λ )(2 b − − λ ) + log( λ − z )(2 a − − λ )(2 b − − λ ) (cid:21) =12( b − a ) (cid:20) π (cid:18) log (cid:18) − a − − z b − − z (cid:19) + η (cid:18) a − − z b − − z (cid:19)(cid:19)(cid:21) + lim (cid:15) → log | a − − z | + log | b − − z | b − a ) log (cid:18) (cid:15) b − a ) (cid:19) , (70)where z / ∈ R , and η ( w ) = (cid:40) πi log( w ) when arg( w ) ∈ [0 , π ) , − πi log( w ) when arg( w ) ∈ [ − π, . (71)Using the above line integral expression and Eq. (65) together with the factorisations (cid:18) λ + 12 (cid:19) + (cid:18) λ − (cid:19) = ( λ + i )( λ − i )2 , (72) (cid:18) λ + 12 (cid:19) + (cid:18) λ − (cid:19) = 3( λ + i/ √ λ − i/ √ , (73) (cid:18) λ + 12 (cid:19) + (cid:18) λ − (cid:19) = (cid:0) λ + i tan π (cid:1) (cid:0) λ + i tan π (cid:1) (cid:0) λ + i tan π (cid:1) (cid:0) λ + i tan π (cid:1) , (74) (cid:18) λ + 12 (cid:19) m + (cid:18) λ − (cid:19) m = 12 m − m − (cid:89) k =1 (cid:18) λ + i tan (2 k − πi m +1 (cid:19) , (75)we can evaluate the integral (64) for α = 2 , , , m and obtain the results stated in Eq.(45). References [1] Calabrese P and Cardy J 2009
J. Phys. A Rev. Mod. Phys. emperature driven quenches in the Ising model [3] Calabrese P, Cardy J and Doyon B 2009 J. Phys. A Rev. Mod. Phys. Phys. Rep.
Rep. Prog. Phys. Preprint arXiv:1608.00614 )[8] Wen X G 2013
ISRN Cond. Mat. Phys.
Rep. Prog. Phys. JSTAT
P08024[11] Brand˜ao F G S L and Horodecki M 2013
Nature Physics Nucl. Phys. B
Phys. Rev. Lett. JSTAT
P06002[15] Korepin V 2004
Phys. Rev. Lett. J. Math. Phys. J. Math. Phys. Proc. Natl. Acad. Sci.
Preprint arXiv:1611.04983 )[20] Groisman B, Popescu S and Winter A 2005
Phys. Rev. A JHEP Phys. Rev. Lett.
JSTAT
P02008[24] Eisler V and Zimbor´as Z 2014
Phys. Rev. A Phys. Lett. A
Phys. Rev. Lett. Phys. Rev. Lett. Int. J. Quant. Inf. Phys. Rev. Lett.
JSTAT
P11001[31] Hastings M B, Gonz´alez I, Kallin A B and Melko R G 2010
Phys. Rev. Lett.
Phys. Rev. Lett.
Phys. Rev. B Phys. Lett. B
Phys. Rev. Lett.
Phys. Rev. Lett.
Phys. Rev. Lett.
Nature
Science
Preprint arXiv:1609.02157 )[41] Melko R G, Kallin A B and Hastings M B 2010
Phys. Rev. B Phys. Rev. B Phys. Rev. Lett.
Phys. Rev. B JSTAT
P01008[46] Alcaraz F C and Rajabpour M A 2014
Phys. Rev. B Phys. Rev. D Phys. Rev. D J. Math. Phys. Comm. Math. Phys.
Comm. Math. Phys. emperature driven quenches in the Ising model [52] Gupta M K and Wilde M M 2015 Comm. Math. Phys.
Comm. Math. Phys.
J. Math. Phys. J. Math. Phys. Phys. Rev. Lett.
Phys. Rev. Lett.
Phys. Rev. E Phys. Rev. Lett. Tr. Mat. Inst. Steklova
J. Stat. Phys.
Phys. Rev. A Phys. Rev. E J. Phys. A Eur. Phys. J. B Phys. Rev. B Phys. Rev. B Phys. Rev. B Phys. Rev. Lett.
J. Phys. A Phys. Rev. B JSTAT
P08006[73] Doyon B 2012 (
Preprint arXiv:1212.1077 )[74] Castro-Alvaredo O, Chen Y, Doyon B and Hoogeveen M 2014
JSTAT
P03011[75] Castro-Alvaredo O A, Doyon B and Yoshimura T 2016
Phys. Rev. X J. Phys. A Nature Physics Annales Henri Poincar´e JSTAT
New J. Phys. Nucl. Phys. B
New J. Phys. Preprint arXiv:1704.03744 )[84] Fagotti M and Essler F H L 2013
Phys. Rev. B J. Stat. Phys.
Phys. Rev. Lett. Phys. Rev. A JSTAT
P04016[89] Eisler V, Karevski D, Platini T and Peschel I 2008
JSTAT
P01023[90] Igl´oi F, Szatmriand Z and Lin Y C 2009
Phys. Rev. B JSTAT
P08019[92] Mosonyi M, Hiai F, Ogawa T and Fannes M 2008
J. Math. Phys. Phys. Rev. A JSTAT
P08029[95] Ares F, Esteve J G, Falceto F and S´anchez-Burillo E 2014
J. Phys. A J. Phys. A Phys. Rev. A J. Phys. A Phys. Rev. A92