Entanglement scaling and spatial correlations of the transverse field Ising model with perturbations

TL;DR

This paper investigates entanglement entropy and spatial correlations in a perturbed one-dimensional transverse field Ising model, exploring non-equilibrium steady states, near-integrability effects, and periodic bond modulations.

Contribution

It provides a comprehensive numerical analysis of entanglement and correlations under various perturbations, including energy currents and periodic bond structures, expanding understanding of phase diagrams and critical behavior.

Findings

Confirmed phase diagram and entanglement entropy in steady states with energy current

Calculated correlations and entanglement properties near integrability

Identified a critical weak bond scaling with the transverse field as (h/J)^n

Abstract

We study numerically the entanglement entropy and spatial correlations of the one dimensional transverse field Ising model with three different perturbations. First, we focus on the out of equilibrium, steady state with an energy current passing through the system. By employing a variety of matrix-product state based methods, we confirm the phase diagram and compute the entanglement entropy. Second, we consider a small perturbation that takes the system away from integrability and calculate the correlations, the central charge and the entanglement entropy. Third, we consider periodically weakened bonds, exploring the phase diagram and entanglement properties first in the situation when the weak and strong bonds alternate (period two-bonds) and then the general situation of a period of n bonds. In the latter case we find a critical weak bond that scales with the transverse field as $J'_c/J$ = $(h/J)^n$, where $J$ is the strength of the strong bond, $J'$ of the weak bond and $h$ the transverse field. We explicitly show that the energy current is not a conserved quantity in this case.

Tables

Table 1. Table 1: Numerical data showing how the wavelength ℛ ℛ \cal{R} and wavenumber k 𝑘 k vary with current L 𝐿 L

$L$	$ℛ$	$k$
1.2	10.31	0.61
1.4	7.67	0.82
1.6	7.29	0.86
1.8	6.10	1.03
2.0	5.98	1.05

Equations

$$H$$

$$H_{Is}$$

$$J_{Is}^E$$

$$H=\sum _k{ϵ}_k\left({\gamma }_k^{†}{\gamma }_k-\frac{1}{2}\right)$$

$${ϵ}_k=\frac{1}{2}\left(\sqrt{J^2+h^2+2Jh\cos ka}+L\sin ka\right).$$

$$C_{zy}=\sum _n⟨{\sigma }_n^z{\sigma }_{n+1}^y-{\sigma }_n^y{\sigma }_{n+1}^z⟩.$$

$$C_{zy}/N=\frac{1}{\pi L^2}(L^2-h^2{)}^{1/2}(L^2-J^2{)}^{1/2}.$$

$$⟨{\sigma }_n^z{\sigma }_{n+R}^z⟩\sim \frac{Q(h,L)}{\sqrt{R}}.\cos (kR)$$

$$Q(h,L)\simeq Q(h,\infty ){\left(\frac{L^2-h^2}{L^2-J^2}\right)}^{1/4}$$

$$k=\arccos (1/L).$$

$$k=(1.01\pm 0.06)\arccos [(0.99\pm 0.05)/L].$$

$$S\sim \frac{c}{6}\log \left(\xi /a\right)$$

$$S=-\sum _{a_n}{\left({\lambda }_{a_n}^{[0]}\right)}^2\ln {\left({\lambda }_{a_n}^{[0]}\right)}^2=-\sum _{a_n}{\left({\lambda }_{a_n}^{[1]}\right)}^2\ln {\left({\lambda }_{a_n}^{[1]}\right)}^2$$

$$\xi \propto {\chi }^{\kappa }$$

$$H_{non}=-\sum _{i=1}^{N}J{\sigma }_i^z{\sigma }_{i+1}^z+h{\sigma }_i^x+D{\sigma }_i^x{\sigma }_{i+1}^x.$$

$${ϵ}_{h=0,k}^{non}=J{\left(1-\left[1-{\left(\frac{1+D/J}{1-D/J}\right)}^2\right]{\sin }^2(k)\right)}^{1/2}$$

$$H_{wb}=-\sum _{i=1}^N(J{\sigma }_i^z{\sigma }_{i+1}^z+h{\sigma }_i^x)-\sum _{j=1}^{[N/2]}(J^{\prime }-J){\sigma }_{2j}^z{\sigma }_{2j+1}^z$$

$${ϵ}_k^{wb}={\left[2J^2+2J^{\prime 2}+4h^2+2\sqrt{(J^2-J^{\prime 2}{)}^2+4h^4(J^2+J^{\prime 2})+8J{J}^{\prime }{h}^2\cos ka}\right]}^{1/2}$$

