Exact hydrodynamic solution of a double domain wall melting in the spin-1/2 XXZ model

We investigate the non-equilibrium dynamics of a one-dimensional spin-1/2 XXZ model at zero-temperature in the regime $|\Delta|<1$, initially prepared in a product state with two domain walls i.e, $|\downarrow\dots\downarrow\uparrow\dots\uparrow\downarrow\dots\downarrow\rangle$. At early times, the two domain walls evolve independently and only after a calculable time a non-trivial interplay between the two emerges and results in the occurrence of a split Fermi sea. For $\Delta=0$, we derive exact asymptotic results for the magnetization and the spin current by using a semi-classical Wigner function approach, and we exactly determine the spreading of entanglement entropy exploiting the recently developed tools of quantum fluctuating hydrodynamics. In the interacting case, we analytically solve the Generalized Hydrodynamics equation providing exact expressions for the conserved quantities. We display some numerical results for the entanglement entropy also in the interacting case and we propose a conjecture for its asymptotic value.


Introduction
A domain wall (DW) state of a quantum spin-1/2 chain is a product state prepared by joining two domains of aligned spins with opposite magnetization, as for instance |↑ . . .↑↓ . . .↓〉.Despite its simplicity, the time evolution of such a state shows non-trivial features of non-equilibrium dynamics and, for this reason, it has been the main character of a large number of studies, including stability analysis [1][2][3], exact computations for the free case e.g.[4][5][6][7][8][9][10][11][12][13][14][15], approximate and numerical results for the interacting integrable e.g.[16][17][18][19][20][21][22][23][24][25][26][27][28][29] and non-integrable e.g.[30][31][32][33][34] chain.For some integrable spin chains such as XXZ with anisotropy parameter |∆| < 1 (see Eq. ( 1) below), this research has led, only very recently, to an asymptotic analytical understanding of the DW dynamics e.g.[35][36][37][38][39][40][41], obtained by means of Generalized Hydrodynamics (GHD) [35,42] and of its quantum fluctuations [43].Though, some aspects of the melting process still lack of fully analytical explanations, such as the subleading corrections to the front dynamics [45] or the diffusive behavior at |∆| ≥ 1 [46][47][48][49].This exact solution substantiates an intuitive hydrodynamic picture at ballistic scales of the non-equilibrium dynamics which follows.Because of the integrability of the XXZ spin chain, stable quasiparticles excitations are initially emitted at the domain wall and propagate in time with a constant velocity, whose value depends on the interaction coupling.The fastest excitations of the spectrum define a fan-shaped spatio-temporal region (called lightcone) around the junction, inside which entanglement spreads and the initially ordered domains melt.Quite interestingly, the physics of the melting process is not modified by the presence of interactions.In particular, it has been shown that, up to a rescaling of the lightcone, the charges profiles in the gapless XXZ model [40] have the same behavior as a free Fermi system [4].This similarity extends even to the half-system entanglement entropy, which shows the same asymptotic growth as S 1 (0, t) ∼ 1  6 log t in the free [7] and interacting [41] case.Notice that this is a peculiarity of DW states, related to the structure of its Bethe Ansatz equations.Indeed, in generic integrable models, interactions are responsible for a dressing of the observables that typically modifies the non-equilibrium dynamics, see e.g.Ref. [44,[50][51][52] for recent reviews.The case of bipartite spin states with initial correlations has also been considered in literature [6,20,21,53], where it is observed that the presence of initial entangled quasiparticles results into a faster spreading of the half-system entanglement.
In this work, we extend the rich analysis on DW states to the case where three domains of aligned spins with different orientations are joined together, i.e., for an initial state like |↓ . . .↓↑ . . .↑↓ . . .↓〉.In the following, we will refer to this state as double domain wall.Intuitively, the non-equilibrium dynamics of a XXZ spin chain prepared in a double DW state can be split into two regimes.At early times, the light cones generated at the two junctions do not intersect and the evolution is that of two independent DW states of Ref. [40].Conversely, when the quasiparticles emitted from the two walls meet, a non-trivial interplay between the physics of the single DWs takes place and characterizes the melting process at later times.In this regime, the half-system entanglement (i.e., with the entangling point in the middle of the central domain, which is initially equal to zero) is fed by the double contribution of left and right moving quasiparticles coming from the two walls and therefore it undergoes a rapid growth.At large times, the entanglement saturates to a constant value S 1 (0, t) ∼ 1  3 log ℓ, which depends on the size 2ℓ of the domain of upwards spins at t = 0. Throughout the rest of the work, we present exact calculations and numerical analysis in order to corroborate these intuitions with quantitative arguments.
Outline.We organize the contents of the paper as follows.In Sec. 2, we introduce the model and we set up the quench protocol considered in this work.Although the physics of the melting process is found to be qualitatively similar for any value of the interaction coupling |∆| < 1, we treat the non-interacting case (∆ = 0) in a dedicated section (Sec.3), for a clearer exposition.In particular, in Sec.3.1, we characterize the semi-classical evolution in phase space in terms of the Wigner function and we derive exact asymptotic results for the magnetization and the spin current profiles.In Sec.3.2 we discuss the behavior of the entanglement.As pointed out in the recent works [6,41,43,54], this study requires a re-quantization of the semi-classical hydrodynamic background in terms of a Luttinger liquid and the use of conformal invariance.Exact numerical lattice calculations are performed to test and complement our findings for the free model.Section 4 contains instead the analysis of the model at finite interactions.After a short summary of the Bethe Ansatz solution in Sec.4.1, in Sec.4.2 we consider the hydrodynamic limit of the spin chain and we characterize the initial state at t = 0 using the local density approximation.The GHD equations are reported in Sec.4.3, together with a detailed derivation of their analytic solution.The exact computation of the profiles follows in Sec.4.4, where one can also find numerical checks based on time-dependent Density Matrix Renormalization Group (tDMRG).The analytic treatment of the entanglement evolution in the interacting case is very challenging as it needs a careful analysis based on quantum GHD [43].This study goes well beyond the goal of this paper.Nevertheless, in Sec.4.5, we present a numerical analysis of the entanglement and we compare our numerical findings with the exact solution of the non-interacting case.Finally, we report our conclusions in Sec. 5 and some technical aspects of the Bethe Ansatz solution in Appendix A.

Setup of the problem
We consider the one-dimensional XXZ model with Hamiltonian where L is the length of the chain and σα=x,y,z x are standard Pauli operators acting on site x.Here, ∆ is the interaction coupling which we set to be in the gapless regime |∆| ≤ 1, where it is customary to write ∆ = cos γ.We do not consider |∆| > 1 because energy arguments show that domain wall states do not melt, see Ref. [46,47].The case |∆| = 1 is very peculiar [48,49] and therefore it will be also excluded from this study.We further focus on the rational case, i.e., on those values of γ that can be written as the ratio γ = πQ/P, with Q, P two co-prime integers, 1 ≤ Q < P. In the rational case, γ admits a continued fraction representation where {ν 1 , . . ., ν q } is a set of numbers satisfying ν 1 , . . ., ν q−1 ≥ 1 and ν q ≥ 2. In the thermodynamic limit L → ∞, the model ( 1) is exactly solved by means of Thermodynamic Bethe Ansatz (TBA).In particular, for large values of the system size L and under the string hypothesis [55], one can describe the excitation spectrum of the spin chain in terms of different species of quasiparticles (generically referred to as strings).A short summary of the TBA solution of the model (1) will be reported in Sec.4.1 while we address the reader to e.g.Ref. [55,56] for a comprehensive discussion.The total number of allowed strings δ is read from the interaction coupling in Eq. ( 2) as Each of the strings has a well-defined length l j , parity u j and sign s j , which are all given in terms of {ν 1 , . . ., ν q } after some manipulations [40,55], see Appendix A. At t = 0, we prepare the system in the product state where ±ℓ are the positions of the domain walls and |↑ x 〉 (resp.|↓ x 〉) is the eigenstate of σz x with eigenvalue +1 (resp.−1).For t > 0 we consider the Hamiltonian dynamics generated by ( 1) during which the central ordered domain of (4) gradually melts into the left and right ferromagnets, see Fig. 1 for an illustration.Our focus will be on the hydrodynamic limit of the model, where several exact results can be derived in the limit of large space and time scales.Specifically, we shall consider the scaling limit x, ℓ, L, t → ∞ at fixed ratios x/t and ℓ/L and we shall investigate the non-equilibrium dynamics of the local magnetization m and of the local spin current J, defined as The initial state (4) has zero entanglement for any bipartition.However, the latter is generated during the melting dynamics through the spreading of quasiparticles from the junctions towards the ordered regions.Therefore, we will investigate the formation and the subsequent growth of the entanglement entropy by focusing on the Rényi entropy of a subsystem A = [x, +∞] and its limit α → 1, i.e., the Von Neumann entanglement entropy In Eqs. ( 7), ( 8), ρA = tr Ā ρ(t) is the reduced density matrix of the subsystem A, obtained by tracing out the degrees of freedom of the interval Ā = [−∞, x) from the full density matrix ρ = |Ψ(t)〉 〈Ψ(t)|.

Analytic solution of the non-interacting case
We first investigate the melting dynamics of the double DW (4) in a non-interacting spin chain, leaving the discussion of the interacting model to Sec. 4. After setting ∆ = 0 in Eq. ( 1), the spin chain Hamiltonian reduces, up to an irrelevant additive constant, to the free Fermi gas where ĉ † x , ĉx are standard lattice Fermi operators, obtained from the ladder Pauli operators σ± x = ( σx x ± i σy x )/2 after a Jordan-Wigner transformation [57].Although the non-equilibrium evolution of the free Fermi gas (9) can be investigated by means of exact lattice calculations [4,5], we shall focus on a hydrodynamic limit.Indeed, a large-scale description gives not only access to asymptotically exact results for several quantities of interest (such as the magnetization profile and the spin current) without difficult computations, but also, it allows us to characterize the entanglement evolution during the melting process, whose lattice derivation is currently out-of-reach with standard techniques.

Semi-classical hydrodynamics
In the thermodynamic limit L, ℓ → ∞ at fixed ratio ℓ/L, we describe the local physics of the model (9) in terms of coarse-grained cells, each containing a large number of lattice sites.Inside each cell, we diagonalize the Hamiltonian (9) in Fourier space, obtaining [6] where the lattice index x is replaced by a continuous variable.Here, η † k,x , ηk,x are the creation and annihilation operators of a fermionic particle with momentum k inside the cell labeled by x.In this limit, it is easy to see that the initial state in Eq. ( 4) corresponds to a gas where fermionic particles with momentum −π ≤ k ≤ π entirely fill the spatial region −ℓ ≤ x ≤ ℓ, leaving empty the rest of the chain.Therefore, in terms of the Wigner function n(x, k) [58,59] (which is essentially the occupation function of the free Fermi gas), the macrostate at t = 0 corresponding to the double DW state (4) is given by n At times t > 0, the fermionic particles propagate independently with constant velocity v(k) = sin k.This implies that the occupation number at time t > 0 can be obtained by tracking backward the trajectory of each particle up to t = 0 and it reads Figure 2: Snapshots of the Wigner function ( 12) during the melting dynamics of the double DW (11).
The result ( 12) might be also viewed as solution of a Euler-like hydrodynamic equation [6, 54] satisfied by the Wigner function, at lowest order in the ∂ x and ∂ k derivatives [60][61][62][63].A microscopic derivation of Eq. ( 13) for non-interacting fermions can be found in Refs.[60,61].Similar studies in higher spatial dimensions d > 1 and at non-zero temperature can be found in Ref. [64,65].In Fig. 2, we show the hydrodynamic evolution of n(x, t, k) during the double DW melting process.
Notice that the Wigner function in Eq. ( 12) shows a particularly simple form as it takes only the values 0 or 1, due to the zero-entropy condition of the local states during the time evolution [66].As a consequence, it is possible to encode all the information about the phasespace dynamics in terms of the contour of the Wigner function, typically called Fermi contour Γ (t), which keeps track of the Fermi points at each position x.
From the inspection of Fig. 2, one can graphically determine the Fermi points at each x and t.In particular, the number n of Fermi points is given by the number of intersections of a vertical line drawn at position x and time t with the contour of n(x, t, k).With this method, we observe that two Fermi points are found at each position x as long as t < ℓ while, for t ≥ ℓ, one enters in a richer landscape where split Fermi seas can be found.This behavior of the Fermi seas can be easily understood considering lightcone regions ||x| − ℓ|/t ≤ 1 around the two junctions x = ±ℓ, determined by the propagation of the fastest modes with k = ±π/2 and velocity v(±π/2) = ±1, see Fig. 3(a) for an illustration.At times t < ℓ, the dynamics of the particles inside the two lightcones is independent and thus we find a connected Fermi sea at any position x.In this regime, the evolution is that of two independent domain walls [4].Conversely, at times t ≥ ℓ, the two lightcones intersect over the spatial region x ∈ [−t + ℓ, t − ℓ] and, as a consequence, a split Fermi sea is found.Precisely, the number n of Fermi points in each regime is Since the propagating modes on the Fermi contour are only those initially located at x = ±ℓ, the analytic expression of the Fermi points follows from the equation of motion Figure 3: (a) Illustration of the particles spreading during the melting dynamics.From the junctions at x = ±ℓ, two lightcones open up and determine the melted region ||x| − ℓ| < t.At early times t < ℓ these lightcones do not intersect and a connected Fermi sea is found at any position x (light-orange region).Conversely, when t ≥ ℓ, the particles emitted from one junction penetrate inside the lightcone of the other, leading to a split Fermi sea for |x| ≤ t −ℓ (dark-orange region).(b) Graphical solution of Eq. ( 15) at fixed time t: In each regime, the values of the roots k (s) F (x, t) are given in Eq. ( 17).Thick red lines denotes the intial position ±ℓ of the two domain walls.which is readily solved as By exploiting the symmetry properties of the contour and with elementary algebra, one can show that the four roots in Eq. ( 16) combine together to give the following Fermi points in each regime: with sign determined by sign(k as illustrated in Fig. 3(b).Denoting such Fermi points with k F in increasing order), and the (split-)Fermi sea Γ (x, t) with one can derive exact asymptotic results by summing up each individual contribution coming from a filled mode at position x and time t, i.e., for each k ∈ Γ (x, t).For instance, the asymptotic result for the magnetization profile is given by Explicitly, using Eq. ( 17) in (20), we obtain Similarly, the spin current in Eq. ( 6) is obtained as the weighted sum In Fig. 4, we show the exact asymptotic formulae for the profiles in Eq. ( 21)-( 22) against exact lattice numerical calculations finding a perfect agreement.The numerical data are obtained from the lattice evolution of the two-point function In particular, we determine G x,x ′ (0) from the exact diagonalization of the Hamiltonian with auxiliary potential which is chosen such that its ground state reproduces the initial state of Eq. ( 4).Subsequently, we evolve the two-point function in the eigenstate basis of the Hamiltonian (9), see Ref. [6,20] for a detailed discussion on the numerical implementation.From the exact evolution of the two-point correlation (24), it follows that

Quantum fluctuations and Entanglement entropy
It is important to notice that the phase-space approach discussed in Sec.3.1 does not include quantum fluctuations of the Fermi contour and therefore, it is not sufficient to characterize the entanglement spreading during the melting dynamics.On the other hand, an exact lattice calculation for such non-homogeneous and non-equilibrium processes is currently not possible with standard techniques, even for the free model (9) considered in this section.To overcome this limit, it has been recently noticed [6,41,43,54] that the asymptotic behavior of the entanglement can effectively described by a quantum hydrodynamic theory, obtained after m(x, t)  21)-( 22) at different times while the symbols are obtained with exact numerical calculations on a lattice of size L = 600 and ℓ = 100.The matching of the profiles with the numerics is extremely good.
the re-quantization of the Fermi contour of Sec.3.1 in terms of a Luttinger-liquid [67].The latter has to reproduce the relevant low-energy quantum physics in the vicinity of the Fermi points, which essentially consists in the formation of particle-hole excitations around the Fermi contour.Following this program, we consider the quantum fluctuations of the fermionic density in terms of a fluctuating-field χ δ ρ(x, t) = and we expand, at leading order in the scaling dimension of the low-energy fields, the timedependent creation and annihilation fermionic operator with standard bosonization method up to a non-universal amplitude and a semi-classical phase which are not important for our scopes (see e.g.Ref. [6,7,54] for more details on the phase and e.g.Ref. [70,71] for the amplitude).In Eq. ( 29), we have introduced the chiral components of the fluctuating field χ = χ+ + χ− which carry the physical interpretations of right (+) and left (−) moving excitations along the Fermi contour.It is therefore useful to introduce a parametrization of the Fermi contour Γ (t) in terms of a coordinate θ along the curve as The dynamics of the chiral fields χ± is then governed by the conformal field theory of a free massless boson or, equivalently, by the Luttinger-liquid Hamiltonian along the Fermi contour (see e.g.[6,41,43,54,71,72]) Figure 5: Parametrization of the initial Fermi contour with the coordinate θ along the unit circle.Notice that the horizontal branches of the contour do not contain propagating modes and therefore are not needed for the description of the entanglement evolution.
where J (θ ) is a jacobian factor and a = ± if k(θ ) ≷ 0. Importantly, in our hydrodynamic description of the problem, quantum fluctuations in the initial state are given by the ground state of ĤLL [Γ (0)] and the effect of the time-evolution is simply to transport such quantum fluctuations along the curve Γ (t), which gets deformed over time according to the semi-classical dynamics of Sec.3.1.
With this framework, we now consider the large-scale contribution to the Rényi entropy in Eq. ( 7) of a subsystem A = [x, +∞] by using a conformal field theory (CFT) approach.For integer Rényi index α, the latter can be expressed as the expectation value of a twist field Tα living at the boundaries of the subsystem A [73-75] Moreover, in our chiral model, Tα admits a decomposition in chiral twist fields Φ− α , Φ+ α that behave under conformal mappings as primary fields with scaling dimension It follows that Sα (x, t) can be written as the n-point correlation function of the chiral fields Φ− α , Φ+ α of the CFT which lives along the Fermi contour at time t.Because of the conservation of momentum during the time-evolution, the Fermi points k (s) F (x, t) at time t can be traced backward to the initial Fermi contour where they can be simply parametrized as as illustrated in Fig. 5. Denoting with θ s the positions of k (s) F (x, t) along the initial Fermi contour, the Rényi entropy reduces to and it can be exactly determined, as we now discuss.From Eq. ( 16) we obtain the Weyl factors while the n-point correlation function is given by where Plugging Eq. ( 36) and (37) into Eq.( 35) and with simple algebra, we arrive to the result where we introduced the shorthands The CFT result in Eq. ( 39) provides only the universal contribution to the Rényi entropies and it has to be complemented with a non-universal cutoff ε(x, t) so that For a non-interacting spin-1/2 chain, the exact expression of the cutoff for the von-Neumann entanglement entropy is known [73,76].In particular, for n = 2, i.e., for a connected Fermi sea, it is given in terms of the magnetization profile in Eq. ( 20) as [6, 7] where C is a known non-universal amplitude, see Ref. [76] and Eq. ( 46) below.In the presence of the split Fermi sea, the computation of ε is more involved but it can be still performed analytically exploiting the Fisher-Hartwig conjecture.As a result, one obtains the generic formula [77] One can check that Eq. ( 43) correctly reproduces the cutoff in Eq. ( 42) for n = 2 while, in the regime with four Fermi points, it gives Hence, including the cutoff contribution ( 42)- (44) in Eq. (39) and taking the replica limit α → 1, we finally obtain the result for the entanglement entropy Figure 6: Evolution of the entanglement profiles in the double domain wall melting of a non-interacting spin-1/2 chain.The symbols are obtained with exact lattice numerical calculations while the full-line is our result in Eq. ( 45).The numerical data are obtained for a chain of size L = 600 and setting ℓ = 100.
where Υ is a non-universal constant [76], related to C as In the regime |t − ℓ| < |x| ≤ t + ℓ where the two lightcones do not intersect, we observe that Eq. ( 45) is in agreement with the exact prediction for the entanglement spreading of a DW state, obtained in Ref. [7].On the contrary, for |x| ≤ t − ℓ, we find a non-trivial behavior of the entanglement that cannot be reduced to just a sum of the two individual contributions coming from each lightcone.We tested our formula (45) for the entanglement entropy against exact lattice numerical calculations finding an excellent agreement, see Fig. 6.The numerical data are obtained from the two-point correlation matrix in Eq. ( 24) restricted to an interval A of length |A| as where ξ j (t), j = 1, . . ., |A| are the eigenvalues of G A (t), see e.g.Ref. [6,20] for more details.
It is interesting to consider the dynamics of the half-system entanglement by setting x = 0 in (45).In this case, the variables in Eq. ( 40) simplifies to ϕ = −φ = arcsin(ℓ/t) and the half-system entanglement entropy reads Moreover, by expanding Eq. ( 45) for t ≫ |x|, ℓ at leading order, it is easy to show that   49) (full-line) is compared with the numerical data (symbols) obtained for a lattice of size L = 350 with ℓ = 25.At times t ≫ ℓ, the half-system entanglement saturates to the constant value 1 3 log ℓ + 2Υ (dashed horizontal line).(b) Asymptotic saturation of the entanglement profiles to the value 1  3 log ℓ + 2Υ for spatial positions |x| ≪ t, see Eq. ( 50).
The saturation of entanglement in Eq. ( 49)-( 50) can understood with a simple CFT argument that we now discuss.At large but fixed time t ≫ |x|, ℓ, the bulk of the system is characterized by a homogeneous magnetization profile (see Eq. ( 21)) and by a correlated region of size (see Fig. 3(a)) This asymptotic setup can be equivalently interpreted as a homogeneous finite-size system of length L(t) containing two species of fermionic particles (one from each junction) with density ϱ = (m + 1/2)/2 ∼ ℓ/(πt), whose half-system entanglement is known [73,84] and reproduces the result in Eq. ( 50) after setting c = 2 and the short distance cutoff to ε(t) = C/ sin(πϱ) ∼ C t/ℓ (see Eq. ( 42) and Ref. [76]).In Fig. 7, we compare the results in Eq. ( 49)-( 50) with exact lattice calculations of the half-system entanglement.As one can see, our exact formula is in perfect agreement with the numerical data.Notice that the result in Eq. ( 50) agrees with that of Ref. [85], obtained with exact lattice calculations.This gives a non-trivial check about the validity of our hydrodynamic approach.

