A Level spacing statistics

# Energy transport between two integrable spin chains

## Abstract

We study the energy transport in a system of two half-infinite XXZ chains initially kept separated at different temperatures, and later connected and let free to evolve unitarily. By changing independently the parameters of the two halves, we highlight, through bosonisation and time-dependent matrix-product-state simulations, the different contributions of low-lying bosonic modes and of fermionic quasi-particles to the energy transport. In the simulations we also observe that the energy current reaches a finite value which only slowly decays to zero. The general pictures that emerges is the following. Since integrability is only locally broken in this model, a pre-equilibration behaviour may appear. In particular, when the sound velocities of the bosonic modes of the two halves match, the low-temperature energy current is almost stationary and described by a formula with a non-universal prefactor interpreted as a transmission coefficient. Thermalisation, characterized by the absence of any energy flow, occurs only on longer time-scales which are not accessible with our numerics.

## I Introduction

Energy transport in one-dimensional systems is a fundamental problem in non-equilibrium physics. Its understanding, besides the clear importance for technological applications (1); (2), is tightly bound to a number of fundamental issues in statistical mechanics. The existence of an anomalous transport regime in opposition to the failure of Fourier’s law, and its connection to the underlying chaotic dynamics were extensively studied in classical interacting systems (3); (4). In quantum systems, the long standing interest in heat transport in nanowires and one-dimensional edge modes (see e.g. (5); (6); (7)) has been recently reinforced by corresponding studies in cold atomic systems (8).

In the spirit of quantum quenches (9), which are typically implemented in cold-atom setups (10); (11); (12), energy transport in one-dimensional quantum many-body systems can be addressed by means of the following partitioning protocol (13); (14); (15); (16). Two (ideally semi-infinite) chains are prepared in thermal equilibrium with different temperatures (see the sketch in Fig. 1). The two chains are then connected and the system evolves unitarily to a steady state which depends on the initial state and on the total Hamiltonian of the two half-chains and the connecting link. A transient regime first appears in which an energy current develops between the two halves. Then, in the thermodynamic limit, a steady state will be reached and this long-time behaviour may be characterized by a finite or zero energy current. Notice that this situation is quite different from the setups usually employed for solid-state devices, where a finite-length wire is connected at both ends to two thermal reservoirs.

In the case of homogeneous interacting chains approximated at low-energy by a conformal field theory (CFT), the low-temperature energy current was predicted to be finite and universal in the steady state (15); (16). The result was further confirmed numerically in spin chains (17) and analytically for free fermions (18). While it is clear that this description must apply in some time regimes, in the stationary state the non-universal features of the specific model under consideration might play a major role. In particular, it is not yet clear whether the energy current can be different from zero in the non-equilibrium steady state (NESS) of general interacting spin-chain models.

Earlier works addressing transport in one-dimensional interacting quantum systems within the framework of linear-response theory have pointed out that integrability and in particular the existence of conserved quantities is crucial for the presence of a persistent energy current (19). There, the persistence of the current was related to an equilibrium Drude weight; afterwards the relation between integrability and transport properties fully emerged and it was intensively scrutinized in numerous important works (see (20); (21); (22); (23); (24) and Refs. therein). These results call for the understanding of the intimate relation, if any, between the properties of the non-equilibrium steady state (NESS) and the universal features of the evolution, such as integrability, also in the partitioning protocol.

In an early work, it was suggested that the crucial property is not mere integrability but the presence of a conserved quantity which has non-zero overlap with the energy current (17). Following the ideas presented in Ref. (25), an analytical ansatz for the steady state of the spin-1/2 XXZ model, where the total energy current commutes with the Hamiltonian, was proposed and found in good agreement with numerical simulations (26). More recently, under reasonable hypotheses it was demonstrated that when the energy current is a conserved quantity, the energy flow must be ballistic and a non-vanishing stationary current is expected (27). For a review of several other peculiar features of the partitioning protocol, see Refs. (28); (29); (30); (31); (33); (34).

In its simplicity, the quench protocol considered here poses a number of very interesting questions, some of them still to be answered. What are the conditions for the steady-state current to be finite? On the other hand, when is the overall system going to equilibrate to a thermal state? Moreover, in one-dimensional spin-1/2 chains bosonization highlights the differences between low-lying bosonic modes and the high-energy fermionic quasi-particles. How do the different scattering properties of these excitations affect the universal low-temperature results for the energy current derived in Ref. (15); (16)? How do these results extend to higher temperatures?

