Energy transport in Z3 chiral clock model

We characterize the energy transport in a one dimensional Z3 chiral clock model. The model generalizes the Z2 symmetric Transverse Field Ising Model (TFIM). The model is parametrized by a chirality parameter θ, in addition to f and J which are analogous to the transverse field and the nearest neighbor spin coupling in the TFIM. Unlike the well-studied TFIM and XYZ models, this model does not transform to a fermionic system. We use a matrix product states implementation of the Lindblad master equation to obtain the Non-Equilibrium Steady State (NESS) in systems of sizes up to 48. We present the estimated NESS current and its scaling exponent γ as a function of θ at different f/J. The estimated γ(f/J, θ) points to ballistic energy transport along a line of integrable points f = J cos 3θ in the parameter space; all other points deviate from ballistic transport. Analysis of finite size effects within the available system sizes suggest diffusive behavior away from the integrable points.


Introduction
Though energy transport has been studied for a long time, a microscopic description of energy transport in interacting quantum and classical systems is still under development, with many recent insights on connections between chaos and transport aided by the improved simulation methods.In classical systems, chaos is neither necessary nor a sufficient condition [1] for diffusive transport.Fermi-Pasta-Ulam problem has a positive Lyapounov exponent, but does not exhibit diffusive heat conduction in any parameter regime.
An extensive amount of work on high temperature transport focusing on spinhalf models in one dimensional (1D) quantum systems [2][3][4][5][6][7][8][9][10][11] have shown that breaking integrability generally leads to diffusive energy transport.It has been analytically argued that integrability in clean systems typically leads to ballistic energy transport [7,12].Interestingly, the relation does not extend to other conserved currents [13][14][15][16][17].The XXZ chain in its zero-magnetization sector shows ballistic energy transport in all phases but spin transport is ballistic in the easy plane phase, diffusive in the easy axis phase, and super-diffusive at the isotropic point.[8,15].Both spin and energy transport are found to be ballistic in other magnetisation sectors.[7,18,19] On the other hand the same model, with a local longitudinal field, is non-integrable but shows ballistic spin transport [17].Disorder further enriches transport physics in such systems [19].
In this work, we step away from the well-studied spin-1/2 model and explore a model with a three dimensional local Hilbert space, namely the Z 3 symmetric chiral clock chain [20][21][22] which generalizes of the Z 2 symmetric TFIM [23].The latter which is mappable to free fermions is integrable and exhibits ballistic energy transport [24].The Z 3 clock model Hamiltonian is integrable in a fine tuned set of parameters but not in general.While the model shares several features with the TFIM, it is not mappable to a free fermionic Hamiltonian.We aim to address the question of how energy transport is affected by the model parameters, in particular how integrability affects transport in this model.Transport through the chain is simulated using the Lindblad master equation (LME) approach implemented using matrix product state (MPS) techniques [11,13,15,25].
Our paper is structured as follows.In Sec. 2, we describe the chiral clock model and present the details of the Lindblad dissipators.We then describe the details for the MPS implementation of the LME in Sec. 4. We find that under a change of basis, the LME and transport properties in one part of the parameter space can be related to that in another part, reducing the parameter space to be studied.This is described in Sec. 3. Results for the simulations are presented in the Sec. 5 and conclude with Sec. 6.

