Ballistic transport and boundary resistances in inhomogeneous quantum spin chains

Transport phenomena are central to physics, and transport in the many-body and fully-quantum regime is attracting an increasing amount of attention. It has been recently revealed that some quantum spin chains support ballistic transport of excitations at all energies. However, when joining two semi-infinite ballistic parts, such as the XX and XXZ spin-1/2 models, our understanding suddenly becomes less established. Employing a matrix-product-state ansatz of the wavefunction, we study the relaxation dynamics in this latter case. Here we show that it takes place inside a light cone, within which two qualitatively different regions coexist: an inner one with a strong tendency towards thermalization, and an outer one supporting ballistic transport. We comment on the possibility that even at infinite time the system supports stationary currents and displays a non-zero Kapitza boundary resistance. Our study paves the way to the analysis of the interplay between transport, integrability, and local defects.


I. INTRODUCTION
In 1941, P. Kapitza reported on the first measurement of the temperature drop near the boundary between helium and a solid when heat flows across the boundary [1].The phenomenon, unrelated to that of conventional contact resistances, is ascribed to the mismatch between the energy carriers of the two materials, and appears even if the interface is perfect at the atomic scale.In the years 1951-1956, Kalatchnikov first, and Mazo and Onsager later, independently developed the so-called acoustic mismatch model, that gives a mathematical formulation to such intuition [2,3].The quantitative comparison between experimental data and theoretical predictions motivated extensive investigations even in other interfaces without helium, e.g., between two solids [2], where the phenomenon has been often recovered.
Owing to the difficulty of performing numerical simulations in order to benchmark the theory, simpler treatable models have been considered.In particular, focusing on junctions of classical one-dimensional (1D) harmonic chains enabled to quantitatively investigate the phenomenon of thermal boundary resistances [4][5][6].The effect has been widely reproduced, and comparisons with suitable adaptations of the theory have been presented.Yet, a numerical investigation of the phenomenon in the fully quantum case has not been performed so far [7].
As a paradigmatic example, we study a setup composed of two different semi-infinite spin-1/2 chains connected through a junction at x = 0 (see Fig. 1) [30].In order to highlight the resistive effects taking place at the junction, we consider two models for which ballistic transport of spin and energy has been unambiguously demonstrated.For x < 0, the spin chain can be mapped to non-interacting fermions: it supports ballistic transport, because the back-scattering phenomena due to interactions are impossible [33,34].For x > 0, we take a class of models that are solvable through Bethe ansatz (and thus integrable), where transport phenomena can be modeled by means of a generalized hydrodynamic theory (GHD) [35,36].The latter predicts ballistic transport, notwithstanding the presence of particle-particle interactions, due to elastic scattering without backscattering [37].While the dynamics of both halves can be described in terms of stable quasiparticles, the microscopic nature of such quasiparticles is different and therefore a nontrivial scattering dynamics is expected at the junction.
We present a detailed numerical study of the long-time behaviour of this setup, by considering the following protocol: (i) we prepare the system in a pure state where the two halves x ≶ 0 have an extensive imbalance of a global An inhomogeneous quantum spin-1/2 chain composed of an XX and of an XXZ half is initialized in a pure state with energy or magnetization imbalance and subsequently evolved in time.The relaxation dynamics inside the light cone (thick black lines, velocities vL,R) witnesses the coexistence of two (intermediate) ballistic regions and of two slow relaxation regions.The slow dynamics initially takes place inside an internal light cone with velocities v L,R , while, at long times, the extension of these regions remains constant in time.Bottom: Space-time profile of the energy current for ∆ = cos(π/4)see Eq. (1).At long enough times, a decay of the current at the junction appears as a consequence of a thermalization process.
conserved quantity; (ii) we unitarily evolve the state and monitor how energy and magnetization are redistributed in time.In particular, we focus on two very representative examples, where either the magnetization or the energy are macroscopically different within the two halves.Pure states are numerically more tractable than density matrices and allow us to explore longer time dynamics.
The emerging scenario is sketched in the top panel of Fig. 1 and summarized below.As already observed in other setups, two wavefronts propagating at different velocities v L and v R emerge from the junction [30,36,38,39]; v L and v R correspond to the velocities of the fastest quasiparticle excitations of the left and right halves, respectively.Outside this light-cone, in agreement with the Lieb-Robinson bound [40], regions 1L and 1R display the initial equilibrium behaviour and are unaffected by the dynamics: these are the original equilibrium regions.Conversely, within the causal region, the system properties are nontrivially affected by the unitary dynamics.
As a first key result we found that, when two different ballistic models are joined together, the relaxation dynamics inside the light-cone exhibits qualitatively different behaviours, depending on whether we are observing the system near or far from the junction: close to the edges of the light-cone (regions 2L and 2R in Fig. 1), a stationary state supporting a stable current flow is rapidly approached; these regions are accordingly dubbed intermediate ballistic regions.Otherwise, around the junction, the current intensity keeps decreasing without reaching any stationary value even at the largest accessible time in our numerics.As a matter of fact, regions 3L and 3R are characterized by a slow relaxation dynamics and are named slow relaxation regions.Let us stress that in the uniform-Hamiltonian setups the relaxation usually occurs first around the junction [36,38,39,41].
Here we observe the opposite behaviour: a slowly decreasing current-flow which is a direct consequence of a sort of thermalization process that originates at the junction.As an example of data supporting our conclusions, in the bottom panel of Fig. 1 we provide the space-time profile of the energy current for a paradigmatic case (see Results for further details).
We can interpret our results using a semiclassical picture of quasiparticles that travel freely inside each half, but undergo scattering events with frequency 1/τ L,R in a region [− L , R ] around the junction.Fast quasiparticles spend less time around the junction and are therefore less affected by these processes; this creates two qualitative space-time regions, separated by interfaces at The regions 2L and 2R are forming because of the fast quasiparticles which pass almost unchanged through the junction: this justifies a quick relaxation to a stationary state on each ray at constant x/t, similarly to the homogeneous case.On the other hand, regions 3L and 3R are characterized by several scattering events, and display a strong tendency towards equilibration.
As a second crucial result, we are able to extrapolate data at long time.We argue that the standard scenario with thermalization in the whole system, leading eventually to vanishing energy and spin currents, is incompatible with the ballistic transport in the bulk of each half.At long times, the system exhibits a local quasistationary state (LQSS) on each ray x/t, as it happens for homogeneous integrable ballistic systems, but with a non-trivial behaviour in a region around the junction, whose width does not scale with time.As a result, the two emerging velocities v L,R tend to zero at long times (i.e. for t | L,R /v L,R |), as depicted in the top panel of Fig. 1.
The emerging asymptotic scenario is corroborated by considering an Ohmic model where the two semi-infinite spin chains are coupled via a diffusive junction of fixed length .The long-time dynamics of such a model can be solved within the hydrodynamic framework that describes ballistic transport in each half, and the resulting stationary profiles are in qualitative agreement with the exact numerical simulations.We conclude by stressing that a sharp discontinuity must appear in several local quantities at x/t → 0 ± .This phenomenology is analogous to that of the thermal Kapitza boundary resistances, which is here witnessed in a fully-quantum scenario without important approximations.

