Front dynamics and entanglement in the XXZ chain with a gradient

We consider the XXZ spin chain with a magnetic field gradient and study the profiles of the magnetization as well as the entanglement entropy. For a slowly varying field it is shown that, by means of a local density approximation, the ground-state magnetization profile can be obtained with standard Bethe ansatz techniques. Furthermore, it is argued that the low-energy description of the theory is given by a Luttinger liquid with slowly varying parameters. This allows us to obtain a very good approximation of the entanglement profile using a recently introduced technique of conformal field theory in curved spacetime. Finally, the front dynamics is also studied after the gradient field has been switched off, following arguments of generalized hydrodynamics for integrable systems. While for the XX chain the hydrodynamic solution can be found analytically, the XXZ case appears to be more complicated and the magnetization profiles are recovered only around the edge of the front via an approximate numerical solution.


I. INTRODUCTION
Nonequilibrium dynamics of quantum many-body systems continues to be one of the most actively developing areas of condensed matter physics [1,2]. Within this vast field, the studies of integrable systems have received a particular attention [3]. To a certain degree, this is due to the powerful analytical tools of integrability, devised during the last century [4,5], which could be further extended to attack a variety of out-of-equilibrium problems. More importantly, integrable quantum systems turn out to be physically very interesting, as they show rather peculiar relaxational behaviour [6]. Last but not least, the spectacular development of cold-atom experiments [7], making the preparation, control and observation of near-integrable systems feasible [8], provided an ultimate boost to the theoretical research activities.
One standard approach to these nonequilibrium studies is a homogeneous global quench, where the initial state is prepared in the ground state of a certain Hamiltonian and time evolved with another one, differing in some global parameter, yet both of them being translational invariant [9]. Indeed, many essential features of anomalous relaxation and thermalization, resulting from the highly nontrivial families of conservation laws in integrable systems [10], could be understood in this simple context. Recently, however, the focus of interest has been shifted towards the study of inhomogeneous initial states, which might lead to the buildup of persistent macroscopic currents. An important milestone in this context has been the introduction of a generalized hydrodynamic picture [11,12], which extends the framework of classical hydrodynamics by taking into account the entire set of conservation laws in integrable systems. Starting from an inhomogeneous initial state, the method has been very successful in describing the time-evolved profiles of various observables (e.g. spin or energy densities and currents) in a hydrodynamic scaling regime, for a number of different situations and model systems [13][14][15][16][17][18][19][20][21].
The simplest example of an inhomogeneity in the context of spin chains is a domain wall, i.e. when two parts of a system are prepared initially in their (otherwise homogeneous) ground states with different magnetizations. The resulting front structures have been widely studied for the XX [22][23][24][25][26][27][28][29], the transverse Ising [30][31][32][33][34] and the XXZ chains [11,27,31,35,36]. In the latter case, a remarkable analytical solution has been found very recently [37], using the method of generalized hydrodynamics (GHD). Another commonly studied example of inhomogeneous initial states is the tensor product of two Gibbs states at different temperatures, which produces nontrivial energy density and current profiles under time evolution, see the recent reviews [38,39] and references therein.
In all of the above mentioned cases the inhomogeneity is sharp and there is no intrinsic length scale involved in the problem. Consequently, the time-evolved profiles of magnetization are scaling functions of x/t, i.e. the distance from the location of the initial domain wall divided by time, reflecting the ballistic nature of the transport. This is not true any more if the inhomogeneity is localised in an interface region of size ξ, with a smoothly varying magnetization profile. Such an initial state can be prepared by applying a magnetic field gradient along the chain, which is then turned off to monitor the subsequent time evolution [40][41][42][43]. In particular, the case of the gradient XX chain with a linearly varying field can be solved exactly [40], and the profiles turn out to be very closely related to the domain-wall problem.
Here we perform a similar study of the gradient XXZ chain, or equivalently, interacting spinless fermions in a linear chemical potential. Although the inhomogeneous potential breaks the integrability of the model, we shall show that, for small gradients, the ground-state magnetization profiles can still be perfectly captured by the Bethe ansatz method combined with a local density approximation. Moreover, it is demonstrated that the gradient also gives rise to nontrivial entanglement profiles. Relying on a recently introduced technique for the study of entanglement in inhomogeneous free-fermion ground states [44,45], we apply a perturbative extension of the method, which is shown to give good results for a broad regime of interactions in the gradient XXZ chain. Finally, we also monitor the time evolution of the profiles after switching off the gradient. Invoking the GHD method, an explicit analytical solution of the hydrodynamical density profile is obtained in the XX case. For the XXZ chain the method is applied with certain simplifications, that are shown to give a good numerical approximation in the edge region of the front.
The paper is structured as follows. In Sec. II we introduce the model and its basic features, continued by an analysis of the ground-state magnetization profiles in Sec. III. The entanglement profiles are studied in Sec. IV, with a brief introduction into the curved-space CFT technique for free fermions and its perturbative extension to the XXZ case. The front dynamics after switching off the gradient is investigated in Sec. V. We conclude the paper in Sec. VI by a discussion of our results against recent developments. Some details of the calculations are enclosed in two Appendices.