In order to contribute to some of these issues, we consider the energy transport between two half-infinite XXZ chains with different (and possibly space-dependent) couplings. In this setup, we can explore how, at low energies, different interaction parameters and different velocities of propagation for the bosonic modes affect the value of the steady-state current. Moreover, in this way we are able to study the effect of local integrability-breaking on the dynamics of the energy current. In most of the paper we will consider the case in which the couplings are uniform in each half-chain (but with a difference between left and right). We also consider a smooth crossover between the two halves which extends over a finite region. In this way we can also investigate the impact of an adiabatic junction on the scattering of the low-lying modes and hence on the energy transport.

What we found is that a pre-equilibration behaviour may set in: the energy current reaches a finite value and then slowly decays toward zero. At low-temperature, this can be justified analytically and is peculiar of the setup we consider. In particular, when the sound velocities of the two halves are matched, the energy current becomes almost stationary and its value is not universal. Such pre-thermal regime is still observable at intermediate temperatures by means of numerical simulations. The asymptotic long-time regime is characterized by the absence of any energy flow. The onset of an effective temperature occurs on longer time-scales not observable with our methods.

The paper is organised as follows. In the next Section we define in detail the model that we consider, the partitioning protocol anticipated in this introduction, and the observables under exam. In Section III we discuss the low-temperature regime, where bosonisation leads to an analytic form of the steady-state current. Numerical results confirming the previous predictions (but also addressing higher temperatures) are presented in Section IV. These calculations are performed with algorithms exploiting the Matrix Product States (MPS) formalism. On the basis of these results, in Section V we present a general scenario for energy transport in inhomogeneous quasi-integrable models. The paper ends with our conclusions.

## Ii Transport between two XXZ spin-1/2 chains with different parameters

Let us begin by introducing the XXZ model and the partitioning protocol. We consider the Hamiltonian (for an introduction to the properties of this model see Refs. (35); (36)):

 ^H=^Hl+^Hr+^hl,0,^Hl=∑n<0^hl,n,^Hr=∑n>0^hr,n, ^hλ,n=Jλ[^Sxn^Sxn+1+^Syn^Syn+1+Δλ^Szn^Szn+1], (1)

being the -th component of the spin- on site () and . We restrict ourselves to the critical phase (as this is most relevant for the energy transport) and assume . At the beginning the system is separated into two independent halves held at different inverse temperatures

 ^ρ0=Z−1e−βl^Hl⊗e−βr^Hr, (2)

where ensures the normalization of the density matrix. For times the state is unitarily evolved with the Hamiltonian defined in Eq. (1), so that the initially-separated reservoirs are put in contact. A heat flow is generated, due to the temperature unbalance in the initial conditions.

 J= limt→∞J(t)=limt→∞Tr[^ρ(t)^j0], (3) ^j0= i2[^H,^Hr−^Hl]. (4)

In order to point out the existence of a stationary energy current in one-dimensional systems, it is customary to employ arguments based on the linear-response theory. These arguments are not applicable here, so that we do not expect any ballistic transport of energy at large times, . In particular, we provide evidence that the system is non-integrable and that no conserved quantity prevents the current from decaying to zero.

We first observe that and describe two XXZ models with different parameters and . They are integrable, as confirmed by the Poissonian level-spacing statistics and the existence of an infinite number of conserved local operators. In particular, the energy current operator

 ^jλ,n=i[^hλ,n−1,^hλ,n] (5)

is a conserved density because of the following continuity equation:

 ∂t^jλ,n=i[^Hλ,^jλ,n]=^kλ,n−1−^kλ,n. (6)

Thus, the integrated energy current from site to site changes in time only because of the boundary terms and . Operators are usually interpreted as longitudinal pressures (24). Additionally, the two Hamiltonians are exactly solvable with Bethe ansatz even in presence of a hard-wall boundary, as in this case. Even if the Hamiltonian in Eq. (1) is obtained by joining two integrable models, whenever or (with or ) the integrability is broken. This can be directly observed by inspection of the level-spacing statistics, which already at finite size suggests an abrupt change to a Wigner-Dyson distribution (see Appendix A for more details). We note that defining integrability is a subtle problem in quantum mechanics (37). Here, we refer to the notion of integrability associated to Bethe ansatz solvability, for which the level spacing statistics has been shown to be Poissonian (9); (38). At the same time, this model is not known to display any non-trivial local conserved quantity beyond the total energy and magnetization. Thus, there seems to be no element in the dynamics that supports a stationary current (19); (9); (28); (39).

However, the original properties of the two halves translate into just a local source of integrability breaking. Let us consider for example the energy-current operator:

 ^ja,b=−1∑n=a^jl,n+2^j0+b∑n=2^jr,n,a<0,b>0. (7)

The cases where or can be easily deduced but are not considered explicitly here. The time derivative of reads:

 i[^H,^ja,b]=^kl,a−^kr,b+^Θ (8)