Analytic solution of the interacting case
We now turn to the discussion of a gapless XXZ model (1) for general |∆| < 1.As we argue below, the presence of non-zero interactions requires some adjustments in our hydrodynamic description but it does not qualitatively affect the physical dynamics obtained from the solution of the non-interacting limit in Sec. 3.

Bethe Ansatz solution
As anticipated in Sec. 2, the Hamiltonian ( 1) is diagonalized by means of the Bethe Ansatz [55,56].In particular, the eigenstates of the model can be written down exactly and are labeled by a set of complex numbers {λ j } (also called rapidities) which generalize the concept of particle momenta of a free Fermi gas to the interacting case.The value of {λ j } is not arbitrary but it is found as solution of a set of non-linear algebraic equations called Bethe Ansatz equations.In the Hilbert space sector with M spins up and assuming periodic boundary conditions of the chain, the Bethe Ansatz equations take the form [55] sinh(λ j + iγ/2) Eq. ( 54) implements non-trivial quantization conditions for the interacting model ( 1) on a ring of length L. In the thermodynamic limit L → ∞ and according to the string hypothesis [55], the solutions of Eq. ( 54) are arranged into strings i.e., in regular patterns in the complex plane composed by l j rapidities having the same real part λ β j and equidistant imaginary parts Here, δ is the total number of strings (given in Eq. ( 3)), l j is the length and u j is the parity of a given string while β is a label for the strings of a given species, see Appendix A for more details and e.g.Ref. [55,56] for a more systematic treatment.Deviations from Eq. ( 55) are expected to vanish exponentially with the system size when L → ∞.The different strings are interpreted as δ different species of quasiparticles, each described by a real rapidity λ β j , corresponding to the string center.Since for L → ∞ the spectrum becomes densely populated, it is useful to introduce a macroscopic spectral distribution of quasiparticles that satisfy the following integral equation [55] where and we recall that s j denotes the sign of the string (see Appendix A).The solution of Eq. ( 57) allows us to characterize the ground state of a translationally-invariant Hamiltonian (1) in the thermodynamic limit.For more generic states, the l.h.s. of Eq. ( 57) is modified as ρ j → ρ j + ρ h j , with ρ h j a distribution function for the unoccupied quasimomenta.The relation between ρ j and ρ h j is not determined by Bethe Ansatz alone but it requires additional thermodynamic arguments [55].For later convenience, we introduce also the occupation function n j (λ) = ρ j (λ) which will play the same role of the Wigner function of Sec.3.1 for the interacting model ( 1) in its hydrodynamic description, see Sec. 4.3.