II. XXZ CHAIN WITH A GRADIENT
The Hamiltonian of the XXZ chain is given by where the spin-1/2 operators S α j = σ α j /2 on site j are given through Pauli matrices with α = x, y, z, and ∆ is the interaction strength. The XXZ chain with a gradient is obtained by adding a linearly varying magnetic field with the slope set by δh. The gradient field is shifted such that it has a zero between sites j = 0 and j = 1, and thus the Hamiltonian is invariant under the spininversion S z j → −S z 1−j . Note that, due to the special geometry of the problem, it is natural to consider open boundary conditions. Performing a Jordan-Wigner transformation, the gradient XXZ chain is equivalent to a fermionic hopping chain with nearest-neighbour interactions and a linear chemical potential.
Our first goal is to characterize the inhomogeneous ground state of H gr via the magnetization and entanglement profiles. The inverse of the gradient defines a length scale ξ = δh −1 . For |j| ≫ ξ, the magnetization is saturated at the values S z j = ±1/2 on the far rightand left-hand sides of the chain, whereas the two parts are connected by a nontrivial interface region. Although H XXZ is solvable by Bethe Ansatz, the gradient field breaks the integrability of the Hamiltonian. Nevertheless, we shall try to capture the properties of the interface by using standard Bethe Ansatz methods in combination with a local density approximation (LDA). The essential argument behind LDA is that, if the field varies slowly enough, δh ≪ 1, then the properties of the inhomogeneous ground state will also change smoothly. In particular, at site j the spins experience a local field strength h(j) = δh(j − 1/2), and the system will locally occupy the ground state corresponding to a homogeneous field h = h(j). In other words, if the length scale ξ associated with the gradient is much larger then the lattice spacing, LDA should give a good approximation of the profiles.

III. MAGNETIZATION PROFILES
Obtaining the magnetization profiles for the gradient problem thus boils down to solving for the ground state of the XXZ chain in a nonzero field h. We will focus on the planar regime |∆| < 1 and use the standard parametrization ∆ = cos γ. Working directly in the limit L → ∞, this is a well-known problem in terms of the thermodynamic Bethe Ansatz (TBA) [46,47]. Indeed, the eigenstates of the XXZ chain are parametrized by rapidities λ, corresponding to magnon excitations created upon the fully polarized state. In particular, the ground state has a simple Fermi sea character, with real rapidities fully occupying an interval λ ∈ [−Λ, Λ]. While for finite L the rapidities satisfy some appropriate quantization conditions, in the thermodynamic limit the roots become continuous, and their density ρ(λ) follows from the TBA equation [46,47] where the derivative of the bare momentum k ′ (λ) ≡ θ ′ 1 (λ) and the integral kernel K(λ) ≡ θ ′ 2 (λ) are given through Since the magnons are created on the ferromagnetic state, the local magnetization can be obtained as Hence, for a fixed Fermi rapidity Λ, the magnetization is a simple function of the root density, which in turn follows from the integral equation (3). However, to extract the profiles one needs the magnetization as a function of the magnetic field h. To obtain the relation h(Λ), one first defines the so-called dressed charge function Z(λ), which is the solution of yet another integral equation In terms of the root density and the dressed charge, the magnetic field is then given by [46,47] h(Λ) = 2πρ(Λ) Z(Λ) sin γ . Unfortunately, the analytical solutions of ρ(λ) and Z(λ) are available only in the limiting cases Λ = ∞, corresponding to S z = h = 0, and Λ = 0 which yields S z = 1/2. In the latter case the value of the critical magnetic field h c follows from the trivial solutions 2πρ(0) = θ ′ 1 (0) and Z(0) = 1 as If the magnetic field exceeds the value h c , the ground state remains fully polarized. Consequently, within the validity of the LDA framework for the gradient chain, the interface connecting the ferromagnetic regions is located between −h c < h(j) < h c , and its half-width is given by ξ(1 + ∆). The complete profile can be obtained via the numerical solutions of the TBA equations (3) and (6), by gradually increasing Λ and calculating the matching pairs of S z and h from (5) and (7). The LDA profile can then be compared to the numerical value of S z j in a finite size chain, obtained via density matrix renormalization group (DMRG) simulations [48,49]. The discarded weight was set to 10 −12 without any restriction on the maximal bond dimension, and we required the groundstate energy to converge within ∆E < 10 −7 for three consecutive DMRG sweeps. The results of these calculations are shown on Fig. 1, for a size L = 200 and field gradient δh = 0.02. One can see that, up to some oscillations, the LDA method gives a perfect description of the profile, both in the antiferromagnetic (∆ > 0, left) as well as in the ferromagnetic (∆ < 0, right) regimes. In the fermionic language, the repulsive or attractive nature of the interactions yields a stretched or squeezed profile as compared to the non-interacting case where the halfwidth of the interface is given by ξ [40]. Interestingly, the match between the LDA profiles and the data remains very good even for ∆ = −0.8, where the interface is rather sharp. On the other hand, the oscillations around LDA are more pronounced for large ∆ > 0.
Finally, one observes in Fig. 1 that the profiles show a nonanalytic behaviour at the edge of the interface. This can be captured by using the approximate solutions of the integral equations (3) and (6) in the limit Λ ≪ 1, which yields for the magnetic field [47] h ≈ h c − 1 tan 2 γ/2 Furthermore, the integral in Eq. (5) can also be approximated as Λθ ′ 1 (0)/π such that the magnetization becomes Substituting h → h(j) = (j −1/2)/ξ, one can see that the profile S z j shows a square-root singularity at the edge, which is independent of the interaction parameter ∆.