where is a local operator with support on a few sites around (see Appendix B for the explicit formula).

When , the fact that and are localized on a few sites which are far from the junction allows for the derivation of a lower bound on which is valid beyond linear response theory. This has been explicitly carried out in a general situation which includes the case when and  (27). The argument links the local energy current to for assuming a sufficiently regular transient behaviour of along the chain. Even if this argument cannot be applied to our case of interest, the locality of the operator may allow the persistence of a quasi-stationary current in the long-time regime. In the following we will see that a pre-equilibration behavior emerges with non-zero energy current, which in some regimes can be quantitatively predicted. Thermalisation and in particular the absence of energy flows () occur only at later stages which are not numerically accessible. We will characterize them by general arguments. This does not happen in the special case and , which has been the subject of extensive numerical and analytical studies (29); (26). In this case there is no occurrence of thermalisation even at very large times, and, in particular, the stationary state depends strongly on the initial condition in the two halves.

## Iii Low-energy limit

We begin by presenting an analytical formula for in the pre-equilibration limit. In the low-temperature regime, one can exploit the low-energy (LE) field-theory and bosonisation to derive a transmission coefficient which originates from the inhomogeneity of the Hamiltonian in Eq. (1) and generalises the formula obtained for the homogeneous case in Refs. (15); (16).

We start rewriting the Hamiltonian density in Eq. (1) in terms of canonical local fermionic operators by standard Jordan-Wigner transformation as

 ^hλ,n=−Jλ2(^c†n+1^cn+h.c.)++JλΔλ(^c†n^cn−12)(^c†n+1^cn+1−12)\leavevmode\nobreak . (9)

Since the inhomogeneous model preserves the total magnetization and for its ground state is still in the zero magnetization sector (), we have . Thus we can rewrite the Hamiltonian (9) as

 ^hλ,n=−Jλ2(^c†n+1^cn+h.c.)+JλΔλ:^c†n^cn::^c†n+1^cn+1:\leavevmode\nobreak . (10)

where the notation above stands for the normal ordering prescription . In order to get a low-energy theory we expand the local fermionic operators in terms of continuous fields. Starting from the non-interacting theory, we describe the excitations near the Fermi momentum , where is the lattice spacing, . Then, it is possible to expand the operator in the following way

 ^cn√a≃eikFx^ψ+(x)+e−ikFx^ψ−(x), (11)

where are the chiral Dirac fields. This procedure is described in details in Ref. (40) for general non-translational invariant Hamiltonians of the type (1) allowing for site-dependent parameters.

Subsequently, the fermionic fields are bosonised as , with bosonic action , where . The effective low-energy theory is in our case specified by the Hamiltonian . The first piece is the inhomogeneous Luttinger liquid (LL) Hamiltonian

 ^HLL=12∫∞−∞dxu(x)\leavevmode\nobreak [(∂x^ϕ)2K(x)+K(x)(∂x^θ)2], (12)

with fields obeying and . The parameters and are the Luttinger parameter and the renormalised Fermi velocity. For an abrupt junction for and for , analogously for negative values of and for . Those parameters are related to the lattice coupling constants in (1) as follows

 Kλ=12[1−1πarccosΔλ]−1, (13)

and with .

The perturbation is the localised backscattering operator, responsible for reflection of fermionic waves through the junction at , and is given by (41); (40)

 ^V=λ(ei√4π^ϕ(x=0)+h.c.). (14)

For an abrupt junction (41); (40), the coupling is proportional to the difference between the renormalized Fermi velocities on the two sides of the chain, i.e. . Its relevance at the fixed point described by the inhomogeneous LL Hamiltonian (12) depends on the parameter

 K=2(1Kl+1Kr)−1, (15)

in close analogy with the seminal work (5). It is relevant when , marginal at , and irrelevant for . Equating the two renormalised Fermi velocities we can directly set the coupling and explore the inhomogeneous LL fixed point. The other fixed point of the renormalization group (RG) flow is characterised by a formally infinite backscattering of the fermionic modes and should be insulating, see Fig. 2. The numerical analysis of Sec. IV shows in this case a slow decay of the energy current.

The complete study of the flow between these two fixed points at finite temperatures with techniques analogous to those in (42) is problematic; in the following we limit ourselves to the derivation of the explicit form of the energy current at the inhomogeneous LL fixed point, relying on the techniques developed in (31). We will also not make any attempt to compute higher cumulants of the energy current as discussed in a similar field theoretical framework for charge transport in (32).