$$j_{2n}+j_{2n+1}=J^{\prime }h{\sigma }_{2n-1}^z{\sigma }_{2n}^y-Jh{\sigma }_{2n}^y{\sigma }_{2n+1}^z+Jh{\sigma }_{2n}^z{\sigma }_{2n+1}^y-J^{\prime }h{\sigma }_{2n+1}^y{\sigma }_{2n+2}^z$$

$${\partial }_t(j_{2n}+j_{2n+1})=-i\left[H,j_{2n}+j_{2n+1}\right]=$$

$$S\sim \frac{c}{6}\log \left(\frac{\xi }{a}\right)$$

$$\xi =1/\ln ({\lambda }_1/|{\lambda }_2|)=-1/\ln |{\lambda }_2|$$

$$\frac{\partial h_{i,i+1}(t)}{\partial t}+\Delta j_i=0$$

$$h_{i,i+1}(t)=U{h}_{i,i+1}{U}^{-1}=e^{iHt}{h}_{i,i+1}{e}^{-iHt}$$

$$\frac{\partial h_{i,i+1}(t)}{\partial t}$$

$$\frac{\partial h_{i,i+1}(t)}{\partial t}=j_i(t)-j_{i+1}(t)$$

$$j_i^E=Jh({\sigma }_{i-1}^z{\sigma }_i^y-{\sigma }_i^y{\sigma }_{i+1}^z)$$

$$j_i^E=h(J^{\prime }{\sigma }_{i-1}^z{\sigma }_i^y-J{\sigma }_i^y{\sigma }_{i+1}^z)$$

Full text

I Introduction

The transverse field Ising model (TFIM) possesses a central role in both quantum statistical and condensed matter physics.Dutta It is used as the benchmark model where new concepts, ideas and theoretical techniques have been derived from or tested. With the surge of activity on quantum information, this model also plays a very important role in simulating interactions among qubits Frontiers ; Qureshi and serves as a quantum paradigm which can be explored by novel means. While the pure TFIM can be solved exactly Pfeuty ; Lieb , relevant physical perturbations (such as longitudinal fields or spin-exchange) prohibit an exact solution. In this case advanced numerical tools have to be employed to understand the physics of the models. In this work, we study the TFIM with ferromagnetic nearest-neighbour interactions with three different perturbations. We use a combination of matrix-product state (MPS) based numerical techniques starting with the time evolving block decimation (TEBD) that allows an efficient time-evolution of matrix-product states in real or imaginary time.Vidal1 This method is similar in spirit to the density matrix renormalisation group (DMRG) method that we also employ White93 ; Schollwoeck ; McCulloch ; Kjall and has been shown to work for the entanglement spectrum near criticality in finite quantum spin chains DeChiara . Ground state entanglement for example of the XY and Heisenberg models shows the emergence of universal scaling behavior at quantum phase transitions. Entanglement is thus controlled by conformal symmetry. Away from the critical point, entanglement gets saturated by a mass scale Latorre . Here, we investigate the entanglement, as well as correlations, of the TFIM with perturbations. The simulations are performed directly in the thermodynamic limit assuming translational invariance with respect to a unit cell.

We start in Sec. II by revisiting the exactly solvable TFIM in a non-equilibrium steady state, with an energy current passing through the system. We confirm the phase diagram obtained in Refs. [Antal, ; Eisler1, ] In addition, we provide new insights for the chiral order parameter, the spin-spin correlations and the entanglement properties. In the case of central charge, we found the change from c=1/2 to c=1 as analytically calculated in Ref. [Kadar, ]. In Sec. III, we add a perturbation that breaks the integrability of the model, seeking to characterise the system using the bipartite entanglement entropy and central charge. Subsequently, in Sec. IV, we weaken the strength of a bond periodically in space, studying the properties of the system first in equilibrium. The critical value of the weak bond scales with the transverse field as $J^{\prime}_{c}/J=(h/J)^{n}$, when $n$ is the spatial period of the weakened bond. We show that the energy current is not a conserved quantity in this case, so there is no out of equilibrium steady state in the sense as in the homogeneous system. We finally summarise and conclude in Sec. V.

II TFIM with energy current

II.1 Diagonalization of Hamiltonian

The Hamiltonian for the TFIM with energy current is given by

[TABLE]