IV. ENTANGLEMENT PROFILES
We continue by analyzing the entanglement profiles in the ground state. We are primarily interested in the behaviour of the von Neumann entropy where ρ A is the reduced density matrix for a bipartition of the ground state |ψ . We consider only the case of intervals A = [−L/2 + 1, r] and B = [r + 1, L/2], where r measures the distance of the cut from the centre of the chain. In particular, r = 0 corresponds to the half-chain case and the entropies for ±r must be exactly equal, due to the symmetry of the ground state under the combined action of reflection and spin-flip transformations.
Although the entanglement entropy is obtained straightforwardly from DMRG simulations, it is much harder to develop an analytical method, based on some LDA argument, that could provide accurate predictions for the entropy profile. A recent breakthrough has been achieved for free-fermion systems by realizing that, in a proper continuum limit, the inhomogeneous problem can be mapped onto a CFT in a curved space [44]. Our goal is to use this approach as a starting point for the more complicated interacting problem of the gradient XXZ chain, we shall thus first review the main ideas and arguments of the method.

A. Dirac field theory in curved space
Let us consider a free-fermion Hamiltonian in 1D continuous space with a slowly varying potential term V (x) where c † (x) and c(x) are the creation/annihilation operators. With V (x) ≡ 0, the constant chemical potential µ sets the ground state to be a Fermi sea between the Fermi momenta ±k F , where from the dispersion ε(k) = k 2 one has k F = √ µ. Within the realms of LDA, the main effect of a slowly varying potential Introducing the time-evolved fermionic operator c(x, y) = e yH c(x)e −yH with imaginary time y = it, a more precise description of the local spacetime behaviour of the theory is encoded in the fermionic propagator where we introduced Here the Fermi velocity is defined as v F (x) = ε ′ (k F (x)) and its space dependence is entirely due to the variation of the Fermi momentum k F (x). Note that the form (13) of the propagator is valid only in the local neighbourhood of spacetime point (x, y), as the derivation involves a linearization of the dispersion ε(k) around k F (x). The main idea of Ref. [44] is to view the two terms of Eq. (13) as the chiral propagators of a 2D massless Dirac field. In a flat Euclidean spacetime, this is defined by the action where z = x + iv F y and the right-and left-moving chiral propagators are given by Comparing to (13), one indeed finds the expected local behaviour, up to some phase factors which can be incorporated by a chiral gauge transformation. The crucial question is, however, whether there exists a proper choice of complex coordinates w(x, y), such that it provides a globally valid Dirac action in a curved space? It turns out that for this to be the case, the underlying Riemannian metric should satisfy [44] so that the chiral propagator behaves locally as Hence, to be able to derive Eq. (13) from a Dirac theory in curved space, one should have It is easy to see that this condition is satisfied by the choice Note that this map involves only the transformation of spatial coordinates, whereas the imaginary time is left untouched. Moreover, it depends on a single parameter, namely the local Fermi velocity v F (x).