Determining the current requires to study the long-time dynamics of two different conformal field theories (free bosons with different compactification radii) at different temperatures when the interface, or defect, preserves conformal symmetry. The Appendix C contains the technical details of the calculation for the interested reader, the result for the energy current is

 JLE=πT12(β−2l−β−2r), (16)

where the transmission coefficient is

 T=4α(1+α)2,α=KlKr. (17)

The coefficient has a simple physical interpretation: It is the transmission coefficient associated to the scattering of the bosonic low lying modes past the step in the coupling constants at the site connecting the two half-chains. It can be obtained by solving the classical equation of motion for the bosonic field (see the sketch in Fig. 3). Such an equation can be obtained from the Lagrangian that follows from (12) having set and reads

 ∂2tϕ(x,t)=K(x)∂x[∂xϕ(x,t)K(x)] (18)

with boundary conditions and

 ∂xϕ(0−,t)Kl=∂xϕ(0+,t)Kr. (19)

The above condition ensures continuity of the momentum density at . We can look for a plane-wave solution of the equation of motion (18) and (19) of the form for and for and solve for the reflection and transmission coefficients and with . The transmission coefficient which is obtained coincides with (17).

Even if the expression for the current is corrected by a transmission coefficient, it still retains the simple functional form

 JLE=f(βl)−f(βr), (20)

firstly observed numerically in Ref. (17). The CFT analysis assumes a spectrum with linear dispersion relation where all the particles travel at the same speed. The curvature of the lattice dispersion, or from a field theory perspective irrelevant perturbations, may alter especially at higher temperatures the intuitive picture of a homogeneous NESS inside the light-cone spreading out ballistically. The major consequence will be the emergence of a truly inhomogeneous current profile (33); (34) possibly describable by a CFT in a curved space-time (43). Conscientiously the prediction (16) holds at fixed spacial coordinate after having waited a time such that .

The mechanism for thermalisation at larger times cannot be captured by the low-energy theory that approximates the original non-integrable Hamiltonian (Eq. (1)) with a free bosonic field theory. Nevertheless, we expect (16) to accurately describe the pre-equilibration behavior near the junction which, as we will see, emerges at times . The robustness of ballistic energy transport is numerically confirmed at short times, see Sec. IV, also in the presence of a relevant backscattering perturbation. Finally notice that for a smooth junction, the anisotropy parameter is constant over a length scale of one single lattice site. Bosonization can follow the same procedure as in the homogenenous case where the backscattering operator (14) is ruled-out precisely by such a symmetry, see for example (44).

## Iv Numerical analysis

In order to complement the results detailed in the previous section and extend them outside the low-temperature limit, we analyze the problem numerically. The algorithms we use are detailed in Appendix D and employ a MPS representation of the purified thermal and time-evolved state (45); (46). By means of a time-evolving block-decimation (TEBD) procedure (47); (48); (49), we are able to compute the time-evolution of the thermal system (23); (50) and thus of the current flowing at the junction.

In Sec. IV.1 we test the predictions of the low-energy effective theory developed in the previous Section, and in particular the validity of the formula (16) when . The more general case is considered in Sec. IV.2, where we also discuss differences and analogies with the previous case.

### iv.1 Fermi-velocity matching case (ul=ur)

In this subsection we want to validate the predictions of the inhomogeneous Luttinger liquid under the Fermi velocity matching condition in Eq. (17). The fact that the factorization of the energy current approximately holds also in the inhomogeneous case has been numerically verified in Ref. (17) in a wide range of temperatures. This allows us to verify Eq. (17) preparing just one of the two halves at low temperature (in our case the right one) irrespectively of the temperature of the other half, a trick which allows us to explore longer times.

We simulate a free-fermion chain () coupled to an interacting one (). The hopping parameter of the right half is chosen in order to fulfill . The left chain is initially prepared at high temperature while the right chain is initialized at low temperature . In Fig. 4 we plot the theoretical prediction of the transmission coefficient . Unfortunately, the study of systems where the transmission coefficient is significantly different from , and thus the presence of a transmission coefficient is most clearly verified, proved to be beyond our numerical possibilities. In particular, in the cases where (see Fig. 4), the time-evolution is significantly slower (this is due to the fact that an increasing mismatch of the microscopic parameters of the two halves makes the state more complex at the junction) and does not allow a reliable assessment of the quasi-stationary behaviour because of a longer transient dynamics.

In Fig. 5 we present our numerical results. In all the cases we observe that after a transient time the energy current displays a quasi-stationary value, reminiscent of the physics of integrable homogeneous models. However it should be stressed that, for the accessible times, the energy current has a residual time dependence (more evident in the cases) which makes the extrapolation of its quasi-stationary value less accurate. For this reason, the quasi-stationary value of the current is obtained after averaging over times (dashed vertical line in the top panels of Fig. 5). The value is then plotted as a function of . The error bars quantify the deviations from the average value, they are evaluated considering the maximum and minimum values of the energy current for . Our numerical data confirm that at low temperatures the energy current scales with and the slope is compatible with the theoretical prediction of the transmission coefficient in Eq. (17).