Model
The Z 3 chiral clock model for a chain of N spins in 1D, is described by the Hamiltonian [20][21][22]26] Each spin has a three dimensional Hilbert space, and the local operators σ and τ have the following matrix representation where ω = exp(2πι/3).We will represent the single site eigenstates of the σ operator as |1 , |ω and |ω .Operators σ and τ satisfy the algebra σ 3 i = τ 3 i = 1, σ i τ i = ωτ i σ i , and σ i τ j = τ j σ i for i = j.This algebra is a Z 3 analog of the algebra of Pauli matrices σ z and σ x .Interplay between f, θ and φ results in a rich ground state phase diagram [27][28][29] hosting trivial, topological and incommensurate phases.
The model has a global Z 3 parity symmetry associated with the operator P = Π i τ i .Apart from the global parity symmetry, the model can have other symmetries [30] namely time reversal T , charge conjugation C, and spatial inversion S depending on the values of parameters θ and φ.Under these symmetry transformations, σ and τ operators transform as T † σT = σ † , T † τ T = τ , C † σC = σ † , and C † τ C = τ † .Charge conjugation swaps the states |ω and |ω .Spatial inversion changes site index i → N − i + 1.All three symmetries are present at θ = φ = 0 while the model has only spatial inversion symmetry when θ = 0 and φ = 0. None of the three symmetries are present when both θ and φ are non-zero.In this work we will focus on the models with φ = 0 for simplicity.For φ = 0 and θ = 0, the individual symmetries C and S are broken but their products are preserved.
At f = 0, all the eigenstates of Hamiltonian can be chosen to be direct products of eigenstates of σ i .Energy of each eigenstate is −2J i cos(θ + α i ), where α i = arg( σ i / σ i+1 ) which take values from {0, ±2π/3}.When θ ∈ (−π/3, π/3), all the spins in ground state are aligned in the same direction, either in 1, ω or ω.Ground state for θ ∈ (±π/3, ±π) has consecutive spins oriented at relative angle of ±2π/3.Parameter f tunes quantum fluctuation in the model.At large f , the ordered phase is destroyed forming a paramagnetic phase.A second order phase transition separates the Z 3 symmetry broken phase (small f ) and Z 3 symmetric phase(large f ).The model was shown to be integrable along the line f = J cos 3θ inside the ordered phase [31].
There has been limited studies of transport properties in the model.Nonequilibrium current in Z 3 chiral clock chain with alternating sites are different temperatures have been studied in Ref [32].At the critical integrable point described by f = J and θ = 0, energy transport between a ground state and high energy state was studied in a generalized hydrodynamics framework in Ref. [33].We will study the energy transport in the ferromagnetic (f < J) regime and at varying values of θ.
A natural framework for investigation of transport properties is to attach baths with different characteristic temperatures at the opposite ends of the chain.This temperature difference creates an energy gradient and energy flow from high to low temperature end.In Ref. [13], Prosen et al. introduced the idea of using few-site jump operators to study transport properties under the dissipative dynamics of LME.This strategy provides computational simplicity and speedup leading to its extensive use for studying spin and fermionic chains [9,11,15,17,19,25,34,35].It has been argued that the local Lindblad approximations cannot faithfully reproduce the coherences produced by coupling to an actual quantum environment [36].The local Lindblad operators we use are intended to maintain local energy densities at the ends of the chain rather than to mimic a realistic quantum bath.We assume that the transport properties are independent of the manner in which the local energy density is realized.
Dissipative dynamics of the system with bath attached at both ends is given by the LME [37] where ρ is the density matrix of the system and D[ρ] is the Lindblad dissipator.The dissipator acts on the two sites at each end of the chain where β L and β R parametrize the inverse temperature for left and right end of the chain respectively.We define two site boundary dissipative term D i,j (β)[ρ] acting on the spin at site i and j using jump operators L a→b = |b a| as Here λ quantifies the coupling strength between the system and the bath.The twosite states |a and |b are the eigenstates, with energy eigenvalues E a and E b , of a two-site Hamiltonian h i,j acting on sites i and j.The transition rates are given by Γ ± = e ∓β(E b −Ea)/2 as shown in Fig. 1(b).h i,j contains the dominant terms of the full Hamiltonian restricted to the ends of the chain.In this manuscript, we have used two types of local boundary dissipators denoted by D 0 i,j and D θ i,j , constructed using two different choice of the two-site Hamiltonians h 0 i,j and h θ i,j .(i) D 0 i,j is defined using the two-site Hamiltonian h 0 i,j = −Jσ i σ † j + H.c. The ground state of h 0 is three fold degenerate (|11 , |ωω , and |ω ω ) and its excited state is six-fold degenerate with an energy gap of 3J between them.We have included Lindblad jump operators only between the non-degenerate eigenstates of h 0 .We note that due to the ferromagnetic nature of h 0 , use of D 0 i,j makes sense only when θ ∈ (−π/3, π/3) where the ferromagnetic states have a lower energy.(ii) The dissipator D θ i,j is constructed using A schematic representation of the jump operators in D 0 i,j and D θ i,j are shown in Fig. 1(b) and Fig. 1(c) respectively.The effective local temperatures generated by the two different dissipators as well as the length scales for thermalization near the boundary will be different for the two choice of dissipators.However we expect that qualitative features of transport will be similar in the two cases if the results are independent of the precise form of the bath.We indeed find this to be the case.
For finite dimensional systems, the LME has at least one fixed point (See Sec 4.2.2 of Ref. [38]).In small systems of upto 5 sites, we diagonalized the Liouvillian and found that it has a unique 0-eigenvalue state.Assuming the uniqueness to be true in larger systems, the time evolution under the above LME should approach a unique non-equilibrium steady state (NESS) defined as To obtain the NESS, we integrated the LME till large t and used saturation of local observables -energy current, energy density and magnetization on each site to check approach to steady state.The local energy density E θ i at site i is chosen to be the three site operator The current operator on the bond between sites i and i + 1 can be written as We evaluate this to be where The energy and current operators satisfy the discrete continuity equation The expectation value I i thermal = Tr(e −βH I i )/Z = 0 of the chosen form of the current operator is zero in the thermal state.This can be seen as follows.It can be checked that the unitary symmetry transformation operator CS introduced in Sec. 2 commutes with the Hamiltonian H(θ, 0) and anticommutes with I θ .Now consider the expectation value of the symmetry transformed current: (10) suggesting that the current as defined is zero in the the thermal state.In the first equality we have used the cyclic property of the trace and the commutation of CS with H.In the second equality, we have used the anticommutation property with I.
Fick's law can be generalized to all transport regimes using an empirical exponent γ as where κ is steady state energy conductance which scales as 1/N γ with system size N .Ballistic and diffusive transport are characterized by γ = 0 and 1 respectively.Systems exhibiting a conduction with 0 < γ < 1 and γ > 1 are said to have super-diffusive and sub-diffusive transport.We characterize the transport in the clock model from the scaling of I with N allowing us to estimate the exponent γ.