B. Inhomogeneous Luttinger liquids
The case of the XXZ chain, or equivalently interacting fermions, turns out to be more complicated. Indeed, the low-energy behaviour of the model is described, after standard bosonization procedure [50], by a Luttinger liquid which is an effective field theory involving the bosonic field φ and its dual θ, satisfying the standard bosonic Clearly, apart from the Fermi velocity v F , we have now the extra Luttinger parameter K that enters the field theory description. It should be noted that, while bosonization yields both of these parameters only within perturbation theory, the exact values of v F and K can be fixed from the Bethe ansatz solution. This procedure is wellknown for the XXZ chain with a homogeneous magnetic field [47,50,51], with the results shown in Fig. 2. The method is summarized in Appendix A. Clearly, when considering the gradient XXZ chain, the variation of the parameters v F (h) and K(h) w.r.t. the magnetic field translates into a spatial variation, with the identification h → x = j/ξ (from here on, the −1/2 shift in the j coordinate will be suppressed for brevity  of notation). Hence, the field-theory description is still given by a Luttinger liquid as in Eq. (21), however with parameters v F (x) and K(x) that vary slowly with position, on a length scale given by ξ. The situation is thus more complicated than the one considered recently in Ref. [45], where the case of variable Fermi velocity v F (x) but uniform Luttinger parameter K was dealt with. This might occur in inhomogeneous problems such as the XXZ chain with a varying coupling J(x). Indeed, in this case the parameter K can essentially be scaled out from the Luttinger liquid Hamiltonian (21) by a canonical transformation φ → √ Kφ and θ → θ/ √ K. On the level of CFT, the rescaling of the bosonic fields amounts to changing the compactification radius of the theory. Nevertheless, the variation of v F (x) can still be incorporated in a curved metric just the same way as was done in (20) for the Dirac action [45].
The inhomogeneous Luttinger liquid problem, with both parameters v F (x) and K(x) being nonuniform, is hard in general. Instead of trying to find an exact solution, we shall rather be interested in uncovering a perturbative regime, where the Dirac field theory in curved space yields a good approximation. Having a look at the behaviour of K(h) in the right of Fig. 2, one observes that the variation is rather strong already for relatively small interactions ∆. At first sight, this would opt against a perturbative approach. However, to better understand the role played by the Luttinger parameter, one should have a look at the fermionic propagator. From standard results of bosonization [50], the propagator of the rightmoving fermion has the form where the proportionality sign indicates that we have omitted the phase factor and normalization constant, cf. Eq. (13). Similarly, the left-moving propagator is obtained by substituting δz → δz.
The local deviation from a Dirac theory is thus encoded in the scaling dimension Obviously, setting K(x) ≡ 1 we obtain D(x) ≡ 0 and recover the local Dirac propagators in (16). Moreover, the result (22) suggests that, whenever the scaling dimension varies slowly and satisfies D(x) ≪ 1 in the full spatial domain, one should expect the curved-space Dirac theory to yield a good approximation to the inhomogeneous Luttinger liquid problem. Remarkably, as shown in the inset on the right of Fig. 2, this condition is indeed satisfied for a broad range of ∆. In particular, one has D → 0 as h → h c for arbitrary values of ∆ [52]. In general, D(x) increases towards the bulk and reaches a maximum deviation of approx. 4% for ∆ = 0.5, resp. 8% for ∆ = −0.5.
Clearly, deviations increase further as |∆| → 1 and one leaves the perturbative regime.

C. Entropy and curved-space CFT
The validity of the above arguments will be tested by comparing the predictions of the curved-space CFT approach for the entanglement entropy to the DMRG results. Once the curved metric is fixed by (20), the entropy can be extracted by calculating expectation values of twist-field operators [53] and applying conformal transformations to simplify the CFT geometry [44]. The steps of this procedure are collected in Appendix B. It should be noted that, while there is no exact analytical expression for v F (h), we found that, for intermediate values of ∆, the behaviour of the Fermi velocity as a function of the magnetic field is very well approximated by the ansatz  (20) and following the recipe in Appendix B, one finds that the entropy of a segment [−L/2 + 1, r] is given by where the conformal distance reads Unfortunately, the curved-space CFT calculation fixes the entropy only up to some non-universal (i.e. depending on the underlying lattice model) term C(r, ∆), which still contains a dependence on the position r of the cut. In fact, this is completely analogous to standard CFT calculations of the ground-state entropy in the XXZ chain with a homogeneous field h. There the entropy of a segment of length ℓ has the form with the universal term given by the chord length From analogy of the expressions (25) and (27), and in spirit of the LDA argument, one infers that the nonuniversal terms are related as In general, the non-universal term must be extracted from the lattice model and an analytical expression is known only in the free-fermion case [54,55] with a numerical value of the constant c 0 ≈ 0.4785. This yields the entropy of the gradient XX chain as in perfect accordance with the results found in [56]. For ∆ = 0, the termC(h, ∆) must be extracted by fitting the ground-state entropy with the ansatz (27), with the results shown in Appendix B. Finally, having obtained both the universal and nonuniversal terms, one can check the results agains the DMRG data. This is shown in Fig. 3, plotted against the variable r/ξ, with the red lines showing the interpolated curved-space CFT ansatz (25). The dots on the lines indicate the parameters h = r/ξ, where the actual values of C(r, ∆) were determined using the relation (29) and data fits for the homogeneous chain. In general, the ansatz seems to give a very good quantitative description of the entropy, although some deviations are visible for larger ∆ values. In particular, one observes that the match between (25) and the DMRG data is best around the edge of the profile, whereas deviations increase as one moves towards the bulk. This completely parallels the behaviour of the scaling dimension D(x) (see inset of Fig. 2) and gives a further confirmation of our perturbative approach. Note that we have also checked the result using the exact form of v F (h) instead of the approximation (24). In this case the conformal distance is given by (B11) in Appendix B and has to be evaluated numerically. The corresponding deviations induced in (25) are of the order of the fitting error in C(r, ∆) and thus the approximate form (26) of the conformal distance is perfectly justified.