The times that we are able to access for this setup are limited and it is very hard to make definite statements about the stationary values of the current. While the time increases, the correlations spread inside a light cone determined by the Fermi velocities. From the numerical point of view, this implies that the bond link needed to faithfully represent the state of the system increases exponentially in time (see Appendix D). For this reason, our numerics is not effective in exploring the time-scale of thermalisation for which the energy current vanishes.

### iv.2 Generic case (ul≠ur)

In this subsection we consider the generic case where the parameters of the two chains are chosen such that . As discussed in Sec. III, the low-energy analysis tells us that in this case the backscattering operator performs a non-trivial renormalization flow, pointing out the prominent importance of the parameter in Eq. (15). Our simulations show that at low temperatures the behaviour of the energy current is qualitatively different from the case.

In Fig. 6 we present our results for three cases characterized by , and . The parameters are chosen so that the Fermi velocities of the simulations are of order , so that results are comparable for the considered time scales. Even if this numerical analysis is not conclusive (because the longest accessible times are not enough to obtain a full picture of the pre-thermalisation and thermalisation dynamics), it is still possible to make some interesting observations. Looking at in Fig. 6 (top) and at its time derivative in Fig. 6 (bottom), we can identify two different behaviours. The former characterizes the cases and and displays steady and slow decay. This result is consistent with the prediction that the system should flow to an insulating fixed point where no energy can be transmitted. The latter appears for and displays a current which increases with time. Again, this is consistent with the renormalization-group prediction which states that the perturbation at the center is irrelevant.

Let us finally consider the fact that the perturbation in Eq. (14) is the result of the abruptness of the junction and would not appear if and would vary smoothly in space (see Sec. III). To this aim, we numerically investigate a system where the anisotropy parameter is space dependent, , and linearly interpolates between the values and in a region of length located in the center of the system. is shown for several values of in Fig. 7 (top). As increases, and thus the junction is smoothened, the time dependence of the energy current transforms from an apparent decay behaviour to a more quasi-steady one. Intuitively, the smooth variation of in space implies a smooth variation of the local Fermi velocity , as shown in Fig. 7 (bottom). Since we have pointed out that matching the Fermi velocities is a sufficient condition to guarantee the irrelevance of the perturbation at the junction, this is consistent with the observation of a non-decaying current for large values of .

## V General scenario

We are now in the position to integrate the previous results into a general framework for the equilibration dynamics of inhomogeneous chains consisting of two different integrable models. We remark that our problem is a particular kind of local quantum quench, where the initial state differs from a steady state of the evolving Hamiltonian only at the junction. The local unitary dynamics induced by abides by a Lieb-Robinson bound, namely by a spreading of correlations which is constrained to be within a light cone with typical velocity , propagating from the center. Outside the light-cone, the state at time is indistinguishable from the initial one apart from exponentially small corrections. Generally speaking, the aim of this section is to describe the dynamics within the light-cone both for the time regimes accessible with our numerics and for the asymptotical stationary properties.

### v.1 Pre-equilibration behaviour

Let us first consider the time regime which is accessible with our numerics. The results in Figs. 5, 6 and 7 show that after a very short transient behaviour of order , the energy current is quasi-stationary, with relative variations which are negligible within the considered time window.

In Fig. 8 (top) we further elaborate on this point and show the current profile (where is the expectation value over the state at time ) for several values of . The propagation is asymmetric due to the different properties of and . Moreover, the current around the junction rapidly becomes homogeneous and quasi-stationary. This is explicitly verified in Fig. 8 (bottom), where the averaged current over sites approaches , while both become almost time-independent.

In between this central region and the unperturbed areas outside the light cone, two transient zones appear: the current here is not nor zero, as at the initial time. Given a time we select ; if we neglect these transient regions, it is possible to draw a connection between and defined in Eq. (7). Indeed, because of the homogeneity of the current profile discussed above, it is possible to make the following approximation: ; since is quasi-stationary the following holds:

 J(t)∼⟨∂t^ja,b⟩tvl+vr. (21)

Interestingly, it is possible to give a more explicit expression to using Eq. (8) and the Lieb-Robinson bound, which lead to the following formula:

 Extra close brace or missing open brace (22)

We remark that Eq. (22) is true up to exponentially small corrections.