Large-scale description of the initial state
Let us now consider a large-scale asymptotic description of the model (1) in terms of periodic boxes ∆x, each containing a large number of lattice sites.In this way, we can determine the local macrostate n j (x, λ) by solving the TBA equations of Sec.4.1 within each coarse-grained cell x.At this point, we are in need of a large-scale description of the macrostate corresponding to the intial state in Eq. ( 4).To achieve this goal, we first investigate the TBA solution of a generic spin-state where h is an external magnetic field and Ŝz = x is the total spin in the z-direction.In this simple instance, the occupation functions are known exactly [40,55].In particular, the first δ − 2 strings are n (h) where { y −1 , y 0 , . . ., y q } and {m 0 , m 1 , . . ., m q } are two set of numbers related to {ν 1 , . . ., ν q }, as reported in Appendix A. The second-to-last and the last strings have instead a different behavior n In the limit of strong magnetic field |h| → ∞, the spin-state (60) becomes fully-polarized and it reduces to the projector where |⇑〉 ≡ x |↑ x 〉 and similarly for |⇓〉.As a consequence, the occupation functions in Eq. ( 61)-( 62) simplify to signaling that only the largest two strings are responsible for the thermodynamics of the polarized state (63) [40,41].This observation is at the basis of our analytical solution of the double domain wall melting problem, as discussed in the following section.
From the solution (64) of a polarized spin state, we can easily derive the occupation functions n j (x, λ) corresponding to the double DW state (4).In particular, we associate to each coarse-grained cell x the polarized state so that the initial occupation functions for our setup read with Θ the Heaviside step-function.