V. FRONT DYNAMICS
We now proceed to investigate the dynamics that follows after switching off the gradient magnetic field. The time-evolved state of the system is then given by As for the static case, we shall study the time-evolved profiles of magnetization and entanglement. The problem has a single characteristic length scale ξ, it is thus natural to work with the rescaled coordinate x = j/ξ and time τ = t/ξ. Since the time-evolution of the state (32) is governed by the integrable XXZ Hamiltonian, we expect the behaviour of local observables to be given via a generalized hydrodynamic (GHD) picture introduced recently [11,12]. Within the GHD framework the effective hydrodynamic state of the system is described by spaceand time-dependent occupation functions of rapidities, that are subject to certain advection equations. Before formulating the problem in general, it is very instructive to start with the noninteracting XX case, where an explicit solution of the GHD equation can be found and compared to the exact lattice solution [40].

A. Exact solution for the XX chain
The XX Hamiltonian in fermion language is given by and can be diagonalized by a Fourier transform. Hence, due to the free-fermion nature of the problem, one can simply use the occupation function n(x, τ ; k) in momentum-space, instead of working with rapidities. The GHD equation [11,12] for the occupation function then reads where the single-particle velocity is v(k) = sin k. We assume that the hydrodynamic state of the system is described, throughout the time-evolution, by an occupation function which is a local Fermi sea between two Fermi points k ± (x, τ ) for any time τ . Since (35) depends on x and τ only via the location of the Fermi sea, the advection equation (34) yields the equation of motion for the Fermi points One is thus left with two independent advection equations which can be solved by the method of characteristics. This gives k ± implicitly, as the solution of the equation where the function F ± describes the initial conditions Note that (38) is simply the LDA result for the gradient XX chain, where we considered δh < 0 in Eq.
(2) such that the fermionic density c † j c j = S z j + 1/2 is zero (one) on the far right (left), and thus the current flows from left to right. Inverting Eq. (37), one should look for the solution of which is a simple quadratic equation with roots Note that real solutions exist only for |x| ≤ x max and thus x max = √ 1 + τ 2 gives the location of the front edge. Using trigonometric identities, the solution can be rewritten in the simple form Comparing to the initial data (38), one sees that, on top of an overall x-independent drift term, the solution simply gets rescaled by x max . With the solution (41) at hand, it is easy to evaluate the expectation values of the particle density ρ and current J in the hydrodynamic state (35). A simple calculation gives Substituting the solution (41) for the Fermi points, one arrives at Note that the drift term in (41) does not enter the expression of the density, which thus has the same functional form as at t = 0, up to a rescaling of the distance by x max . It is easy to check that the particle density and current in (43) satisfy the usual continuity equation The hydrodynamic result (43) can now be compared to the exact lattice solution of the front dynamics. The lattice density and current are given via the fermionic correlators C j,l (t) = c † j (t)c l (t) as ρ(j, t) = C j,j (t) , J (j, t) = Im C j,j+1 (t) .
The explicit form of the correlation matrix was obtained in Ref. [40] and reads where ϕ = atan(t/ξ) and (47) is the ground-state correlation matrix with Bessel functions J j (ξ). Thus, up to a phase factor, the correlators have the exact same form as in the ground state, with an effective length scale ξ 2 + t 2 . In fact, this is the simplest manifestation of a more general phenomenon of an emergent eigenstate solution, discussed in Ref. [43]. Furthermore, the solution can also be related to the case, where the initial state is given by a sharp domain wall (i.e. the initial density is a step function) by setting ξ = 0 in (46). This yields ϕ = π/2 and the initial correlations (47) have to be evaluated with an argument t. For the domain-wall case the density and current are well-known scaling functions of the variable v = j/t and read [22] Then, in the general case ξ = 0, the scaling form of ρ simply follows by substituting t → ξ 2 + t 2 , since the density is unaffected by the phase ϕ. For the current J one has to multiply the result by an additional factor of sin ϕ, as is clear from (45) and (46). These manipulations lead immediately to the hydrodynamic result (43), which is thus exact in the appropriate scaling limit.