In Fig. 10 we plot and : the comparison shows interesting analogies which seem to hold even though the transient regions, plotted in Fig. 8 (top), are significant. In general, Eq. (21) is well justified whenever the transient regions grow as with , introducing only subleading corrections (see for example the case discussed in Ref. (51), but note that for free fermions this is not verified (34)).

Let us compare these results with those presented in the literature for homogeneous integrable systems (15); (16). If and  (26), the operator and thus , with slope given by . This slope can be shown to give an exact lower bound to the stationary energy current at the junction  (27). This simple case corresponds to a perfectly transmitting junction, as expected because the system is homogeneous.

In general, the behaviour of allows us to discern among more complicated scenarios in which transmission might not be perfect. In particular, whenever is quasi-stationary in time and different from , it is intriguing to interpret it as an impurity determining a finite transmission because the matching between the modes of the halves is not perfect (31). This is for instance what is discussed in Sec. III, where the inhomogeneous low-energy theory of the system introduces a transmission coefficient and thus a non-optimal transmission of energy.

### v.2 Asymptotic long-time regime

Let us now move to the study of the system in the asymptotic long-time regime which is not accessible with our numerics. Although nobody has yet identified the minimal conditions for thermalisation, the expectation is that generic systems (e.g. exhibiting Wigner-Dyson level spacing statistics) should thermalize. We have argued that at long times the current at the junction decays to zero (see Sec. II). Whereas this steady-state property is a direct consequence of the assumption of thermalisation, it is possible to further elaborate on how the system approaches this situation. In particular, we will now elaborate on the existence of two domains of growing size which (i) are characterized by a non-zero energy current and (ii) propagate from the junction towards the two extremities of the chain (see the sketch in Fig. 9).

We begin by showing that must grow linearly in time for and . Starting from the definition, the time-derivative reads:

 Extra close brace or missing open brace (23)

In this equation, by choosing properly, can be taken large enough to ensure to have stationary and thermal behaviour. Assuming to be at times so long that thermalisation in the region around has occurred,

 Extra close brace or missing open brace (24)

We postpone the detailed analysis on the exact value of to Sec. V.3, here it suffices to assume the reasonable inequality . Recalling that is still outside the light cone and thus is a thermal expectation value with inverse temperature , the following holds:

 Extra close brace or missing open brace (25)

and the integrated current increases linearly in time as long as the wavefront has not reached . A simple analysis shows that this happens also on the right half of the system in a similarly defined - region.

How is it possible to make this result compatible with the fact that for long times ? The simplest scenario entails two domains which propagate in opposite directions, respectively in the right and in the left half. The size of each domain naturally grows because its two edges propagate due to different physical phenomena. The most external edge moves with the wavefront velocity at which energy is spreading in the chain. The most internal one is related to the spread of thermalisation originating from the junction. Thus, even if at short times the quenched dynamics resembles that of the homogeneous case, as discussed in Sec. V.1, the integrability-breaking induces thermalisation and accordingly a slow decay to zero of . This feature first appears at the junction and then affects all the other sites of chain.

A last comment is in order. Since , using the previous results we obtain that also grows linearly in time. Indeed, is a bounded operator that cannot compensate the linear growth in time of the other two terms. Since it has been shown that grows linearly in time also for the homogeneous chain, where , this indicates that there is no general connection between and in the asymptotic long-time limit. The discussion related to Eq. (21) is valid only for intermediate times where quasi-stationarity appears and is still different from zero.

### v.3 The problem of the effective temperature

We now consider the problem of identifying . In highly inhomogeneous scenarios, such as the one considered in this article, different regions of the system may experience thermalisation both (i) at different times and (ii) at different effective temperatures. Here, for example, the asymptotic extrema of the chain are always thermal at their initial temperatures; additionally, we expect thermalisation at to occur first at the junction and then propagate towards the edges. In similar situations with well-defined initial temperatures and unitary dynamics, the problem has been approached within the framework on Onsager kinetic equation (54). In order to link these methods to our microscopic model, it would be necessary to derive, at a coarse-grained level, a hydrodynamical equation governing the temperature dynamics in our system. Here we limit ourselves to some general considerations related to some fundamental constraints for a microscopic theory of thermalisation in a inhomogeneous system.

In typical quench problems, the conservation of energy (the time evolution is unitary) makes it natural to define through the following implicit equation:

 Tr[^Hle−βl^Hl]Tr[e−βl^Hl]+Tr[^Hre−βr^Hr]Tr% [e−βr^Hr]=Tr[^He−β∗^H]Tr[e−β∗^H]. (26)