NESS currents at θ, θ + 2π/3 and −θ
In this section, we show that, under the time evolution (Eq. 3) with the dissipator D θ , the NESS current at θ is same as that at −θ and θ + 2π/3.Using this equivalence of transport behavior at different θ, we can reduce the parameter region to be studied from θ ∈ [0, 2π) to [0, π/3].To see the equivalence, we consider the unitary operators U 1 = Π i τ i i and U 2 = Π i C i .These transform the Hamiltonian as follows Transformation of the dissipator D θ [ρ] under the unitaries U 1 and U 2 is given by With these, it can be checked that the stationary solution ρ NESS to the LME (Eq. 3) at θ + 2π/3 and at −θ are related to the solution at θ by Note that we have implicitly assumed that there is only one NESS at each θ.The energy density E θ i and current I θ i transform similarly to H(θ, 0) under U 1 and U 2 .The thermal expectation value of the current at −θ is given by Similarly, we find that . These symmetries in the current and energy as a function of θ were verified in our numerical implementation of the LME.In Fig. 2 symmetry in NESS current is shown for system size N = 14 and f /J = 0.4 using the dissipator D θ .These results allow us to use the transport properties evaluated in θ ∈ [0, π/3] to infer the transport properties in the whole range [0, 2π].
Similar arguments for the case of the dissipator D 0 shows that the energy and currents at θ and −θ are equal to each other.

Numerical Implementation
Evolution under the LME (Eq. 3) was implemented using the Matrix Product State (MPS) formalism where we represent ρ as an MPS of the form Each tensor A has physical indices of dimension 9.The MPS is normalized such that the density matrix satisfies the trace preserving condition: In the LME (Eq. 3) operators can act on the density matrix ρ either from the left or right.Equivalent matrix product operator for the right and left action of operators on ρ contracts with non primed and primed indices respectively.We can write Eq. 3 in the super-operator form where L is a time independent super-operator called Liouvillian operator.The solution to Eq. 18 which can be formally written as |ρ(t) = e Lt |ρ(0) can be evaluated using a fourth-order approximant to MPO similar to those used in Refs.[39,40].Matrix exponential approximant of any order can be expressed as product of several first order approximants W II (τ ) = I + τ L as where τ 's are complex numbers proportional to t.To obtain an approximant correct till order p, we match coefficients of t of each order up to p on both sides of Eq. 19.The τ i are solutions of these p simultaneous nonlinear equations.Assuming that the fixed point is unique, the choice of initial state should not affect the NESS.For completeness we describe the initial state preparation.We started with an infinite temperature state and time evolved it under the following Liouvillian L with Lindblad dissipators D i,i+1 (Eq.5) connected to all sites with inverse temperature at each site linearly varying with site number between β L and β R .The steady state of the time evolution under L (t) is later used as an initial state for the actual time evolution.The initial state as well as the time evolved states are in equal mixtures of the three Z 3 parity quantum numbers.
The inverse temperatures at the left and right ends are β L = 0.133 and β R = 0.266 respectively.The spin coupling is set to be J = 1 and the coupling to the Lindblad dissipators is set to λ = 0.05.Simulations were performed for systems with f = 0.4 and for a set of θ in the range [0, π/3].Calculations were separately performed using the two different Lindblad dissipators D 0 and D θ .In all of our calculations, bond dimension χ being used is 200.For a select set of parameters we increased the bond dimension to 800, and no significant change was observed beyond 200 in the local observables.In panel (c) we compare the estimated γ obtained using the same dissipator D 0 but for f /J = 0.2 and f /J = 0.4.