B. Numerical results for ∆ = 0
We start the discussion of the interacting case by summarizing the GHD method for a general integrable system. These models possess an extensive set of local (or quasi-local) conserved charges [10], each of them yielding an appropriate continuity equation that must be respected by the dynamics. Furthermore, within the TBA description, each of the charge and respective current densities are fixed by occupation functions n l (λ) where λ is the rapidity parameter and l indexes the families of quasi-particles in the integrable model at hand [57,58]. In fact, as it was shown in the pioneering works [11,12], these continuity equations can be reformulated as advective flow equations for the occupation functions Although formally very similar to the free-fermion case Eq. (34), the essential difficulty of the above advection equation is due to the term v l ({n}; λ), which is the dressed quasi-particle velocity. Indeed, as the index {n} suggests, it depends in general on the full hydrodynamic state of the system, therefore coupling the complete set of advection equations (49). The precise form of the dressed velocity is fixed within the TBA formalism for the corresponding model via the dressing equation, see Refs. [11,12] for details.
In the present case, the ground-state results of Sec. III and a naive analogy to the XX case in (35) would suggest that the hydrodynamical description of the dynamics could be specified by a single occupation function of the form where the entire spacetime-dependence is via the Fermi rapidities Λ ± (x, τ ). There is, however, an essential difference between the momentum and rapidity parametrization, which is already transparent in the ∆ = 0 case. Namely, the real rapidities yield only the low-momentum particles k ≤ π/2, whereas those with k > π/2 correspond to complex rapidities along the iπ axis [46]. In TBA language, the latter correspond to a different family of quasi-particles with negative parity, and thus the modes are not smoothly connected to each other. In fact, this is also apparent for the ground-state at τ = 0, where for x → 0 one has Λ ± → ±∞ and the sea of rapidities in (50) is already completely filled. While in the XX case these complications can be avoided by working in momentum space, for ∆ = 0 one is bound to run into difficulties when considering the solution of the advection equations (49). In general, the quasi-particle structure can be constructed for some particular values of γ (i.e. rational multiples of π) and is given in terms of string-solutions for the rapidities [46]. In analogy to the XX case, the real rapidities yield only quasi-particles with bare momentum restricted to the interval k ∈ [−γ, γ], and the particle content at τ = 0 and x < 0 (i.e. in the negative magnetization sector) must be obtained by looking at the nontrivial action of spin-flip on the quasi-particle structure [59,60]. However, since the resulting occupation functions n l (x, τ ; λ) are all coupled together via the dressed velocities, even the numerical solution of (49) becomes rather complicated.
Instead of trying to find a full numerical solution to the advection equations, we shall try to simplify the problem by dealing with only one particle type of real rapidities (l = 1), with occupation function as in Eq. (50). Obviously, such an ansatz will immediately fail around x → 0 + for any τ > 0, as the quasi-particles from the left penetrate the right half of the system. Nevertheless, it is reasonable to expect that (50) gives a good approximation of the hydrodynamics around the front edge, at least for small times τ , where all the other quasi-particle occupations could be neglected. The equation of motion for the Fermi rapidities then reads where the dressed quasi-particle velocity is defined via the dressing equation Note that the bare momentum derivative k ′ (λ) = θ ′ 1 (λ) is given by Eq. (4), whereas the bare energy ǫ(λ) is given by (A2) of Appendix A.
In a restricted spatial domain, an approximate solution for Λ ± can be found again via the method of characteristics, by solving with the initial condition F ± (x) = Λ ± (x, 0). Inverting this function, we get the initial position x(Λ) = x(−Λ) = F −1 + (Λ) as a function of the Fermi rapidity Λ ≥ 0. With the identification x → h, this is nothing else but the function in Eq. (7), we are thus looking for the solutions of Since the left hand side of (55) is strictly positive, two distinct solutions exist only in an interval x min < x < x max for any fixed τ > 0. In particular, the front position x max (τ ) denotes the coordinate where the solutions Λ − = Λ + = Λ * of (55) coalesce, whereas x min (τ ) corresponds to Λ + → ∞. In general, the solutions Λ ± can be found numerically by an iterative procedure and are shown on Fig. 4, with an upward drift of the Fermi rapidities for increasing times.
With the solutions Λ ± (x, τ ) at hand, it remains to extract the magnetization of the hydrodynamic state with occupation function (50). Analogously to the groundstate in Eq. (5), it is now given by where x = j/ξ and the root densityρ(λ) follows from an integral equation similar to Eq. magnetization profile (red lines) show a very good agreement with the DMRG data. For larger times, however, there are visible deviations which tend to increase for larger ∆. The deviations from the approximate solutions might have different sources. First of all, we have completely neglected the contributions from the quasi-particles l = 1 emanating from the left half of the chain. Indeed, since the evolution of their occupation functions is also advective, the spatial domain where n l (x, τ ; λ) = 0 grows for large τ and eventually penetrates even into the edge region. For small enough τ , however, due to the boundedness of the corresponding speed of advection v l , these particles should not contribute to the magnetization around the edge. Nonetheless, in case of an overlap region of nonvanishing occupations l = 1 around x min , they do couple the advection equations (51) via the quasi-particle velocity. The coupling introduces small changes into the solution Λ ± (x, t) far away from x max , which could nevertheless propagate into the edge region for large enough times. In fact, as remarked in [18,19], the proper way to integrate (49) would be to apply the method of characteristics in small backwards time steps always updating the velocities v l ({n}; λ) using the proper TBA prescription [18,19] that takes into account all the occupations {n}, evaluated at x and τ . The effects discussed above can also be visualized by looking at the entropy profiles, shown in the bottom row of Fig. 5. In fact, there is a clear separation between the edge region, where the profile resembles that of the ground state (see Fig. 3), and a bulk part with an enhanced amount of entropy. The latter should correspond to the domain where the quasi-particle excitations from the left half have already arrived. As time increases, the edge region gradually shrinks and the separation from the bulk gets more pronounced for larger ∆. Moreover, for large τ and ∆ = 0.5 one has additional nonmonotonous structures in the profile. We conclude by providing an approximate formula for the front size x max in the small time limit τ ≪ 1. As remarked earlier, this is determined by the condition Λ − = Λ + = Λ * , such that the magnetization from (56) equals 1/2. Since the occupation is identically zero, one has no dressing and the velocity is given by where in the second step we used that Λ * ≪ 1 for small times. Using the approximate form of the ground-state magnetization curve in (9), the equation (55) becomes quadratic and can be solved for Λ * . Requiring the solution to be unique yields the relation whereas the leading x → x max behaviour of the magnetization is given by The front thus expands very slowly for early times, in complete analogy with the XX results is (43). In contrast, for times τ ≫ 1 much larger than the inherent length scale of the problem, one expects the usual ballistic behaviour x max ∝ τ .