II. THE MODEL
We consider a quantum spin-1/2 chain described by the Hamiltonian Ĥ = n ĥn : where Ŝα n is the operator associated to the α component of the n-th spin (hereafter we adopt units of = k B = 1).The anisotropy parameter is space dependent: for n ≤ 0, ∆ n = 0 whereas for n > 0, ∆ n = ∆.This model describes a perfect junction between a XX model and a XXZ model; a Jordan-Wigner transformation maps the former into non-interacting spinless fermions and the latter into interacting ones.We fix ∆ = cos(γ) ∈ [0, 1) so as to have a gapless model whose transport properties are well understood in the uniform regime.
At t = 0, the system is initialized in a pure state.We consider two kinds of initial states.The first one is the ground state |Ψ 1 of the Hamiltonian Ĥ = n λ n ĥn , where λ n = −1 for n ≤ 0 and λ n = 1 for n > 0. In this way, for n ≤ 0 the state locally approximates the maximally excited state of the XX chain, whereas for n > 0 it approximates the ground state of the XXZ chain.The state |Ψ 1 is thus the starting point of a partitioning protocol with two different inverse temperatures, β L and β R , in the limit β L,R → ∓∞.In this way, the energy imbalance is maximal and the state is pure.The second one is the domain-wall state [39,42,43] |Ψ 2 = ⊗ n |ψ n , with |ψ n = |↑ for n ≤ 0 and |ψ n = |↓ for n > 0, which corresponds to maximal magnetization imbalance.Thereafter, the system evolves unitarily as |Ψ 1,2 (t) = e −i Ĥt |Ψ 1,2 , and allows to use matrix-product-statesbased algorithms (see App. A) so as to speed up the computation and explore large enough times to have access to novel equilibration regimes [44].We focus on two different sets of observables: in the first case, we look at the energy transport, through the local energy density, E n = ĥn , and energy current, Ĵ(E)  [ ĥn , ĥn−1 ] .In the second one, we focus on the spin transport and measure the local magnetization, S z n = Ŝz n , and the associated current Ĵ(M) Ŝy n .Throughout this article, time is measured in units of J −1 , energy in units of J, energy current in units of J 2 and magnetization current in units of J.
Since the dynamics takes place inside a light cone (see Fig. 1), we will mostly present our numerical data using rescaled space-time units x/t.As an aside, this also permits to translate the problem in the standard framework of transport between two reservoirs held at fixed distance [10,45].