Results
In this section we report the main results of the numerical simulations.The estimated current and energy density in the NESS, and the scaling exponent γ of the current as a function of system size are presented.In addition, we also present the level spacing statistics and the operator space entanglement entropy in the NESS.

NESS Current and Conductance
The mean NESS energy current  chain is used instead) Panel (a) shows the NESS current obtained using the dissipator D θ for parameter f /J = 0.4.Panels (b) and (c) show the same for the dissipator D 0 for model parameter f /J = 0.4 and 0.2 respectively.In all cases we find a peak current at the θ where we expect the system to be integrable.When the model parameters are changed from f /J = 0.4 to f /J = 0.2, the θ at which the model is integrable changes.Accordingly the location of the peak current also changes.The NESS current is independent of the system size at the integrable point, consistent with it exhibiting a ballistic transport.The current decreases with the system size at other θ.These qualitative features are the same for both choice of dissipators.
At each value of θ, the system size dependence of the NESS current can be parametrized using γ obtained by fitting the NESS current measured in different system of sizes from N = 14 to N = 48 to the form . N eff is the effective length of the chain which is N − 4 as two spins from each end is associated with the Lindblad dissipators.Figure 4(a) shows the current as a function of system size for a representative set of values of θ.Within the range of system sizes accessible, we are able to fit the data to the power law form.
We present the estimated γ(θ) as a function of θ in panels (b) and (c) of Fig. 4. The estimates suggest a clear ballistic energy transport only at the integrable point where γ ≈ 0. In Panel (b) of Fig. 4, scaling exponents computed using the two different dissipators show qualitatively the same behavior, and the two estimates quantitatively agree except in a region near small θ.We suspect that the difference at the small θ may be a consequence of different length scales associated with thermalization at the boundary, resulting in different effective lengths for the chain.In Fig. 4(c) γ is plotted for NESS obtained using the dissipator D 0 for f /J = 0.2 and 0.4, showing ballistic transport at the expected value of θ = 1 3 cos −1 (f /J).Level spacing statistics (within a symmetry sector of Z 3 parity) computed in a finite system of size N = 11 (Fig. 5) show Poisson statistics at the integrable point and a mixture of GOE and Poisson distributions at other values of θ.The distribution is closer to GOE away from the integrable points.Consistent with this, the estimates of γ increase away from the integrable points, however it does not indicate fully diffusive behavior in any region of θ.Studies in disordered spin-1/2 systems have suggested large length scales at weak disorder leading to super-diffusive behavior being observed in finite size calculations [8,9].We cannot rule out a similar possibility -that a diffusive behavior emerges in larger systems -with the results from the currently accessible system sizes.Spatial profiles of the energy density and current in the NESS for the super-diffusive and ballistic cases are shown in Fig. 6.As expected the energy density is independent of the position in the bulk in the case of the ballistic system.The analysis in this section relies on the scaling of the current with system size.This yields γ provided that the energy densities at the ends of the chain are independent of the system sizes (such that conductance is proportional to the current).In very large systems this can be true, but in small systems the energy densities can be affected by the bath at the other end, resulting in an energy difference that is system size dependent.An estimate of the local energy density that will be realized at the ends if there were local equilibration near the bath can be obtained by attaching only bath to the system.We performed this calculation for each of the two baths.Figure 7 presents the results one of these calculations.Figure 8 shows examples of energy densities as a function of position for different system sizes and parameter regimes (sites very close to the baths have been excluded).The estimates of the expected energy densities if the baths had locally equilibrated with the ends of the chain are shown in dotted lines.At the θ very close to the integrable point (Fig. 8 panel (b)), the energy densities are midway between the bath energy densities (dotted lines).The energies are approximately independent of the position and system size.In the case of the θ larger than the integrable value (panel (d) of Fig 8), the energy densities realized in the chain are very close to the bath energy densities (indicated by the dotted lines).In the case of θ smaller than that of the integrable point, the energy densities are position dependent but are far from the estimated bath energy densities.The system size dependence of these energy density difference may then need to be taken into account to make a correct estimate of γ.
In Fig. 9 we show the results of the γ estimated from the scaling with system size of the conductance.In order to define the conductance, we have assumed that the energy density differences are proportional to temperature differences, taking the ratio of the current to the energy density difference between the 4th site from either ends of the chain, distance between them being N − 7. The scaling exponent obtained by fitting the conductance to N −γ(θ) in the panels (c) and (d).The results indicate a larger value of γ than what was obtained from scaling of current.
For θ larger than that of the integrable point, the sites near the ends appear to have nearly equilibrated with the bath (Fig. 8(d)).In these cases we find the scaling γ to be very close to that of a diffusive system.For smaller θ, where the energy gradients are smaller and much larger system sizes may be needed in order to reliably estimate the true scaling properties.We have not shown the conductance scaling in the vicinity of the integrable points as the energy gradients are nearly zero and numerically estimated conductances show wild variations.
We now discuss a broader range of f values.For not too small system sizes, we expect the peak current and conductance κ to occur at the θ values exhibiting ballistic transport.We may therefore use the peak conductance at each f as a proxy to identify the values of θ at each f exhibiting ballistic transport.Figure 10 shows the estimated current re-scaled and shifted by f -dependent constants chosen such that for each f , the maximum value of I rescaled is 1 and minimum is 0. Within the numerical uncertainties due to the finite resolution of θ values, we find that the peak current occurs along the expected line f /J = cos(3θ) of integrable points [26,31].