where $J$ is the strength of the nearest-neighbor interaction, $\vec{\sigma}$ are the Pauli matrices, $h$ is the transverse applied field, $J_{E}$ the energy current and $\lambda_{1}$ is a Lagrange multiplier. The energy current $J^{E}_{Is}$ is derived from the continuity equation for the conserved local energy operator Zotos (the derivation is presented in Appendix B).

This Hamiltonian is exactly solvable through a standard process by mapping it to a quadratic fermionic model. After a Bogoliubov transformation and Fourier transform to momentum space, the Hamiltonian in terms of fermionic operators $\gamma_{k}$ reads:

[TABLE]

with energy relation Antal

[TABLE]

It is convenient to use the notation $L=Jh\lambda_{1}$. For small $L\leq J$, the energy spectrum is gapped and no current flows. Then the zero-current phase is extended until $L>J$ where the imposed current density is strong enough to destroy the gap and mix the excited states with the ground state. Then we enter the region of finite energy current flow. The phase boundaries are calculated by requiring the dispersion relation and its derivative to vanish with respect to the wavenumber $k$ Antal . The values of the critical $L_{c}$ are: $L_{c}=h$ if $h\geq J$, or $J$ if $h<J$; the energy current is non-zero for values of $L\geq L_{c}$.

The transition into the current carrying region can be seen directly through the computation of the chiral order parameter, essentially the average energy current over all bonds:

[TABLE]

The expectation value of $C_{zy}$ is:

[TABLE]

At the phase boundaries, the critical behavior of $C_{zy}$ is: $C_{zy}\propto(L-h)^{1/2}$ if $h\neq J$ and $C_{zy}\propto(L-J)$ if $h=J$. This has been also verified by the numerical calculations in the present study. For example $C_{zy}$ as a function of h is shown in Fig. 1.

II.2 Spin-spin Correlations

The behavior of the spin-spin correlation function $\langle\sigma^{z}_{n}\sigma^{z}_{n+R}\rangle$ reveals distinct properties within the current region ($L=2J$ in Fig. 2). We confirm the oscillatory behavior accompanied by a power-law decay (as $1/\sqrt{R}$), where $R$ is the distance (in number of sites) between two correlated spins, as obtained in Ref.[Antal, ]. The specific power-law is considered to govern the critical region of non-equilibrium steady-state models in general Torre . The correlation function is written as:

[TABLE]

The value of $Q(h,L)$ can be calculated exactly in the limit $L\rightarrow\infty$, as the limit of the amplitude of the same correlations of the XY-model in 1D Baruch , with the result $Q(h,L=\infty)=e^{1/2}2^{-4/3}A^{-6}\simeq 0.147$, where $A\simeq 1.282$ is the Glaisher’s constant. For different values of $L$, Ref. [Antal, ] approximated $Q(h,L)$ away from the phase boundaries with:

[TABLE]

while, close to the boundary $h=J$, $Q(h=J,L)\simeq Q(h,\infty)(\frac{L^{2}}{L^{2}-J^{2}})^{1/8}$. The wavenumber $k$ of the spatial dependence of the correlations is independent of the magnetic field Antal and as $L\rightarrow\infty$ it is given by:

[TABLE]

Numerically, the correlations are computed using the iTEBD algorithm and the critical correlations in the energy current-carrying region has been verified. Using Fast Fourier Transform (FFT), the peak of the oscillations in space, determines the wavelength ${\cal{R}}=2\pi/k$. This is shown in Fig. 2).Table I lists results for current, wavelength and wavenumber. A least squares fit of the data indicates that Eq. (10) is in general reliable and we find that:

[TABLE]

II.3 Entanglement scaling and central charge

The universal entanglement properties of the TFIM are accessible within conformal field theory Calabrese . A universal formula for the entanglement entropy

[TABLE]

depends on the central charge $c$ and the correlation length $\xi$ ($a$ is the lattice spacing). In this expression for the entanglement entropy the prefactor is 1/6 (and not 1/3), due to the fact that we consider an infinite DMRG algorithm in the calculations, therefore there is only one contact point between the two parts of the spin chain, to be taken into account. This is contrasted to the case of a chain with periodic boundary condition.

To study the properties of the system, at a quantum critical point, where the correlation length is infinite, within the MPS framework, we can perform an entanglement scaling, as $\chi\rightarrow\infty$ pollmann09 ; tagliacozzo , in the spirit of the work of Ref. [White93, ]. More specifically, the entanglement entropy is calculated from the Schmidt decompositions singular values using the von Neumann entropy

[TABLE]

The values are identical for either half of the chain (odd or even bond) and the bipartite splitting of the MPS has entropy that is maximally entangled at $\log{\chi}$.