Generalized Hydrodynamics
The hydrodynamic evolution of the initial state ( 66) is established by means of Generalized Hydrodynamics (GHD),which provides the following set of continuity equations for the occupation functions [35,42] This set of equations has the same structure of Eq. ( 12), obtained in Sec.3.1 for the noninteracting model, but it is characterized by an effective velocity v eff j = v eff j (x, t, λ) which depends self-consistently on the state n j (x, t, λ) through the dressing operation [35,40] where q j is a generic function 1 .For the specific case of Eq. ( 68), the energy e(λ) and the derivative of the momentum eigenvalue p ′ (λ) are e(λ) = −π sin(γ)a j (λ) , p ′ (λ) = 2πa j (λ) (70) and their dressed value is obtained as solution of Eq. ( 69).
The GHD equations ( 67) are exact only in the limit of large space and time scales.Nevertheless, they proved to give a very accurate description of the dynamics even for relatively small times and distances, as confirmed by several numerical tests e.g.[35,37,40,42,[89][90][91][92][93] and even some experimental studies [52,94,95].Therefore, we will make use of the GHD equations as a basis for our analysis, providing further numerical tests based on time-dependent Density Matrix Renormalization Group (tDMRG) simulations of our results.
A formal solution of GHD equations ( 67) is obtained with method of characteristics, yielding where Plugging the initial state (66) in Eq. ( 71), it is easy to see that only on the largest two strings are relevant for the dynamics since Moreover, the structure of the interaction kernels in Eq. ( 58) reveals the following symmetry properties of the last and the second-to-last strings [40] and the signs s δ = −s δ−1 .This implies that the dressing operation (69) of the largest two strings gives and, with the same calculation, one can show that e dr′ δ = −e dr′ δ−1 .It follows that the largest two strings share the same effective velocity (68) v eff δ = v eff δ−1 , which in turns implies from Eq. ( 73) that the occupation functions From this observation, using Eq. ( 77) in Eq. ( 75) (and similarly for the dressed energy), we conclude that the dressing operation becomes trivial and the effective velocity v eff δ reduces to where we defined ζ 0 = sin γ/ sin(π/P) .

Magnetization and Spin-current profiles
The analytic solution (80) of GHD allows us to derive exact asymptotic results for the conserved charges and currents of the model.Proceeding similarly to what done in Sec.3.1, we define the (split-) Fermi sea of the interacting model as This equation of motion can be easily solved by noticing that the quantity s δ p δ (λ) is a monotonic function in the interval [−π/P, π/P].Therefore, in analogy with the free case, we define the function k(λ) ≡ s δ p δ (λ) and we obtain the roots that arrange together to give the following Fermi points with sign(k (s) F (x, t)) = sign(x(|x| − ℓ)).We use then the Fermi points (83) to construct the (split-)Fermi sea Γ (x, t) as in Eq. (19).Notice that the structure of Γ (x, t) is very similar to that found in the non-interacting case (cf.Eq. ( 17)) although the presence of the interactions is responsible for a rescaling of the lightcones with ζ 0 .At this point, we are ready to compute the non-equilibrium dynamics of the charges profiles.In particular, the profile of the magnetization is obtained as the weighted sum over the single quasiparticle contributions that, in the case of a double DW melting dynamics, reduces to Figure 10: (a) Numerical analysis of the half-system entanglement at large times t ≫ ℓ for a lattice of size L = 120 with ∆ = 0.5 and different domain sizes ℓ.We observe a saturation of the entanglement to a constant value that we extracted with a best fit of our data (dashed-lines).(b) Comparison of the numerical results with the Ansatz in Eq. ( 88) using the proxy (89) for ∆ = 0.5, 0.7.Although finite-size effects affect our data, we see a good agreement, which we expect to improve with larger values of L, ℓ and t.In the plot, error bars are extracted from the covariance matrix of the fitting plateaus.
it has been shown [41] that the low-energy field theory is characterized by a constant Luttinger parameter, which displays a fractal dependence on ∆.In this simple instance, conformal invariance is not broken and an asymptotic prediction for the asymptotic growth of the halfsystem entanglement can be derived with CFT scaling arguments [41].Though, in our double DW setting, this is not the case and we find the exact computation of the entanglement highly non-trivial.Nevertheless, from the numerical study in Fig. 9, we observe that the behavior of the entanglement profiles is not qualitatively influenced by the presence of interactions.In particular, in the regime |ζ 0 t − ℓ| < |x| ≤ ζ 0 t + ℓ, we observe that the entanglement spreading is given by two independent contributions that are developed around each junction.Conversely, when |x| ≤ tζ 0 − ℓ, we find an interplay of the quasiparticles emitted from the two walls, which results into a fast growth of the half-system entanglement entropy, similarly to what observed in Sec.3.2 for a free model.Motivated by this strict analogy with the free case, we expect an asymptotic saturation of entanglement and we conjecture the following form where c(γ) is a non-universal constant depending on the parameter γ that we are unable to determine.However, by considering the entanglement evolution of two initial domains of upwards spins with different sizes (respectively ℓ and ℓ ′ ), we can test our conjecture (88) considering the difference We show the results of our simulations in Fig. 10.In particular in Fig. 10(a), we observe an asymptotic saturation of the half-system entanglement for different values of ℓ.The panel 10(b) is instead a numerical test of our Ansatz (88), where we construct our proxy in Eq. ( 89) with the plateaus obtained as in Fig. 10(a).We find a good agreement for different values of ∆ already for modest values of L and ℓ.

Summary and conclusions
We investigated the non-equilibrium dynamics of a spin-1/2 XXZ model (1) at zero temperature, initially prepared in a state with two domain walls (4).This quench setup is a generalization of the widely studied bipartitioning protocols (see e.g.[4-8, 35, 37, 40-42]), where now a central domain of upwards spins x ∈ [−ℓ, ℓ] is joined to two semi-infinite down-oriented ferromagnets, respectively on its left [−L/2, −ℓ − 1] and right [ℓ + 1, L/2] sides.In particular, we focused on the rational case, where ∆ = cos γ and γ = πQ/P, with Q, P two co-prime integers.In this case, as first noticed in Ref. [40], the structure of the Bethe Ansatz solution allows us to analytically solve the GHD equations ( 67) and, as a consequence, to determine the exact asymptotic evolution of the charges and currents profiles during the melting dynamics (see Eq. ( 85) and ( 87)).The non-trivial result of this paper is that such exact solvability of the model extends also to a late-time regime characterized by the presence of correlations between the two domain walls, and therefore, by the emergence of split Fermi-seas.To our best knowledge, this is the first case where a fully-analytical GHD solution is found with multiple Fermi points.In Sec. 3, we also treated the case of a non-interacting spin chain (i.e. a free Fermi gas initially confined in a segment of size 2ℓ) in detail.For the free gas, after the derivation of the charges profile with standard hydrodynamic techniques (see Sec. 3.1), we applied recent developments in quantum fluctuating hydrodynamics [6,54] to give an exact description of the entanglement entropy (see Eq. ( 45)).Moreover, we found an asymptotic saturation of the half-system entanglement entropy as S 1 (0, t) ∼ 1 3 log ℓ.Motivated by the strict analogy of the melting dynamics in the free and interacting case, we expect a similar behavior for the interacting chain (88), as confirmed by several numerical simulations (see Fig. 10).
We mention that the exact solution of the GHD equations reported in this work can be, in principle, extended to a whole class of spin chain polarized states such as multiple domain wall configurations and spin-helix states [97].Nonetheless, it could be interesting to investigate the entanglement evolution of a double domain wall state with quantum GHD [41,43] and, in this way, to prove our conjecture (88) about the half-system entanglement growth.In terms of these series of numbers we have the following relations for the length l j l j = y i−1 + ( j − m i ) y i for m i < j < m i+1 , and l δ ≡ y q ; (91) the parity u j u j = (−1) (l j −1)Q/P for m i < j < m i+1 , and u δ ≡ (−1) q ; (92) and the sign s j s j = sign(ω j ) , for m i ≤ j < m i+1 , where ω j ≡ (−1) i (p i − ( j − m i )p i+1 ) and the series {p 0 , . . ., p q } is defined by As a concrete example, let us consider the case with Q ≡ 1 (hence, γ = π/P), known also as root of unity point.In this instance, one finds δ ≡ P strings (see Eq. ( 3)) with the following properties l j = j , u j = 1 , ω j = P − l j j ∈ {1, . . ., P − 1}; (95)

Figure 1 :
Figure 1: Illustration of the melting dynamics of the double domain wall state in Eq. (4).At t = 0 the system is prepared in the product state of three ordered domains, resulting into a step-like shape of the magnetization profile m(x, t) (left panel).(a) At early times, the central domain of aligned spin begins to melt into the left and right ferromagnets generating non-homogeneous propagation fronts around the junctions.After the initial transient regime, the two propagating fronts merge (b) and contribute together to the relaxation (c)-(d).

Figure 4 :
Figure 4: (a) Magnetization and (b) spin current profiles for the non-interacting double domain wall melting problem.The different curves show the behavior of the analytic solutions (21)-(22) at different times while the symbols are obtained with exact numerical calculations on a lattice of size L = 600 and ℓ = 100.The matching of the profiles with the numerics is extremely good.

Figure 7 :
Figure7: (a) Half-system entanglement entropy for a non-interacting spin-1/2 chain initially prepared in a double domain wall state (4) as function of time.The exact prediction in Eq. (49) (full-line) is compared with the numerical data (symbols) obtained for a lattice of size L = 350 with ℓ = 25.At times t ≫ ℓ, the half-system entanglement saturates to the constant value1  3 log ℓ + 2Υ (dashed horizontal line).(b) Asymptotic saturation of the entanglement profiles to the value1  3 log ℓ + 2Υ for spatial positions |x| ≪ t, see Eq. (50).

Figure 8 :
Figure 8: (Top) Magnetization and (bottom) spin current profiles during the double domain wall melting dynamics for several values of the interaction coupling γ (corresponding to ∆ = 0.223, 0.309, 0.5, 0.707, from the left to the right panels).The different curves show the behavior of the analytic solutions (21)-(22) at different times while the symbols are obtained with tDMRG numerical simulations obtained for a lattice of size L = 225 and ℓ = 35.The matching of the profiles with the numerics is seen to be very good although finite-size effects are still visible.