However, if one considers a situation where the ratio between the lengths of the left and right halves is not , as assumed until now, but a different fraction, one realizes that the estimation of from Eq. (26) changes. This is absurd as the local dynamics does not involve the boundary and thus is not affected by this rescaling. Intuitively, this equation assumes that the stationary state inside the light cone is affected by the total amount of energy, and thus also by the fraction which lies outside it. However, this latter fraction could not play any role in the equilibration dynamics and cannot affect the value of the stationary temperature. This discussion highlights the fact that thermalisation is a local process which only in fine-tuned situations (e.g. completely homogeneous systems) can be addressed via global conserved quantities like Eq. (26).

From a practical point of view, one can envision to define through the long-time dynamics of a local observable sitting at the junction. If the stationary state is thermal-like, all local observables should agree with the thermal expectation value in the long-time regime.

## Vi Conclusions

In this work we studied energy transport between two (ideally) semi-infinite XXZ chains characterised by different coupling constants. The partitioning protocol considered here is identical to the one discussed in several papers recently appeared in the literature: the two half-chains are initially prepared in two separate thermal states and kept at different temperature. The chains are then connected, the system is let unitarily evolve and the energy current flowing through the system is analysed as a function of time.

Most of the paper is devoted to the case where each half-chain has uniform couplings; there is, however, a difference in their value between the left and right halves. Each half-chain is integrable, the whole system is not. This setup allowed us to study in details the scattering mechanisms regulating the energy flow in the system, as well as their relation to the integrability of the underlying model and to its conservation laws. In order to dwell more into the scattering of bosonic modes and quasi-particles at the junction, we also considered a situation in which the couplings are (adiabatically) varied through a finite region. The results that we presented are obtained both with a bosonisation approach and with numerical simulations based on time-dependent matrix product states.

The general scenario that emerges is quite rich including a pre-thermalisation regime and a final thermalisation. Since integrability is only locally broken in this model, a pre-equilibration behaviour may appear. In particular, when the sound velocities of the bosonic modes of the two halves match, the low-temperature energy current is almost stationary and is described by a formula with a non-universal prefactor interpreted as a transmission coefficient. Thermalisation, characterised by the absence of any energy flow, occurs only on longer time-scales which are not accessible with our numerics.

Quite interesting in this respect would be the investigation of the transition between the pre-equilibration and the asymptotic long-time regime: What is the typical time-scale of the thermalisation onset? How does it depend on the system parameters and on the initial temperatures? Another attractive direction would be the development of an analytical approach to this problem at high temperatures where the results obtained exploiting the low-energy field-theory are not reliable. It is not to be excluded that also the investigation of higher-order cumulants of the full counting statistics may provide valuable information (55).

###### Acknowledgements.
We acknowledge enlightening discussions with: J. Dubail, M. Haque, J. E. Moore, J-.M. Stephan and X. Zotos. AB, DR, LM, and RF acknowledge the Kavli Institute for Theoretical Physics, University of California, Santa Barbara (USA) for the hospitality and support during the completion of this work. This work was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915, by “Investissements d’Avenir” LabEX PALM (ANR-10-LABX-0039-PALM) (AD), by LabEX ENS-ICFP: ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL* (LM), by EU (IP-SIQS) and (STREP-QUIC) (RF), by the Italian MIUR (FIRB project RBFR12NLNA) (DR), and by Scuola Normale Superiore through the internal project “Non-equilibrium dynamics of one-dimensional quantum systems” (RF, DR, and LM).

## Appendix A Level spacing statistics

In this section we study how the level spacing statistics (LSS) of the XXZ chain is affected by the inhomogeneity discussed in the main text. We perform an exact diagonalization of systems made of spins in the zero magnetization sector ( states), and we analyzed the distribution of the normalized level spacings , where are the eigenvalues of (1) in ascending order (we discarded the lowest and the higher ), and is the average level spacing.

In the top panel of Fig. 11 we compare the results for the homogeneous case ( and ) with the inhomogeneous case with choosing and (which corresponds to the case). From our data we observe that the level spacing distribution moves from a Poissonian shape, , to a Wigner-Dyson distribution, , very quickly as soon as we switch on the anisotropy term . A Wigner-Dyson distribution of the LSS clearly signals the tendency toward level repulsion (52), a typical feature of non-integrable systems which seems to be present in this context at any value of the anisotropy strength .

In the lower panel of Fig. 11 we perform the same analysis for the same anisotropy parameters but tuning in order to get and we observe that the integrability is broken in the same way. In order to explore also the cases we set the anisotropy parameters as in Fig. 6. The results shown in Fig 12 indicate, also in these cases, a clear tendency toward a Wigner-Dyson distribution.

In conclusion, our LSS data indicate that in our system (1) the integrability is always broken in the inhomogeneous case , regardless of the value of and whether the renormalized Fermi velocities are matched or not.