The information on the critical properties of the system can be obtained by employing finite entanglement scaling. In a critical chain the scaling equation

[TABLE]

holds, where $\chi$ is the number of states kept and $\kappa$ is a constant that depends only on the central charge of the model pollmann09 . This method to extract the central charge is followed in the next sections as well and more details are presented in Appendix A.

From the calculation of the entanglement entropy, we extract the central charge $c$ in the current-carrying region which is critical with an expected central charge $c=1$. A change in $c$ naturally occurs at the boundary to the gapped region where $c=0$. For example at fixed current $L=1.2$ and non-critical field $h=0.5J$ the system sits in the current region. The typical scaling behavior of the entanglement entropy is presented in in the inset of Fig. 3).

The system remains critical in the entire current carrying region. The transverse field dependence of the scaling is shown in Fig. 3. At a small growing perturbation in $h\ll L$ is observed a denser set of points on the left-hand side of Fig. 3, and shows a central charge $c=1$. Further increase of the magnetic field strength into the non-critical paramagnetic (PM) region hits the boundary at $h=L$ where the central charge is $c=0$. Multiple points around $L=1.2J$ show numerical noise which is due to instabilities near the transition.

III TFIM with a perturbation that breaks integrability.

Integrability breaking can be achieved by introducing certain perturbation terms in the TFIM. One way that this can be achieved in the TFIM is through the introduction of an interaction $D\sigma_{i}^{x}\sigma_{i+1}^{x}$ longitudinal to the magnetic field direction $x$ and transverse to the original spin ordering interaction $z$. The Hamiltonian reads

[TABLE]

We choose to work with this Hamiltonian, as it is to large extent unexplored compared to other ways to break the integrability of the system. If the spin couplings have equivalent interaction strength $J$=$D$=1, the Hamiltonian is the isotropic XY model with an integrability breaking $x$-directed magnetic field. Additionally, in the absence of magnetic field $h$=0, and $D\geq J$, the anisotropic XY model emerges This has a known exact solution, with energy dispersion relation

[TABLE]

As a consequence, there are two limits where the system is integrable with known critical entanglement properties.

III.1 Spin-spin correlations

The system passes through a critical point as $D$ increases, the $zz$-correlations are presented in Fig. 4. At low values of $D$ there is long-range order until the critical point is reached at $D_{c}$, correlations with power law decay emerge. The non-integrable term acts similarly to a transverse field, and at $D_{c}$, the system experiences a phase transition from a $Z_{2}$ symmetry breaking FM to PM phase.

III.2 Entanglement scaling and central charge

There are two known critical points at $h=J,D=0$ and $h=0,D=J$; the TFIM and the isotropic XY (or XX) limits, with central charges $0.5$ and $1$, respectively Vidal2 . A critical line connects these two points. The full phase diagram is found numerically using a combination of the iDMRG algorithm and finite entanglement scaling as shown in Fig. 5. As demonstrated in Appendix A, the critical line that connects the two known critical points is located by observing a peak in the effective correlation length $\xi$ as it increases with $\chi$. The entanglement scaling along the line maintains $c=0.5$ until the XY point is reached, where $c=1$. This is in complete analogy with the anisotropic XY model which only at the point of isotropy (XX model) the universality class is of the free boson with central charge c=1 but away from that point, in the critical region, the universality class is free fermions with c=1/2 Latorre2 . We see the same physics by perturbing the system differently in this case.

IV TFIM with a periodically weakened bond.

IV.1 Diagonalization of Hamiltonian

So far, the two variations of the TFIM describe homogeneous two-site interactions across the length of the spin chain. If we consider weakening the coupling at every even site, the Hamiltonian can be written as:

[TABLE]

This is equivalent to defining two Majorana fermion populations with two independent particle number operators. The Hamiltonian can be diagonalised with a four-dimensional Bogoliubov transformation Jafari with the dispersion relation:

[TABLE]

Integrating this equation over $k$ for $N\rightarrow\infty$, the exact ground state energy is calculated. The second derivative of the energy with respect to $J^{\prime}$ diverges at sufficiently weak $J^{\prime}$, Fig. 6. A second-order phase transition occurs as the system passes through this point. This becomes clear from the form of the order parameter, which changes continuously to zero upon approaching the transition bond strength. The behavior is universal and identical to the situation where the system passes through the critical point by varying the magnetic field in the TFIM.

IV.2 Spin-spin correlations