Operator space entanglement
Analogous to the notion of entanglement between different bipartitions of many body states, one can define an operator space entanglement entropy (OSEE) [41,42] from the MPS representation of the density operator.From the Schmidt decomposition of the state across a partition located at bond i, the entropy can be computed as   Empirically we find that at the integrable points, away from the edges, the OSEE is independent of the location of the partition, and for the non-integrable points, S i shows weak position dependence.The singular values from which the OSEE was constructed also is weakly position dependent in the case of the integrable points.Translation invariance of the entropy as well as of the expectation values of the local operators -energy density and current -at the integrable point suggest the possibility of a translation invariant MPS approximation for the NESS at the integrable points similar to Ref [43].

Conclusion
A large body of studies on quantum transport in spin chains performed primarily on spin-half models have indicated that integrable systems show a ballistic energy transport and deviations from integrability generally lead to a diffusive behavior [3,5,6] with possible exceptions [17].
In this work we have studied the transport properties of the Z 3 clock model that goes beyond the spin half chains.At the integrable points in the model parameter space, NESS shows a system size independent current, suggesting a ballistic energy transport.At all other values of the parameters the current decreases with the system size.The transport scaling exponent γ estimated from scaling of the current alone shows indicates a super-diffusive behavior.Careful analysis of the energy density profiles suggests that this is likely to be a consequence of finite size effects in the system.System size dependence of the energy gradient also needs to be taken into account.The scaling exponent inferred from the conductance instead shows the values closer to diffusive behavior.The results demonstrate the connection between integrability and ballistic transport in a larger class of models beyond the well-studied spin half chains.
We have used local Lindblad coupling to the edges of a finite chain of chiral Z 3 clock to approximately model the coupling of the system to the bath.Within this approach, we obtained similar results when different dissipator models were used at the edge, suggesting a robustness of the results to the precise nature of the coupling of the system to the bath.Direct computation of the Drude weights can be an independent approach to verify the characterization of transport properties in the model [44][45][46][47][48].