## Appendix B Details on the ^Θ operator

In this Appendix we provide the explicit form of the operators and appearing in Eq. (8). Let us consider the integrated energy current operator as defined in Eq. (7). By computing the time derivative of such operator one gets Eq. (8), which is here reproduced for the sake of readability:

 i[^H,^ja,b]=^kl,a−^kr,b+^Θ. (27)

The operators and are linear combinations of two- and four-point operators with support on few sites around and respectively; is the sum of two-point operators localized around the junction.

Using the definition (5) of the local energy current and exploiting the canonical commutation relations (where is the Kronecker delta and is the Levi-Civita symbol), after some algebra we get:

 ^kl,a = Missing or unrecognized delimiter for \bigl (28) −^Sxa−2^Sya−1^Sya^Sxa+1−^Sya−2^Sxa−1^Sxa^Sya+1+^Sya−2^Sxa−1^Sya^Sxa+1+^Sxa−2^Sza−1^Sxa^Sza+1+^Sza−2^Sxa−1^Sza^Sxa+1+ +^Sya−2^Sza−1^Sya^Sza+1+^Sza−2^Sya−1^Sza^Sya+1)+^Sza−1^Sza)−4(^Sxa−2^Sza−1^Sza^Sxa+1+^Sya−2^Sza−1^Sza^Sya+1)]; ^kr,b = Missing or unrecognized delimiter for \bigl (29) −^Sxb−1^Syb^Syb+1^Sxb+2−^Syb−1^Sxb^Sxb+1^Syb+2+^Syb−1^Sxb^Syb+1^Sxb+2+^Sxb−1^Szb^Sxb+1^Szb+2+^Szb−1^Sxb^Szb+1^Sxb+2+ +^Syb−1^Szb^Syb+1^Szb+2+^Szb−1^Syb^Szb+1^Syb+2)+^Szb^Szb+1)−4(^Sxb−1^Szb^Szb+1^Sxb+2+^Syb−1^Szb^Szb+1^Syb+2)]; ^Θ = Missing or unrecognized delimiter for \bigr (30) Missing or unrecognized delimiter for \bigl +2JlJr^Sy0^Sy2(ΔlJl−ΔrJr)+2JlJr^Sz0^Sz2(ΔrJl−ΔlJr)+2ΔrJr(Jr−Jl)(Jl+Jr)^Sz1^Sz2].

The cumbersome expression for is significantly simplified when :

 ^Θ= −J3r4(Δl−Δr)[(Δl+Δr)(^Sx0^Sx1+^Sx1^Sx2+ +^Sy0^Sy1+^Sy1^Sy2)−2(^Sx0^Sx2+^Sy0^Sy2−^Sz0^Sz2)]. (31)

Additionally, when and (homogeneous case) we obtain .

## Appendix C Details on analytical calculations

Here we detail the computation leading to the result (16). The equation (19) is the action of the left-to-right transfer matrix on the fields and . In order to have the same normalization for the stress-energy tensor on the left and on the right one has to redefine (53) the bosonic fields as

 ϕ(x)={√Klφ(x)x<0√Krφ(x)x≥0 (32)

Further, introducing the notation , we can rewrite (19) as

 [∂xφl\par∂tφl]=[α1/200α−1/2][∂xφr\par∂tφr]. (33)

Notice that the equation above ensures the continuity of the momentum density through the junction , with . Total momentum is however not conserved by the dynamics.

To connect with the CFT formalism we must introduce the chiral waves and and recast (33) in the form a scattering matrix, see Fig. 13.

One finds

 [∂φr\par¯∂φl]=[cosγsinγ−sinγcosγ][∂φl\par¯∂φr]. (34)

with . The matrix in (34) can be interpreted as the representation of the action of the dual defect map on the incoming waves (31). To compute the energy current one needs to apply the composition of maps to the momentum operator and evaluate the result on the initial state , where the two free bosonic CFTs are in a thermal Gibbs ensemble at inverse temperature . We recall that the stress-energy tensor components are and and that the map acts as a pure reflection on the chiral incoming waves

 Ω0[∂φr\par¯∂φl]=[01−10][∂φl\par¯∂φr]. (35)

Since we assumed the asymptotic value of the current -independent (see however the discussion at the end of Sec. III) we can can take and obtain by simple manipulation and

 Ω−10∘Ω[¯Tl]=cos2γ\leavevmode\nobreak Tr+sin2γ\leavevmode\nobreak ¯Tl+2cosγsinγ\leavevmode\nobreak ∂ϕr¯∂ϕl. (36)

Then we can compute the energy current as (31)

 JLE=Tr[ρ0Ω−10∘Ω[Tl−