VI. DISCUSSION
We have studied the magnetization and entanglement profiles for the inhomogeneous ground state of the XXZ chain with a field gradient. For the description of the magnetization profiles we combined the LDA argument with the TBA solution in a homogeneous field, and found perfect agreement for small gradients. The entanglement profiles are well approximated, for a wide range of interactions ∆, by a technique relying on a mapping to a Dirac field theory in curved space. The dynamics of the interface was also investigated after switching off the field gradient, and compared to the predictions of generalized hydrodynamics. For the XX chain the GHD equations can be solved analytically and reproduce the exact lattice solution in a proper scaling regime. In the XXZ case we have considered only an approximate form of the GHD equations which was argued to be applicable in the edge regime of the front. The numerical solution of these equations were indeed found to give a good description of the edge profile for small enough times.
It is very interesting to compare our results to a recent study of front propagation in the XXZ chain from a domain-wall initial state [37]. Due to the simple structure and relations between the quasi-particle occupations in the left-and right-hand side, an analytical solution to the GHD equations, and hence the corresponding magnetization and current profiles could be found [37]. In particular, it has been shown that, for any ∆ = 0, the magnetization profile around the front edge changes from a square-root to a linear behaviour. This is in contrast with our results for the front dynamics from an initial state with a smooth gradient, where the edge of the profiles still show a characteristic square-root singularity. For very short times τ ≪ 1, this behaviour is also supported by the approximate solution (60) via the GHD equations. Note, however, that for larger times there is a strong indication of a shrinking edge regime (see Fig. 5) from the DMRG results. Thus it seems plausible, that for very large times the above regime is entirely washed out and replaced by a different edge behaviour. To check this one should push the simulations further in time, which requires considerable efforts.
It would also be desirable to obtain a complete numerical solution of the GHD equations (49), instead of the ansatz (50) which is only applicable around the front edge for small enough times. For an initial state with a slowly varying temperature profile [18,19], the integration of the advection equations was achieved by a backwards Euler method, Eq. (57), which is essentially the method of characteristics applied repeatedly with an infinitesimal time step. This could be generalized to the present case by taking into account the proper initial conditions n l (x, 0; λ) for the full set of occupation functions. The question is, however, if the occupation functions characterized by a single Fermi sea would be stable under such an advective evolution. In fact, in a recent study of the Lieb-Liniger model with an inhomogeneous initial density profile [17], it has been pointed out that instabilities are likely to occur, leading to a breakup of the Fermi sea into several disconnected parts. Whether this would also occur for the gradient XXZ chain at hand is clearly a very interesting open question which deserves further studies.
Another open issue is the change in the transport properties as one approaches ∆ → 1, i.e. the boundary of the gapless regime. Indeed, for the evolution from a sharp domain wall, it was observed in earlier numerical studies that the ballistic behaviour is replaced by a slower, but still super-diffusive transport [35]. Although the superdiffusivity seemed to be confirmed by further numerical studies [61,62], very recently it has been argued that the dynamics of the domain-wall melting could also be compatible with simple diffusion [63,64]. Studying the melting of a smooth interface in the gradient XXX chain could shed further light to this question.
Regarding the entanglement profiles, it should be stressed that our curved-space CFT arguments for the calculation of the ground-state entropy apply only in a perturbative regime. However, we identified that the relevant small parameter D(x) ≪ 1 is the scaling dimension in Eq. (23), whereas the Luttinger parameter itself might still show considerable deviations from the free-fermion point K = 1. Nevertheless, it would be of great interest to devise a method which could incorporate an arbitrary smooth variation of the Luttinger parameter K(x) (i.e. the compactification radius of the CFT). Note, however, that even if this problem could be tackled, a complete analytical solution for the entropy profile seems out of reach for ∆ = 0, as it would require knowledge of the non-universal constant term in (27) for a homogeneous XXZ chain. While this term is easy to extract from data fits (see Fig. 6 in Appendix B), we are not aware of any analytical approach to this problem.
Finally, understanding the dynamical entropy profiles also remains an open question. Our numerical results make it clear that there is strong qualitative difference between the XX and XXZ chains in this respect. Indeed, in the non-interacting case one deals essentially with a ground-state problem, with the full entropy profile given as in Eq. (31) under substitution of an effective length scale ξ 2 + t 2 . In contrast, for ∆ = 0 there is a clear separation between the bulk and edge regimes. While the latter still has the qualitative features of a ground-state profile, the bulk shows a more rapid growth of entanglement, with the separation becoming more pronounced for larger ∆. In view of the recent analytic results on the entropy evolution in XXZ chains [65,66], the question naturally emerges whether similar arguments could be applied for the present case.
The Fermi velocity is defined as the derivative of the dressed energy w.r.t. the dressed momentum evaluated at the Fermi rapidity. Introducing v dr = (ǫ ′ ) dr with the prime denoting derivative w.r.t. rapidity, the dressed velocity follows from the integral equation Moreover, noticing that the derivative of the bare momentum is given by k ′ (λ) = θ ′ 1 (λ) as defined in Eq. (4), it follows from (3) that the derivative of the dressed momentum is simply the root density multiplied by a factor of 2π. Hence, the Fermi velocity can be obtained as Finally, one should also fix the Luttinger parameter. The easiest way is to consider the magnetic susceptibility which can be given as [51] Together with (5) and (A6), this could be used as the definition of K. However, using some variational properties of the Bethe ansatz root density [47], it is also possible to show that the Luttinger parameter is simply related to the dressed charge as The formulas (A6) and (A8) can now be evaluated, for a fixed value of Λ, via the numerical solution of integral equations (3), (6) and (A5). Varying the value of Λ, one obtains the Fermi velocity v F (Λ) and Luttinger parameter K(Λ) as functions of the Fermi rapidity. The functions v F (h) and K(h) can then be found by inverting the function h(Λ) in (7). The result is shown in Fig. 2 for various values of the interaction parameter ∆. The Fermi velocity (shown on the left) is found to be very well approximated by the square root expression (24) reported in the main text, although deviations increase with increasing |∆|. Finally, the expectation value on the upper half plane is simply given by T n (g(w 0 )) uhp ∝ (Im g(w 0 )) −∆n , where the proportionality sign indicates the presence of non-universal multiplicative factors.
Putting everything together, we arrive at the following expression for the Rényi entropy S n = c 12 (1 + 1/n) ln L + C n , where the conformal distance is given by Im g(w 0 ) .
Substituting e σ = v F and using (B6), this can be rewritten as Now the term Re w 0 is given by an integral expression in terms of the Fermi velocity, for which we do not have an exact analytical expression. However, we have found that the simple square-root formula in Eq. (24) gives a very good approximation of the actual data obtained from TBA. Using this approximation, the integral can be carried out and one finds Inserting this into (B11), we obtain the conformal distance as Using R = h c ξ, we arrive at the formula (26) reported in the main text. Unfortunately, due to the presence of the non-universal term C n , the expression in (B9) does not fix the full analytical form of the Rényi entropies. In fact, this term still depends on both the location of the cut r, as well as the interaction parameter ∆. It is, however, possible to fix C n from ground-state calculations of the entropy on the lattice. Indeed, considering a finite XXZ chain of length L in a homogeneous magnetic field h < h c , the ground state corresponds to a c = 1 CFT on a strip. Using the mapping (B6) with the substitution W → L/2, the CFT result for the entropy of a segment of length ℓ reads S n = c 12 (1 + 1/n) lnL +C n , where the conformal distance is given by the chord length in Eq. (28). Hence, setting the ground-state magnetic field equal to h = r/ξ, one finds from analogy of Eqs.
The variation of C n as a function of the cut position can thus be inferred from the ground-state behaviour C n as a function of the magnetic field. The latter can be obtained by analyzing the lattice data obtained from DMRG for various ℓ ≤ L/2, with fixed h, ∆ and L, and fitting the ansatz (B14) for each data set. Repeating the procedure for various h and ∆, we obtained the nonuniversal term for the n = 1 case. This is shown on Fig.  6, as a function of the scaled magnetic field h/h c and various values of ∆. It should be noted that an analytical expression ofC 1 is known only for the free-fermion case ∆ = 0, where it is given by Eq. (30). Numerical results for the case h = 0 were also obtained in Ref. [67].