III. RESULTS
We begin with the study of energy transport in the first situation, where the system is initially prepared in |Ψ 1 .In Fig. 2a, we plot the local energy density profile E n (t) for γ = π/4.For all the displayed times, the energy profile is significantly different from that of the homogeneous case (γ = π/2).This is particularly evi-dent when analyzing the profile for rays |n|/t 0.5 close to the junction: energy is redistributed in the form of two approximately flat plateaus, which develop at different values, depending on ∆, for n/t ≶ 0. Such a result is systematically present as soon as γ = π/2 (∆ = 0).This is shown in Fig. 2b, where we plot the energy profile at the largest accessible time, t = 40, for several values of γ: the plateaus appear whenever ∆ = 0 and their width increases monotonically with ∆ > 0.
A similar study can be performed for the domain-wall initial state |Ψ 2 , when focusing on the spin transport; as a matter of fact, the results are analogous to those obtained above for the energy transport.The magnetization profile is displayed in Fig. 2c for γ = π/4 and several values of time, and in Fig. 2d for several values of γ and t = 40.Fig. 2d undoubtedly shows that, as soon as γ = π/2, and thus the right half is described by a Hamiltonian which is different from that of the left half, the magnetization of the left chain deviates from the behaviour of the homogeneous system (γ = π/2) and exhibits a slow relaxation in the region −0.5 n/t < 0.
Let us start our analysis with a discussion of the re-gions 0.5 |n|/t < 1, that are identified with the intermediate ballistic regions 2L and 2R.In Figs.2a and 2c it is shown that here the energy and magnetization profiles reach a stationary value.In Fig. 2b and 2d, we compare the long-time behaviour for different values of γ in the right half.We observe that in region 2L the energy and magnetization densities are unchanged and coincide with that of the free homogeneous case.This fact is rather surprising, because the energy flow in this region originated by particles emitted from the state in the right half, whose nature does depend on γ.
We can address this phenomenon by means of a semiclassical and intuitive picture: since this γ-independent region emerges for sufficiently large rays, it has to be ascribed to the properties of the fastest quasiparticles, that because of their velocity, are essentially transmitted by the junction.However, the fact that the energy profile does not depend on γ, means that the junction turns them into the fastest carriers of the left half with the same speed.If this interpretation is correct, all the properties of the region n/t −0.5 (which corresponds to the regions 2L and 1L) should be well-described by the dynamics of the homogeneous case ∆ = 0, for which the system is free and exactly solvable.
In order to verify the latter statement, in Fig. 3a we present our results for the energy-current profile at some specific values of time, and γ = π/4.A more complete color-plot of the full space-time of the energy current is provided in Fig. 1, bottom panel.We also superimpose the energy-current profile of the free homogeneous system with the same initial temperatures β L → −∞ and β R → ∞ (dashed lines).The predictive power of this simple interpretation of the region 2L is remarkable.Moreover, it is not restricted to the left noninteracting half: considering the homogeneous situation with ∆ = cos(π/4) on the whole chain gives a good description of the current dynamics for n/t 0.5 as well (also called regions 2R and 1R).The same statement is true for the case of spin transport, as well.In Fig. 3b, the two 2L and 2R regions are once again well described by the corresponding homogeneous problems initialized in the domain-wall state |Ψ 2 (dashed lines).
We now move to the central regions −0.5 n/t < 0 and 0 < n/t 0.5, that are identified as the slow relaxation regions 3L and 3R, respectively.In Figs.3a and 3b, we observe that around the junction the relaxation dynamics is slow: the current intensity maintains a monotonically decreasing trend even at the longest accessible times of our numerics.We stress that the rate of such decay is much slower than the microscopic energy scale J in Eq. (1).
We thus analyze how the slow-relaxation regions expand in time.We look at the energy transport setup, but our findings remain valid when looking at the magnetization imbalance as well.As a quantifier of its extension, we track the time-evolution of the maximum of the energy current profile in Fig. 3a, which is plotted in Fig. 3c.The emergence of an internal light cone after a transient time is evident.In Fig. 3d, we show the dependence on ∆ of the typical velocity v R in the interacting right half; the data are compatible with a linear behaviour passing through the origin.In other words, when ∆ → 0, there is no region of slow dynamics, as expected for the uniform setup.This constitutes an important consistency check of our analysis.We remark that v L,R can be defined only at intermediate time scales while v L,R → 0 at asymptotically long times.

A. Absence of complete thermalization at the junction
In order to explain the slow dynamics around the junction, a natural possibility is that the system is undergoing thermalization; such a scenario would be motivated by the fact that integrability is always broken for γ = π/2 [30].However, it is undeniable that, at the time scales accessible with our numerics, relaxation did not yet fully take place.We now discuss how the numerically observed scenario may evolve at later times and whether it could be consistent with thermalization.In the following discussion, the thermodynamic limit is assumed, or anal- , so that the boundary plays no role.
The relaxation dynamics of closed many-body quantum systems has been thoroughly explored in the last decade [8][9][10][11][12].Several studies in homogeneous settings have shown that equilibration generically occurs, and that the stationary state is a generalized Gibbs ensemble (GGE).When the model is integrable, it must take into account the full set of local conserved quantities of the model.Since our model is not integrable, the corresponding ensemble is solely characterized by temperature and magnetic field.However, deep in the bulk of each half and far from the junction, the former approach could be applied also to our model.We now show that this possibility forbids a complete thermalization of the system at the junction and spreading from the junction.In particular, we demonstrate that it is not possible to fully describe the junction with a single well-defined temperature, and thus to assume that for instance no current flows there.
We introduce the concept of rapidity ζ x/t and consider the limit of infinite times at fixed rapidity ζ.We study the expectation value of a generic local observable Ôn : Here ρ ζ characterizes the LQSS, the density matrix on each ray of the space-time with rapidity ζ.Since the state is integrable within each half, we assume that for any ζ < 0 (ζ > 0), ρ ζ takes the form of a GGE expressed in terms of the left (right) conserved quantities.This observation is not sufficient to identify the specific form of GGE; however, using the integrability of each half, we obtain that the expectation value in Eq. ( 2) should be a smooth function of ζ.This allows us to uniquely identify ρ ζ once the appropriate boundary conditions for ζ → 0 ± and for ζ → ±∞ are found (see App. C for more details).
Clearly, if ζ < v L or ζ > v R , we are inspecting the state outside the causal region; since here no dynamics can occur, ρ ζ must coincide with the initial state.If we consider ζ = 0, namely at any fixed but arbitrary distance n from the junction, and assume that the junction lets the whole system thermalize, a thermal state must emerge at large times: Here, Z is the partition function while β and h are the uniquely defined effective inverse temperature and magnetic field, respectively.In practice, the density matrix ρ ζ must interpolate smoothly between the initial states outside the light cone and the thermal state at the junction.In particular: We stress that the two limits in Eq. ( 4) are different because the local Hamiltonian ĥn defined in Eq. ( 1) depends on the sign of n.
Yet, one such solution is incompatible with the free propagation of quasiparticles within each half, as it can be seen from the following argument.Let us consider for simplicity the left half and take a rapidity ζ which is negative and small.At large times on this ray, the left movers will be coming from the junction described by ρ TH ; instead, the right movers will have propagated almost freely from the initial state on the left.One such LQSS is not thermal because in the thermal case left and right movers should be characterized by the same temperatures.This remains true taking the limit ζ → 0 − .Thus, he corresponding LQSS at ζ → 0 − will not be a thermal state and Eq.(4a) is clearly violated (see App. C for more details).
We conclude that Eqs.(3) cannot hold in this setting.As a consequence, we expect that at larger times, even though integrability is broken, the system will equilibrate to a non-trivial steady-state.Our results point at the existence of a finite region around the junction that does not spread in time and that connects the ζ = 0 − and ζ = 0 + GGEs without widening in time.This scenario is consistent with those of Refs.[32,46] in a related setting.