Fig. 7 shows the correlations as every even bond is weakened through the transition point at $J^{\prime}_{c}=0.25J$ when $h=0.5J$. This coincides with the peak in Fig. 6. Long-range order is maintained until $J^{\prime}$ is tuned to $J_{c}$ and power-law decay appears. On the other side of the transition $J_{c}<0.25J$ the region experiences exponential decay. When the bond strength is completely diminished $J^{\prime}=0$ and the chain is isolated into pairs. The correlations are instantaneously cut-off and the chain ceases to correlate beyond the first neighbour.

Short-range interactions are modified from the usual equivalent bond strength TFIM. The odd bond preserves a modulated pairing between spins experiencing the full strength $J$. Correlations show an almost sawtooth decay instead of a smooth curve.

IV.3 Entanglement scaling and central charge

The entanglement scaling properties are investigated at the divergent points in Fig. 6. The critical nature at $(J^{\prime}_{c}/J)=(h/J)^{2}$ shares universal properties with the TFIM. This suggests the entanglement scaling should exhibit similar behavior.

Fig. 8 presents a plot of the central charge at the critical points. The central charge is identical to the TFIM. Otherwise when $J^{\prime}/J\neq(h/J)^{2}$ the system is tuned away from criticality, the scaling disappears and the central charge is zero.

IV.4 In the presence of energy current

In the situation where an energy current is introduced in the system, in a similar way as in the homogeneous TFIM, then it can be proved that there is no energy current conservation. The system then is not in a steady state. To show this we shall take the energy current over two consecutive bonds, as schematically seen in Fig. 9 and find its time derivative. The energy current density over two sites sites $2n,2n+1$ that involves three bonds as indicated in Fig. 9 reads:

[TABLE]

The continuity equation then takes the form:

[TABLE]

The form of the above equation indicates that upon summation over all sites, the last term on the right hand side, proportional to $J^{\prime}h^{2}$, will vanish except from the boundaries but the first term on the right hand side will only vanish if $J^{\prime}=J$. Thus there is no energy current conservation in the system, contrary to the homogeneous case.

IV.5 Generalization to $n$-site period

The previous discussion considered the weakening of every second bond. The generalisation to a periodically weakened bond of a period $n$-sites reveals a relationship for the critical $J^{\prime}$ that scales with the transverse field as $J^{\prime}_{c}/J=(h/J)^{n}$ . Figure 10 (a) shows the critical lines for $n$ from 2 to 6.

The physics behind the above result, is that as the distance between the weakened bonds becomes larger, then lower values of $J^{\prime}$ are required to produce the same effect and, in the limit of very large $n$ then $J^{\prime}_{c}$ approaches [math] which indicates that it effectively requires cutting the chain in order to produce the same result. The analytical treatment and proof of this relation is beyond the scope of the present work and will be presented elsewhere as it involves the non-trivial diagonalisation of a chain of $n$ spins with periodically modulated couplings and periodic boundary conditions Feldman . The case of a single site with transverse field in an $n$-periodic chain has been investigated analytically, e.g., in Ref. [Genovese, ]; the treatment there is simplified by the fact that the defect is on the site rather than on the bond.

V Discussion

In this work we present a detailed study of entanglement entropy and critical correlations of the one-dimensional TFIM in three cases using mainly numerical calculations. We have revisited the TFIM with an energy current confirming the phase diagram and correlations Antal and the entanglement properties where we found that the central charge $c$ takes the value 1 instead of 1/2 which was also analytically calculated Kadar . Similar change of the central charge was seen in the entanglement scaling for the XX chain in the current carrying region Eisler2 . The reason is that the Fermi sea is now doubled. For symmetrically arranged multiple Fermi seas it has been shown Keating that the entanglement is proportional to the number of the Fermi seas. This has been generalised to arbitrary arrangement of the Fermi seas and to cases without particle-number conservation. It is worth mentioning that in Ref. [Kadar, ], it was pointed out that there exists a duality relation between the TI chain with current and XX chain with current.

Then as a natural continuation, we added an integrability-breaking term (an interaction of the form $\sigma^{z}\sigma^{z}$) and studied the entanglement properties. We found a critical line that connects two known critical points; the pure TFIM and the isotropic XY (or XX). Away from the point with XX symmetry, the central charge is c=1/2 along the critical line.