Figure 1 :
Figure 1: Pictorial representation of the bath system setup is shown in panel (a).Action of the jump operators L |11 →|1ω and L |1ω →|11 are shown for the ferromagnetic regime for dissipators D 0 and D θ is shown in (b) and (c).

D θ f /J = 0. 40 Figure 2 :
Figure 2: Plot of the NESS current I θ as a function of θ for its full range of values from 0 to 2π showing the equivalence of transport properties at θ, −θ and θ + 2π/3.Data is shown for the system size N = 14, with the dissipator D θ at model parameters f /J = 0.4.Vertical lines show multiples of π/3.

Figure 3 :
Figure 3: NESS current I θ as a function of θ in the ferromagnetic regime.Panel (a) shows the current when the Liouvillian is defined using the dissipator D θ and f /J = 0.4.Panel (b) and (c) show the current for the case of the dissipator D 0 with f /J = 0.4 and f /J = 0.2 respectively.Different lines indicate different system sizes.The peak current appears at the integrable point θ = cos −1 (f /J)/3 shown by vertical dashed lines in all cases.

Figure 4 :
Figure 4: Panel (a) presents I θ (N ) (rescaled by the current in the smallest system size) vs N − 4 in log-scale for fixed values of θ for f /J = 0.4 and using the dissipator D θ .Panels (b) and (c) show the scaling exponent γ estimated from the system size dependence of the NESS current I θ (N ).Estimated γ is plotted as function of θ in panels (b) and (c).In panel (b), we compare the exponent γ obtained from the two different choice of dissipators D 0 and D θ .In panel (c) we compare the estimated γ obtained using the same dissipator D 0 but for f /J = 0.2 and f /J = 0.4.
as a function of the chiral parameter θ is shown in Fig 3 (results do not change if the current at the center of the

Figure 5 :
Figure 5: Level spacing distributions for θ = 0.1 and θ = 0.456 (close to the ballistic point) are shown in panels (a) and (b) respectively.In (c), variation of mean level spacing ratio r is plotted versus θ showing change in level spacing statistics from GOE to Poissonian at the integrable point cos −1 (f /J)/3.

Figure 6 :
Figure 6: Spatial profile of energy density and current for system sizes N =14, 20, 24, and 28 with position shown on the x-axis rescaled by a factor of 1/(N − 2).Profiles for Hamiltonian parameters θ = 0.2 in panel (a) and (d), θ = 0.38 in panel (b) and (e), and θ = 0.2 in panel (c) and (f).θ ∼ 0.38 is close to the integrable point.

Figure 7 :
Figure 7: Energy density as a function of the position in the NESS obtained after attaching only one bath to a chain.The two different lines indicate the energy densities realized upon attaching baths with parameters β R and β L .Different overlapping lines of different thicknesses show the data for different system sizes.

Figure 8 :
Figure 8: Each panel shows the energy densities as a function of the position for different system sizes.Position on the x-axis has been rescaled and shifted that center of the chain is at 0 and the 4 th spin from the ends are at ±0.5.The two dotted lines show the expected energy densities had the each one of the baths fully equilibrated with the chain (See Fig 7).The four panels show the data for four different cases.Panels (a) and (c) show results for θ less than that of the integrable point.Panel (b) shows the data at a θ very close to the integrable point.Panel (d) shows the same at θ larger than that of the integrable point.

Figure 9 :
Figure 9: In panel (a) and (b), log(I θ (N )/∆E) vs log(N − 7) is shown for both dissipators D 0 and D θ .Scaling exponent γ is obtained by linearly fitting log(I θ (N )/∆E) vs log(N − 7) data and is plotted as function of θ in panels (c) and (d).θ in the vicinity of integrable points (vertical dashed line) are not shown as the numerically obtained conductance κ show wild oscillations due to vanishing energy gradient.

Figure 10 :
Figure 10: (a) Rescaled current I rescaled of NESS as a function of θ and f /J is plotted.Comparison of numerically estimated θ ballistic and the integrable line θ = cos −1 (f /J)/3 is shown in (b).

Figure 11 :
Figure 11: Operator space entanglement entropy S i plotted as a function of the bond location i.All data are for system size N = 32 and at f /J = 0.4.The two panels show the entropy for the NESS obtained under the dissipators D 0 and D θ .