B. A solvable model for the junction
In the previous section, we argued that the microscopic details of the junction determine how the GGE on the left is connected to the one on the right, and that this connection takes place in a finite region that does not scale with time.Before discussing whether such picture is consistent with the numerical simulations presented so far, in this section we consider a simpler model where a quantitative inspection of the behaviour of the junction at long times is possible.
Instead of connecting the left and right Hamiltonians with a two-site junction, we insert a finite region of length centered around n = 0, whose Hamiltonian Ĥchaos is fully chaotic (see Fig. 4).Such Hamiltonian is not integrable, is disordered and features diffusive transport, as in Ref. [47].The idea is that locally the system relaxes in finite (and not extensive) time scales τ rel , and that τ rel 2 /D, with D the diffusion constant.We still assume that the total magnetization is conserved, i.e. [ Ĥchaos , Ŝz ] = 0, though for simplicity we will only focus on energy transport.
The argument of the previous section still applies, so that no overall thermalization shall be expected.Nevertheless, if the value of is large, we can once again employ a hydrodynamic description.Taking a coarse-grained coordinate x ∈ [− /2, /2], we assume local equilibration to a thermal ensemble characterized by a space-time dependent temperature T (x, t).Then, the two fundamental equations describing temperature and energy-current dynamics in this region follow from the heat equation and Fourier's law, and read: with α(T ) the thermal diffusivity of the junction and κ(T ) the thermal conductivity at temperature T .In this hydrodynamic picture, the microscopic details of Ĥchaos are encoded in α(T ) and κ(T ); for weak temperature variations, we can assume that they are constant, namely α(T ) = α 0 and κ(T ) = κ 0 .In the stationary limit, the temperature profile depends linearly on x and the energy current is uniform.The stationary value of J E depends on the temperatures at the edge of the junction T ± : This last equation provides the definition of a thermal Kapitza boundary resistance: We stress that, although we use the notion of temperature to describe the region around the junction, this is not invalidating the results of the previous section.Here, a current is flowing and thus thermalization did not completely take place; in this context, the notion of a temperature is thus an approximation which becomes more and more accurate as → ∞.We use GHD to describe energy transport in the two halves for any non-vanishing rapidity ζ = x/t and bridge the original properties outside the light cone with the obtained solution around the junction.In practice, one looks for a solution in the left and right half that (i) has the proper behaviour outside the light-cone, and (ii) has the same energy current [see Eq. (6b)] in the limit ζ → 0 ± .The problem still remains unconstrained: R is a free parameter of the junction model; once it is fixed, the temperatures T ± are determined.
Although this model based on a chaotic junction has not been introduced in order to explain the numerical data in Figs. 2 and 3, we investigate whether it can provide an effective description of such data.In Fig. 4b, we plot the energy and current profile for different values of the resistance R for the initial case |Ψ 1 obtained with the chaotic junction model.The energy current is continuous by construction, but the energy profile has a discontinuity.This discontinuity is reminiscent of the behaviour observed in Fig. 2; for a better readability, the numerical data for γ = π/4 at the longest accessible time are superimposed.A qualitative agreement for the energy profile is observed.We observe that the analytical solution washes out completely the distinction between regions 2 and 3, merging them in a unique region inside the light-cone similar to regions 2.