Subsequently by taking the TFIM and altering every other bond (from $J$ to $J^{\prime}$), we analyzed the energy dispersion relation and found the critical points. By adding a current we showed that there is no energy current conservation, which is only recovered if $J=J^{\prime}$. We generalized the model to one with $n$-site periodicity and found a critical $J^{\prime}$ that scales as $J^{\prime}_{c}/J=(h/J)^{n}$, where $n$ is the distance of two adjacent weakened bonds. This demands further analytical treatment to be reported elsewhere. It should be noted that there have been studies of isolated impurities of quantum spin chains, via bosonization Eggert and most recently using numerical techniques Biella .

Overall, taking a TFI chain as a basic model with scientific and also practical (such as in quantum information) consequences, we provide new results in different variants by using matrix-product state based methods. The entanglement properties of quantum XY spin chains of arbitrary length has been investigated by using the notion of global measure Wei1 ; Wei2 . In earlier works it was found that the field derivative of the entanglement density becomes singular along the critical line. The form of this singularity is dictated by the universality class which controls the quantum phase transition. Moreover, it was pointed out that there is a deeper connection between the global entanglement and the correlations among quantum fluctuations. This is another direction to be clarified.

Acknowledgements

We would like to thank Sasha Balanov, Claudio Castelnovo, Fabian Essler, Alex Zagoskin for useful discussions and Mike Hemsley for early collaborations on related problems. The work has been supported partly by EPSRC through grant EP/M006581/1 (JJB) and by Research Unit FOR 1807 through grants no. PO 1370/2-1 (FP).

Appendix A: Finite entanglement scaling

In infinite matrix product state methods the length of the system is no longer a limiting factor. The finite size scaling is traded for finite entanglement of the state with the bond dimension $\chi$ used as the new scaling resource, which now limits the representation of the state. The universal formula for the entanglement entropy is

[TABLE]

and depends on the correlation length $\xi$ calculated from the two dominant eigenvalues of the quantum transfer matrix

[TABLE]

As the MPS has been built in canonical form with orthonormal left and right eigenvectors, the first eigenvalue in the transfer matrix is 1, thus we need to calculate the second eigenvalue. The correlation length scales inversely with energy gap until it diverges at the critical point when the gap vanishes.

As explained in Sec. II, the entanglement entropy is calculated from the Schmidt decomposition’s singular values using the von Neumann entropy, followed by a finite entanglement scaling. In a critical regime $\xi\propto\chi^{\kappa}$. Although the perfect description of the state at a critical point requires $\chi$ to diverge to infinity, the fact that $\chi$ is finite leads to an effective finite correlation length $\xi_{\chi}$ at the critical point. The phase transition is seen at the peak in a curve as the control parameter passes through the critical point. Depending on the number of states retained the peak can be broad for low values of $\chi$ and give an less accurate representation of the exact critical point. For larger values of $\chi$ the peak becomes sharper and more reliable. This is seen in Fig. 11 which is an example of how the critical points in Section III were calculated.

Furthermore, in Fig.12 we illustrate by using examples, how the central charge is determined from the scaling of entanglement entropy.

Appendix B: Energy current derivation

Application of energy current creates a non-equilibrium steady-state if the time derivative of the total energy current is zero. In a steady-state energy flows through the system at a constant rate continuously, without any impulses causing discontinuous energy transfer. The energy is conserved and thus characterised by a conservation law. If the Hamiltonian of the system can be written as sum of terms that depend on two sites: $H=h_{1,2}+\cdots+h_{i,i+1}+\cdots+h_{N,N+1}$ and $[h_{i,i+1},h_{i-1,i}]\neq 0$ then the conservation of energy is expressed in the form of a continuity equation Zotos

[TABLE]

$\Delta j_{i}=j_{i+1}-j_{i}$ is the discrete divergence of the current.

Applying a unitary transformation (working in the Heisenberg picture)

[TABLE]

the time derivative is:

[TABLE]

where we have used the relations $dU/dt=iHU$ and $[H,U]=[H,U^{-1}]=0$. As the two commutation relations on the right hand side depend on the previous $i-1$ and next $i+1$ sites, they can be identified as the local current operators $j_{i}$ and $j_{i+1}$ respectively, meaning the continuity equation is

[TABLE]

with $j_{i}(t)=i[h_{i-1,i},h_{i,i+1}(t)]$. We now apply the above for the two cases of TFIM we have considered, finding:

1. Energy current in TFIM

[TABLE]

2. Energy current in TFIM with alternating bond strength

[TABLE]

if the strength of the interaction (bond) is $J^{\prime}$ between sites ${i-1,i}$ and $J$ between ${i,i+1}$. Note that in this case, to derive the continuity equation we need to take into account the flow of energy density passing through two adjacent sites, due to the doubling of the unit cell.