C. Thermalization dynamics and stationary transport
In the spirit of the previous section, we investigate whether the state around the junction can be successfully described as a thermal state, ρ ∝ e −β( Ĥ−h Ŝz ) .We focus for simplicity only on the energy-transport scenario, leaving for brevity the data for spin transport to the App.E. In this case, thanks to the spin-flip invariance S n z → −S n z of the initial state, the Gibbs state is characterized by a zero magnetic field, and only a temperature has to be identified.We associate an effective inverse temperature β eff n (t) to each set of three neighbouring sites by considering the local energy density ε n for the sites n, n+1, n+2, and inverting the thermodynamic relation β(ε) that associates a temperature to the corresponding energy density ε (see App. B).We stress that by doing so, the inverse temperature β eff n (t) is both space-and time-dependent, as in the hydrodynamic picture presented above.The result is plotted in Fig. 5, panels (a) and (b); we observe that the inverse-temperature profiles recall in their qualitative features those of the energy profiles plotted in Fig. 2.
To quantify how well the temperature β eff n (t) captures the whole set of local properties of the state, which is what happens when thermalization occurs, we compute the operator distance between the three-site reduced density matrix of |Ψ 1 (t) , dubbed ρ n,n+2 (t), and the three-site reduced density matrix of the thermal state e −β eff n (t) Ĥ /Z, dubbed ρ n,n+2 [β eff n (t)] [48].As an indicator for the distance, we use the trace norm (or Schatten norm for p = 1): d(ρ 1 , ρ 2 ) = ρ 1 − ρ 2 1 , that coincides with the sum of the singular values of the difference ρ 1 − ρ 2 .
In Fig. 5c, we analyze the distance for several values of time and γ = π/4.The plot displays the same qualitative features of the energy profiles in Fig. 2a.Although a stationary state has not been reached, the region around the junction displays a clear tendency towards equilibration, that is signaled by a pronounced drop of the value of the distance in time.This is more clearly observed in Figs.6a and 6b, where we plot the time-dependence of Ĵ(E) n=0 (t) and Ĵ(M) n=0 (t) in the two situations; different values of ∆ are considered.In both cases, the current reaches a steady value for ∆ = 0 (γ = π/2) but, as soon as ∆ = 0, it shows a slow decreasing trend; the effect becomes more pronounced with increasing ∆.On the other hand, the regions 0.5 |n|/t 1 do not display any tendency towards equilibrium at all.The original equilibrium regions 1L and 1R have the minimal value of such distance, and as expected it does not depend on time (not shown).In Fig. 5d, we analyze the time dependence of such distance for four fixed values of n, which are taken close to the junction.The plot shows that the relaxation behaviour depends on which side of the junction is considered.
Inspecting the local effective temperatures in Fig. 5, panels (a) and (b), one sees that plateau values are reached around the junction which are different in the left and right parts.The interpolation between the two values takes place on a length scale of few sites and does not depend on time.This agrees with the simple picture of a boundary Kapitza resistance introduced in Eqs. ( 6) and (7).As a consequence, a stable discontinuity builds up, thus guaranteeing the energy current at the junction not to vanish.Let us stress that a similar situation takes place in the case of magnetization transport, and in that case we can speak of a magnetic Kapitza boundary resistance.This phenomenology has been widely observed in classical 1D systems and several mesoscopic transport experiments [5,6], while here it is first presented in a fully-quantum description.

D. Entanglement production
The equilibrating phenomenology that emerges from the junction is reflected by other measurable properties of the system, such as the bipartite entanglement [49,50].We take advantage of the fact that our system is pure, so that the von Neumann entropy of the reduced density matrix of the left part In Fig. 6c and 6d, we plot S(t) for several values of γ.Both in the case of energy and spin transport, the growth of S(t) is logarithmic for the homogeneous case γ = π/2: this is a consequence of integrability which allows evolving the initial state with a weak generation of entanglement on each ray [51][52][53][54].Indeed, as soon as the system becomes inhomogeneous, the behavior of the entropy qualitatively changes and grows superlogarithmically, testifying the process that is taking place at the junction.The data are particularly clear for the case of spin transport, where the initial state is an uncorrelated product state.

IV. CONCLUSIONS
Ballistic transport of particles and energy in 1D quantum systems is a fundamental concept in non-equilibrium physics.This phenomenology has been observed experimentally in a wide spectrum of physical platforms ranging from nanowires and carbon nanotubes to the edge states of Hall bars and topological insulators.In this work we addressed the problem of studying the emergence of Kapitza boundary resistances when two 1D quantum ballistic conductors are joined together.
Recent works have pointed out that transport phenomena are tightly bound to a number of fundamental issues in statistical mechanics.The study of equilibration in 1D closed quantum systems has recently assessed that homogeneous integrable models have a peculiar long-time dynamics because in these systems quasi-particles interact without suffering from backscattering.This allows for the presence of stationary ballistic currents and is accounted for by the GHD description of ballistic transport.
The situation becomes less trivial when the quasiparticles of the two conductors have different nature [30,55].From a theoretical viewpoint the system becomes nonintegrable and the known framework does not apply any longer.Yet, ballistic transport (a consequence of integrability) is supported in each separate half of the chain.In general, the community is currently considering the presence of forms of local integrability breaking terms, such as localized defects [30,32,46,[56][57][58]; our work fits within this research effort.
Although our system is not integrable, our conclusions are in sharp contrast with what one would expect for generic chaotic and diffusive systems, even when an inhomogeneous background is considered [59].For instance, the light cone characterized by velocities v L,R is due to the integrability of each half and is not expected to appear in fully non-integrable models, where rather a diffusive scenario is the anticipated situation.Moreover, our numerics witnesses in real time the appearance of a slow-relaxing region that expands from the center.This region is the direct consequence of the integrabilitybreaking properties of the junction.At asymptotically long times, we argue that the system eventually reaches a state which supports stationary currents even though integrability is broken.The two halves equilibrate to distinct GGEs, which are connected in a non-trivial way by the junction.
We conclude by stating that systems where integrability is broken only in a finite region of space are therefore special, even though they look non-integrable under conventional indicators such as the level spacing statistics [30,32].This demands a novel characterization of integrability (or the lack thereof) appropriate for inhomogeneous settings.
This work may be the object of experimental investigations in several platforms where coherent Hamiltonian dynamics can be easily realized, and the partitioning protocol employed in this work might reveal crucial for studying transport phenomena in systems where leads are not easily implemented.A natural possibility are ultra-cold atoms, where the realization of 1D lattice systems is possible and recent in-situ microscopy techniques allow for a site-resolved study of coherent time evolution [60].The recent developments of assembled quantum simulators motivate similar investigations in arrays of Rydberg atoms [61].Superconducting circuits have recently attained significant coherence times, and can realize one-dimensional systems with tunable (and in particular inhomogeneous) non-linearities [62].Here too our findings could be investigated.
Summation over repeated α indexes is assumed.The matrices A [σj ] have two indexes which label the bond link and can take at most χ values (the first and the last one are exceptional and have only one index, in order to enforce open boundary conditions).This latter parameter sets the precision of the ansatz and of the associated numerical simulation.Roughly speaking, the larger the bond link, the more correlations within the state can be faithfully described; a small bond link, instead, describes an almost uncorrelated state (χ = 1 is a product state).
The first step requires the preparation of one of the initial states |Ψ 1 or |Ψ 2 , detailed in the text.The former is prepared by a standard iterative ground-state search [16] of Hamiltonian (1) with λ n = −1 for n ≤ 0 and λ n = 1 for n > 0. The latter is a MPS with bond-link dimension χ = 1 which is constructed explicitly.
The system is subsequently evolved in time with the Hamiltonian Ĥ = n ĥn , where ĥn is given in Eq. (1).We employed the time-evolving blockdecimation (TEBD) algorithm with a fourth-order Trotter decomposition and time step dt = 0.02J −1 .We have checked that numerical inaccuracies due to this choice are negligible on the scales of all the figures shown here.
The bond-link dimension χ of the MPS representation increases exponentially in time because of the spreading of correlations.Furthermore, the fact that the system is inhomogeneous increases the complexity of the state, as already observed by a direct comparison of the results published in Refs.[30,63].Since a higher bond-link dimension implies a larger computational cost to evolve the MPS, an upper bound χ max of the representation needs to be set.Remarkably, in all the simulations performed here, convergence for observables under consideration for t ≤ 40 has been reached for χ max ≤ 500.In order to faithfully explore the large times considered in the paper we simulated a chain made of L = 140 sites with open boundary conditions.This guarantees that, for the parameters considered in this work, the dynamics is not affected by the boundaries.

Appendix B: Energy-temperature and magnetization-chemical potential relations
In order to ascribe an effective local temperature β eff n to the time-evolved state, one needs to compare the local energy of the system at time t with that of a system at thermal equilibrium.This quantity is both space (since our setup is not homogeneous) and time dependent.Let us stress that the state |Ψ 1 does not involve any spin transport since Ψ 1 | Ŝz n |Ψ 1 = 0, ∀n and the Hamiltonian specified by Eq. ( 1) conserves the total magnetization.As a consequence, we can restrict our analysis to states with zero effective local chemical potential.
At each time t, we may define the local energy at the site n by considering a subsystem of three contiguous sites {n, n + 1, n + 2}.The average energy is thus given by Ēn (t) = ĥn t + ĥn+1 t /2, where • t denotes the expectation value over the time-dependent state.Next, we consider a homogeneous system with the parameters of the left (if n ≤ 0) or right half (if n > 0) at thermal equilibrium, and compute the average energy density in the bulk E(β); note that β can be both positive and negative.Inverting this relation, we can associate to the system a local effective inverse temperature β eff n , for each time of the time evolution.
The reconstruction of the three-site reduced density matrix ρ n,n+2 (t) requires the calculation of 64 observables and has been performed only for γ = π/4 and γ = π/2.The latter is the homogeneous and noninteracting case and is used as a benchmark to highlight the peculiarities of our inhomogeneous setups (see App. F for details).For all other values of γ plotted in Fig. 5b, the average energy is simply computed as E n (t) = ĥn t .
The non-equilibrium dynamics generated from the domain-wall state |Ψ 2 is characterized by a non vanishing local energy E n (t) and magnetization S z n (t) = S z n t .Therefore, in addition to the local inverse temperature β eff n (t), we may identify an effective local chemical potential µ eff n (t) as well.For the gapless XXZ chain at finite temperature and chemical potential, these quantities can easily worked out by standard thermodynamic Bethe ansatz (TBA) [36,64].
We thus proceeded as follow: (i) for fixed β and µ, we associated to each point in space a TBA description in terms of root densities {ϑ n (λ)} depending on the local value of the anisotropy ∆. (ii) we computed, in that local TBA, both the energy density e n (β, µ) and the magnetization density s n (β, µ) for β ∈ [−5, 5], µ ∈ [−5, 5] and discretization step dβ = dµ = 0.1; (iii) after interpolating those functions, we were able to extract the time-dependent local effective temperature and chemical potential by numerically solving the two equations e n (β, µ) = E n (t) and s n (β, µ) = S z n (t).
Appendix C: Details about the impossibility of a complete thermalization We now provide some technical details on the proof of the absence of complete thermalization.We first focus on the left half of the system.As the model is integrable in the bulk, in the thermodynamic limit it admits an infinite and complete set of conserved densities {q j µ } ∞ µ=1 [65], whose support is localized around the site j (possibly with exponential tails [66]).Each density satisfies the operatorial continuity equation provided that j ≤ j * , where j * negative and sufficiently distant from the junction.Now consider the total amount of a given conserved quantity in the region of space Using (C1), the time derivative ∂ t Q µ (t; ζ) involves the expectation value of the current Tr[ Ĵµ j ρ 0 (t)] at the endpoints j = j * and j = −ζt, which can be computed using ( 2) and (3) For the XX model, this is enough to deduce Eq. (4a), because completeness of densities {q j µ } µ implies the one of currents [67].On the contrary, in the presence of interactions, such an inference is subtler.For integrable models supporting relativistic invariance, densities and currents are related by Lorentz transformation, Eq. ( C3) is suffcient to deduce Eq. ( 4).Here, following the arguments which led to GHD in [35], we assume that Eq. ( 4) must hold whenever Eq. (C3) holds.
However, the ballistic propagation of quasiparticles within each half is generically incompatible with Eq. ( 4).For the sake of simplicity we analyse the XX case.After a Jordan-Wigner transformation, the model is diagonalized by fermionic operators in Fourier space ĉk , ĉ † k , satisfying {ĉ † k , ĉk } = δ k,k .Then, a GGE state is completely characterized by the occupation number n(k) = Tr[ĉ † k c k ρ 0 (t)] at momentum k and we indicate with n (L) (k) the state corresponding to ρ L .Then, according to GHD, the occupation number is promoted to a space-time dependent function satisfying [31,68,69] Here, the velocity takes the simple form v(k) = d (k)/dk with the dispersion relation (k) = −J cos(k).Eq. (C4) needs to be solved in the domain x ∈ (−∞, 0), with the boundary conditions: n(k; x < 0, t = 0) = n (L) (k) and n(k; x = 0, t) = (1 + e β( (k)− h) ) −1 ≡ n TH (k), the equilibrium Fermi-Dirac distribution at temperature β and chemical potential h.The solution in the limit x/t = ζ → 0 − is easily found to be: Eq. ( 4) is clearly not satisfied, except for the trivial case n (L) (k) = n TH (k) where no dynamics occurs.In a similar way, one can extend the same argument to the interacting integrable case for x > 0 [70,71].This shows that Eq. ( 3) is inconsistent with the ballistic propagation of quasiparticles.
Appendix D: The chaotic model of the junction In this section, we provide additional details about the model of diffusive junction and how the data in Fig. 4 were obtained.As explained in the previous section, integrability of the two halves implies that they admit a description in terms of GHD.For the left half, Eq. (C4) applies.Instead, for the right half, the XXZ model with ∆ = 0 requires Eq. (C4) to be generalized as [36] ∂ Comparing with (C4), we changed the notation for the occupation number n → ϑ α and we employed a parametrization in terms of the rapidity λ instead of the momentum k in agreement with standard literature.For ∆ = cos(π/P ), the index α = 1, . . ., P labels the particle type and is due to the presence of stable boundstates (see [64] for generic values of ∆).A fundamental difference is that the velocity is dressed by the interactions and becomes a state-dependent quantity v α (λ; ϑ) = [ ] dr α (λ)/[k ] dr α (λ) [37], where α (λ) and k α (λ) are respectively the bare single-particle energy and momentum.The dressing q(λ) → [q] dr (λ) is a linearoperation defined for any set of functions q α (λ) of rapidities via D2) The kernel T α,β (λ, µ) depends on the interaction and we refer the reader to [36,38] for the explicit formulas as a function of ∆.Every GGE density matrix ρ GGE is characterized by the set of functions [ϑ α (λ)] P α=1 [38,[64][65][66].For instance, the expectation value of a conserved density qj and current Ĵj are written as [k ] dr α (λ)ϑ α (λ)q α (λ), (D3a) where q α (λ) are the single-particle eigenvalues corresponding to the conserved density qj .
In the study of thermal transport, Eqs.(C4, D1) in the region |x| > /2 can be used in combination with Eq. ( 5) in order to have a full hydrodynamic description.The initial state |Ψ 1 corresponds to the left/right GGE states [64] n (L) (k) = θ(|k| − π/2) , ϑ (R) α (λ) = δ α,1 . (D4) where θ(x) is the Heaviside step function.In practice, Eqs.(C4, D1) can be solved in the domain |x| > /2, while the junction acts as a source of thermal particles at x = ± /2 respectively at the instantaneous temperatures T (± /2, t).In the simplified case of constant diffusivity α 0 and conductivity κ 0 , we can directly look for a stationary solution at t → ∞.In this case, the stationary temperatures T ± in Eq. (6a) at the edges of the junction are the only unknowns.One can then solve Eq.Then, one can determine the values of T ± by imposing that the left/right energy currents computed via (D3b) at ζ = 0 ± match with Eq. (6b) and the only phenomenological parameter is the Kapitza resistance in Eq. ( 7).
The spin trasport from the state |Ψ 2 can be modeled in a similar way.However, in this case the initial magnetization imbalance induces also a non-vanishing energy current.A model of a diffusive junction in this case would require a 2 × 2 Onsager matrx of phenomenological coefficients.For simplicity, we do not consider this case here.

Appendix E: Thermalization dynamics for spin transport
In this section we show the thermalization dynamics for the spin transport scenario.In analogy to the analysis done in the main text for the energy transport scenario, we study the dynamics of the local effective chemical potential µ eff n (t) (see App. B for details).In the panel (a) of Fig. 7 we show the effective chemical potential profiles at different times for a fixed anisotropy.The data show the progressive emergence of a discontinuity at the junction.Again, as in the case of energy transport, this behavior can be understood in terms of a magnetic Kapitza boundary resistence.In panel (b) of Fig. 7 we show the profile of µ eff n for the longest time accessible with our numerics (i.e.t = 40) as a function of the anisotropy parameter γ.While in the homogenous case γ = π/2 the profile is smooth, for γ = π/2 the discontinuity increases as the anisotropy is increased.
Appendix F: Absence of thermalizing tendency in a homogeneous integrable model In this section we show the the absence of thermalizing tendency when γ = π/2 (homogeneous case).In analogy with the analysis performed in the main text for inhomogeneous systems, we study the distance between the local reduced three-site density matrix ρ n,n+2 (t) and the thermal density matrix ρ n,n+2  systems do not display any tendency towards thermalization at the junction (n/t = 0).The rescaled profiles of d n (t) at different times collapse for t 10.This also shows the absence of slow relaxation regions in homogenous systems which relax to a non-equilibrium stationary state on a faster timescale.

FIG. 1 .
FIG.1.Transport in inhomogeneous spin chains.Top: An inhomogeneous quantum spin-1/2 chain composed of an XX and of an XXZ half is initialized in a pure state with energy or magnetization imbalance and subsequently evolved in time.The relaxation dynamics inside the light cone (thick black lines, velocities vL,R) witnesses the coexistence of two (intermediate) ballistic regions and of two slow relaxation regions.The slow dynamics initially takes place inside an internal light cone with velocities v L,R , while, at long times, the extension of these regions remains constant in time.Bottom: Space-time profile of the energy current for ∆ = cos(π/4)see Eq. (1).At long enough times, a decay of the current at the junction appears as a consequence of a thermalization process.

FIG. 2 .
FIG. 2. Energy and magnetization profiles during the relaxation process.Energy-transport scenario: panel (a), energy profile for γ = π/4 and several values of time; panel (b), energy profile at the longest accessible time, t = 40, and several values of γ.Spin-transport scenario: panel (c), magnetization profile for γ = π/4 and several values of time; panel (d), magnetization profile at the longest accessible time, t = 40, and several values of γ.

FIG. 3 .
FIG. 3. Profiles of the energy and magnetization currents during the relaxation process.Energy-transport scenario: panel (a), energy current profile for γ = π/4 and several values of time.Spin-transport scenario: panel (b), spin current profile for γ = π/4 and several value of time.In both the panels (a) and (b), the dashed lines for negative/positive values of n/t are the results of the simulation of an homogeneous systems with the parameters of the left/right halves for t = 40, respectively.Internal light cone: panel (c), space-time plot of the left/right maxima n L/R of the energy current profile for several values of γ (when only one maximum is detected we set n L = n R = 0); panel (d), extrapolated velocity v R of the internal light cone computed by fitting the data of panel (c) for 30 ≤ t ≤ 40.The red line is obtained by fitting the data with v R = A ∆ + B which gives A = 1.33 ± 0.06 and B = −0.03± 0.03.

FIG. 4 .
FIG. 4. Extended chaotic junction.Top: cartoon of the setting where the two integrable Hamiltonians are joined via a chaotic Hamiltonian Ĥchaos on sites.Panels (a),(b): Analytic results from GHD for γ = π/4 in the energy-transport scenario for different values of the boundary resistance R, as indicated in the legend.Energy profile (a) and energy current (b).The black line shows the numerical data for t = 40.

FIG. 5 .
FIG. 5. Central equilibrium regions in the energy-transport scenario.Effective inverse temperatures: panel (a), profile of β eff n for γ = π/4 and several values of time; panel (b), profile of β eff n at the longest accessible time, t = 40, and several values of γ.Equilibration tendency: panel (c), distance between the three-site reduced density matrices ρn,n+2(t) and ρn,n+2[β eff n (t)] for several times and γ = π/4; panel (d), full time-evolution of the distance for four representative sites taken in both halves for γ = π/4.
is a good measure of the entanglement between the left and the right half at each time.Its mathematical expression reads S(t) = −Tr [ρ left ln(ρ left )].

FIG. 6 .
FIG. 6. Behaviour at the junction.Decay of the current at the junction: panel (a), energy-transport scenario; panel (b), spin-transport scenario.Time-dependence of the entanglement entropy of a system bipartition: panel (c), energytransport scenario; panel (d), spin-transport scenario.

FIG. 7 .
FIG. 7. Chemical potential profiles.Effective chemical potential: panel (a), profile of µ eff n for γ = π/4 and several values of time; panel (b), profile of µ eff n at the longest accessible time, t = 40, and several values of γ.
. A mismatch between the stationary currents computed with ρ TH and ρ ζ would imply a finite rate of growth of the charge Q µ (t; ζ) in the space region [−ζt, j * ], which however is always bounded by ||q µ || ∞ × ζ with || • || ∞ the operator norm and ζ arbitrarily small.As Eq. (C3) must hold for any µ, it imposes an infinite set of constraints on the GGE density matrices lim ζ→0 ± ρ ζ .