Scaling of fronts and entanglement spreading during a domain wall melting

We revisit the out-of-equilibrium physics arising during the unitary evolution of a one-dimensional XXZ spin chain initially prepared in a domain wall state |ψ0⟩=|⋯↑↑↓↓⋯⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vert \psi _0\rangle =\vert \dots \uparrow \uparrow \downarrow \downarrow \dots \rangle$$\end{document}. In absence of interactions, we review the exact lattice calculation of several conserved quantities, including e.g. the magnetization and the spin current profiles. At large distances x and times t, we show how these quantities allow for a ballistic scaling behavior in terms of the scaling variable ζ=x/t\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta = x/t$$\end{document}, with exactly computable scaling functions. In such a limit of large space-time scales, we show that the asymptotic behavior of the system is suitably captured by the local occupation function of spinless fermionic modes, whose semi-classical evolution in phase space is given by a Euler hydrodynamic equation. Similarly, analytical results for the asymptotic fronts dynamics are obtained for the interacting chain via Generalized Hydrodynamics. In the last part of the work, we include large-scale quantum fluctuations on top of the semi-classical hydrodynamic background in the form of a conformal field theory that lives along the evolving Fermi contour. With this procedure, dubbed quantum generalized hydrodynamics, it is possible to obtain exact asymptotic results for the entanglement spreading during the melting dynamics.


Introduction
In the last two decades or so, the study of transport phenomena in one-dimensional quantum systems is experiencing renewed interest due to the recently developed theoretical and experimental technologies.Special attention has been dedicated to one-dimensional integrable systems, where unusual transport properties emerge due to the presence of an infinite set of conserved charges.In particular, such conservation laws are responsible for a lack of thermalization of the system at large times and, consequently, for the emergence of current-carrying steady states, denoted as generalized Gibbs ensembles (GGE) [1].
The scope of this short review is to provide a concise (and yet comprehensive) discussion of the physics of the DW melting by revisiting early results and modern GHD approaches to the problem, thus giving a handbook of the main available results and of their derivation.The paper is organized as follows.In Sec. 2, we introduce the model, the quench protocol, and the quantities that we wish to characterize during the dynamics.The contents are organized in two main sections: i) Sec. 3 dedicated to the study of the profiles of conserved charges and currents by means of exact lattice calculations and hydrodynamics; ii) Sec. 4, where we employ a requantization of the hydrodynamic theory for the study of entanglement.In both sections, the study of the non-interacting case is treated separately for a better exposition.In Sec. 5, we give a short summary and some concluding remarks.

The domain wall melting problem
We consider the one-dimensional XXZ model with Hamiltonian where N is the number of lattice sites (which we assume to be even for future convenience), σa j , a = x, y, z, denotes the standard Pauli operators acting on site j and ∆ is the anisotropy parameter.The hopping amplitude J defines the overall energy scale and will be set to unity from hereafter.We shall focus on the gapless regime of the system |∆| < 1, with the usual parametrization ∆ = cos γ.We do not consider the case |∆| > 1 because it can be shown with energy arguments that initial domain wall configurations do not melt, see Refs.[87,88].The case |∆| = 1 is very peculiar [89,90] and therefore it will be also excluded in this short review.
At t = 0, we prepare the system in the product state where | ↑ j (resp.| ↓ j ) is the eigenstate of σz j with eigenvalue +1 (resp.−1), see Fig. 1 for an illustration.For t > 0 we consider the Hamiltonian dynamics generated by (1) during which the initial DW state (2) gradually melts.During this process, we wish to characterize the transport properties of the system.To this end, we will investigate the non-equilibrium dynamics of the profiles of conserved charges and currents of the integrable model.For the sake of concreteness, we will highlight the main features related to the melting dynamics by considering the magnetization profile and the spin current Figure 1 Setup of the domain wall melting problem.At t = 0 the spin chain is prepared in the product state (2), corresponding to a step-like profile of magnetization between the values ±1/2.For t > 0, this configuration is let to unitarily evolve with Hamiltonian (1).
The latter can be deduced from the fact that the total magnetization is conserved by the unitary dynamics and, therefore, a local magnetization current can be defined using a lattice continuity equation.This extends to any conserved quantity and associated current in the model, which are all obtainable with simple hydrodynamic arguments as reviewed in the next section.
In the last part of this short review, we shall focus on the entanglement spreading.In particular, we will characterize the entanglement of a bipartition of the onedimensional chain as for which the reduced density matrix is obtained as and, from it, the Rényi entropies are defined as with index p ∈ N. The entanglement entropy is then obtained in the limit p → 1 and reads as

Transport properties during the DW melting
As discussed below, the physics of the DW melting is qualitatively similar in the whole region |∆| < 1.However, we find convenient to discuss the non-interacting case ∆ = 0 separately since, in this case, the system reduces to a set of free fermions and the transport properties can be determined with exact lattice calculations.From this, an emergent hydrodynamic picture can be rigorously determined.In the interacting case instead, the characterization of the transport during the DW melting is possible only by means of GHD and it requires the review of the Bethe ansatz solution of the interacting model.

Non-interacting case
For ∆ = 0, the Hamiltonian (1) reduces to the XX model.The latter can be mapped to a system of spinless non-interacting fermions with tight-binding Hamiltonian up to a boundary term which is negligible when N → ∞, as considered below.Here ĉ † j , ĉj stand for the creation and annihilation operators of a lattice spinless fermion at site j and they satisfy canonical anticommutation relations {ĉ † i , ĉj } = δ i,j .The two forms of the Hamiltonian (9) are related through Jordan-Wigner transformation [91,92] where σ± j are the ladder combination of Pauli operators σ± j = (σ x j ± iσ y j )/2.In Fourier space, the model ( 9) is readily diagonalized in terms of the operators with single-particle wavefunction χ q,j = √ 2N −1 sin(k q j) and quantized momenta k q = πq/N , q = 1, . . ., N .
Taking the limit N → ∞, it is possible to replace the quantized momentum k q with a continuous variable k ∈ [−π, π] and use the single-particle wavefunctions χ k,j = exp(ikj)/ √ N .

Lattice description of the dynamics
From the exact solution of the non-interacting spin chain (11), one can determine the time evolution of the operators b † k in Heisenberg picture as and obtain, via inverse Fourier transform, the time evolution of the lattice fermionic operators as where J ν (z) are Bessel functions of the first kind.From this observation, the expectation value of lattice operators can be obtained by elementary calculations.For instance, the magnetization profile is given by The initial domain wall state (2) forces j ≡ l in the previous double sum leading to [47,53,57,61,93] where the last expression is obtained using the properties In Fig. 2, we compare the exact result in Eq. ( 15) with numerical calculations for the finite lattice model, finding an excellent agreement.
In general, given the non-interacting nature of the problem, one can determine exact results for the transport properties during the DW melting from the knowledge of the two-point function which can be simplified, thanks to the summation theorem [94] for Bessel functions, to [47,53,57,61,93] (18) For instance, the magnetization profile in Eq. ( 15) is recovered from the diagonal of G n,m as and the spin current can be obtained as Similar results for other quantities can be worked out following the same lines.

Scaling limit of fronts
In the limit n → ∞, t → ∞ at fixed ratio ζ = n/t, the magnetization profile admits a scaling form with a scaling function Φ(ζ) that can be analytically determined from Eq. ( 15) using the asymptotic expansion of J n (n/ζ) when n → ∞.From this calculation, one obtains Eq. ( 22) highlights the presence of a non-equilibrium steady state (NESS) for this quench protocol, characterized by a non-homogeneous profile of magnetization.Since σz j is a conserved quantity for the spin chain (1), it follows that the NESS is characterized by the presence of a non-vanishing spin current j M n (t) that can be determined from the continuity equation where ∇ n is a discrete derivative defined on the lattice.Using Eq. ( 22) in (23), one finds the scaling form for the NESS spin current The same result is obtained from Eq. ( 20) by taking the asymptotic behavior of the Bessel functions.

Emergent hydrodynamic description
In the limit N → ∞ and for a translationally-invariant Fermi gas, the occupation function of fermionic modes can be obtained from the Fourier transform of the twopoint correlation function in real space where the site j is arbitrary because of translational invariance.Given the noninteracting nature of the problem, n k is expected to fully characterize the properties of the fermionic gas.In order to investigate the DW melting problem then, one can consider j = ζt and take the limit t → ∞ before performing the summation over the index n [47].In this way, the local properties of the fermionic gas are characterized by the steady-state value of a region which is shifted from the orgin of an amount ζt.Using the expression ( 18) for the two-point function, and taking the asymptotic expansion of the Bessel functions, one arrives at the expression dp e ipn (26) and after a Fourier transform, it finally leads to [47] We argue that this exact asymptotic result for the occupation function n k (ζ) allows for an intuitive hydrodynamic interpretation.Indeed consider again the tight-binding fermionic Hamiltonian (9).The continuum limit is then attained by dividing the chain into a set of equally-spaced interval of size ∆x = M a, where a is the lattice spacing and M 1 is the number of sites in each interval ∆x.By taking the scaling limit where a, ∆x → 0 at fixed ∆x/a = M , the lattice site j ∈ Z is replaced by the continuous variable x ≡ aj ∈ R, and the Hamiltonian ( 9) is written as [39,95] with continuous fermionic operators ĉ † x ≡ ĉ † ja = ĉj , ĉx = (ĉ † x ) † .Furthermore, around each coarse-grained point, the system is assumed to be in an eigenstate of a system in a periodic box of size ∆x.Therefore, the Hamiltonian (28) can be diagonalized in Fourier space as light cone region particle evoling along Left -The initial state (2) can be described in terms of the local fermionic occupation function (32), here illustrated in the position-momentum plane, i.e., in the phase-space of the free Fermi gas.Right -At times t > 0, the evolution of the occupation function is obtained from the equation of motion of each noninteracting mode and it displays a non-trivial phase-space configuration in the region |x| t (light cone region), where the correlations during the DW melting process are spread from the junction, see main text.
where now the local excitations in x are created by b † k,x , defined through At this point, one can define the local fermionic occupation function for the DW initial state as For free fermionic systems, one can notice from Eq. ( 31) that the local occupation function of fermionic modes corresponds to the Wigner function of the system.In the hydrodynamic limit that we are considering, the Wigner function takes only positive values n k (x) ∈ [0, 1] and it carries the physical interpretation of probability of finding a particle of momentum in a coarse-grained position (x, k) of the phase space.It is then easy to see that the macrostate corresponding to the initial state ( 2) is given by The form of Eq. ( 32) for the hydrodynamic description of the DW state is self-evident as it fills the l.h.s. of the system with modes −π k π, leaving the right part empty.Due to the non-interacting nature of the particles, each mode of momentum k, initially at position x 0 , evolves ballistically with constant velocity v(k) = sin k according to the equation of motion It follows that the time-evolved fermionic occupation function takes the form The latter can be also interpreted as the solution of the following Euler equation which is satisfied by the Wigner function at lowest order in ∂ x and ∂ k derivatives.
Moreover, since our initial state (32) has zero entropy and since the Euler equation ( 35) is known to preserve the entropy of a given initial state at any time during the time evolution, the solution for n k (x, t) can be reconstructed by following the time evolution of its contour.For intuitive reasons, the latter is typically referred to as Fermi contour Γ t and it contains the information about the local Fermi points k ± F of the model at any space-time position (x, t) [17,40]: By noticing that the Fermi contour in the initial state is given by those particle living at the interface at x = 0, the expression for k ± F (x, t) is easily obtained from Eq. ( 33) (37) Using this, one can write the time-evolved fermionic occupation function as hence recovering the result in Eq. ( 27) via hydrodynamic arguments.From Eq. ( 38) one can obtain the asymptotic behavior of conserved charges and associated currents straightforwardly.Indeed, given a local conserved charge Q of the model ( 9), the expectation value of the charge density is given as the integral sum over the available modes at space-time position (x, t) weighted with the single-particle eigenvalue h Similarly, the associated current is obtained as  2).The full lines show the scaling function in Eq. ( 22) (resp.Eq.( 24)) for the magnetization (resp.the spin current).Bottom panels -Same plots as function of ζ = x /t.The data collapse at different times signals the presence of ballistic transport in our quench setting.
For instance, by setting h one can recover the fermionic density n x (t) ≡ m x (t) + 1/2 and the associated current in Eqs. ( 22)- (24).Notice that using Eq. ( 39)-( 40), the conservation law for each conserved quantity follows directly from the evolution of the filling function in Eq. (35).In Fig. 4, we compare the asymptotic results for the profiles of m x (t) and j M x (t) with exact numerical calculations for a spin chain of N = 600 sites, finding an excellent agreement and a perfect data collapse in the scaling variable ζ = x /t.Before turning to the analysis of the interacting spin chain, we wish to conclude this paragraph with a final remark.From both analytical and numerical calculations, we clearly observe a dependence of the asymptotic fronts on the scaling variable ζ = x /t, rather than on x and t separately.This feature is related to the Euler hydrodynamic description of the DW melting problem, and for integrable models, it gives rise to ballistic transport.Along a ray of fixed ζ, the system keeps its steady-state configuration during the whole time evolution.For ζ → ±∞ (correspoding to either |x| → ∞ or t → 0), one inevitably moves outside the correlated region and hence recovers the initial DW configuration of the spin chain.Indeed, in our problem, the modes that are responsible for the non-equilibrium dynamics are originated at the junction x = 0 and spread throughout the left and right part of the spin chain with a finite velocity, bounded by the maximum value max k [v(k)] = 1.Therefore, the system develops a non-trivial profile in the so-called light cone region −t x t, which is determined by the equation of motion of the fastest excitations k = ±π/2, while Figure 5 The DW melting problem in the x-t plane.The initial junction at x = 0 acts as a source of propating particles, leading to the gradual melting of the initial ordered chain inside the light cone region −t x t.Outside the correlated region, the system keeps its initial DW configuration.Each ray of constant ζ = x /t corresponds to a specific NESS for the dynamics (as shown in the figure by the uniform coloring of rays inside the light cone).We illustrate this feature by plotting the magnetization profile (22) as a color plot .
it keeps its original configuration in those space-time positions where the propagating modes have not arrived yet.In general, a light-cone effect is expected whenever a maximal velocity of propagation exists (e.g.due to the Lieb-Robinson bound [98]).We show this by plotting the profile of magnetization (22) in a (x, t) plane, see Fig. 5.

Interacting case
We now turn to the analysis of transport during the DW melting of the interacting spin chain (1) , which we re-write for convenience As anticipated, we are interested in the regime |∆| < 1, where it is customary to write ∆ = cos(γ).We further focus on the rational case, i.e., on those values of γ that can be written as ratio γ = πQ/P with Q, P two-coprime integers 1 Q < P .
In this case, the interaction γ admits a continued fraction representation where {ν 1 , . . ., ν q } is a set of numbers satisfying ν 1 , . . ., ν q−1 1 and ν q 2. In the limit N → ∞, the model ( 42) is exactly solved by means of Thermodynamic Bethe Ansatz (TBA), see e.g.Refs.[99,100].In particular, for large N and according to the string hypothesis [99], the excitation spectrum of the spin chain is described by different species of quasiparticles, generically referred as strings.The total number δ of strings species is determined by the interaction parameter γ through [8,17] As we will see in the following paragraphs, the asymptotic physics of the DW melting in the interacting case is found to be qualitatively similar to the non-interacting one.

Thermodynamic Bethe ansatz solution
In the following, we provide a short summary of the TBA solution of the model ( 42) tailored for the introduction of the objects which are needed for the DW melting problem discussed below.The interested reader can find a comprehensive treatment e.g. in Refs.[99,100].
In the presence of interaction ∆ = 0, the diagonalization of the spin chain (42) can be performed via Bethe ansatz.In particular, in the Hilbert space sector with M spins up and assuming periodic boundary condition for the spin chain, one can write down the exact eigenstates of the model in terms of some complex parameters {λ j } M j=1 (usually called rapidities), which generalize the concept of particle momenta of a free Fermi gas to the interacting case.The rapidities λ j are obtained through the solution of non-linear algebraic equations implementing nontrivial quantization condition for the interacting model, see e.g.Refs.[8,17,99].As mentioned above, for sufficiently large N , the rapidities of the model are arranged in symmetric patterns around the real axis called strings S j Each string specifies a quasiparticle species in the model.The parameter λ α j ∈ R is called string center and the index α = 1, . . ., j identifies a specific quasiparticle belonging to the string S j .The number j determines the total number of quasiparticles of a given species.The quantities j and F j can be determined using Bethe ansatz.Their precise expressions can be read in e.g.Refs.[8,17,99].
In the limit N → ∞ the spectrum of the model becomes densely populated and the state of the system can be suitably described in terms of a spectral distribution of quasiparticles determined by the string center only.The spectral function ρ (j) (λ) can be obtained directly from the solution of the integral equation where s j = ±1 is called string sign and a j (x), T j,i (x) are interaction kernels, see e.g.Refs.[8,17,99] for their expressions.The knowledge of ρ (j) (λ) fully characterize the thermodynamic properties of the zero-entropic states of the spin chain (42).
For later purposes, we introduce also the occupation function for zero-entropic states as For more generic states having a non-zero entropy (e.g. a thermal state), one has to introduce a density distribution for the unoccupied rapidities (or holes) ρ (j) h (λ) and complement the Bethe ansatz solution with additional thermodynamic arguments that relate the functions ρ (j) (λ) and ρ (j) h (λ).

Generalized Hydrodynamics
In the presence of an inhomogeneity (such as the kink state (2)), the Bethe ansatz solution of the model ( 42) breaks down due to the absence of translational invariance.Nevertheless, by considering a hydrodynamic limit for the interacting spin chain similarly to Sec. 3.1.3,one can describe the asymptotic properties of the model as a collection of periodic boxes of size ∆x, each containing a large number of lattice sites.In this way, the local properties of the initial non-homogeneous state can be determined in terms of a local occupation function n (j) 0 (x, λ), obtained from the TBA solution of the model within the cell ∆x.The large-scale dynamics of the occupation function is then established by the GHD [2,3] This set of equations has the same structure of Eq. ( 35) for the non-interacting model, but it is characterized by an effective velocity v eff dressed by the interactions, which depends self-consistently on the macrostate n (j) t (x, λ).A formal solution of the GHD equations ( 49) is obtained with the method of characteristics, yielding where We mention also Ref. [101], where a geometric approach for the solution of GHD equations has been derived.In general, the solution of Eq. ( 49) is obtained by a numerical implementation.However, in the specific case of non-homogeneous initial macrostates that are locally described by a fully-polarized spin configuration (as it is the DW initial state under analysis), the Bethe ansatz solution of the model greatly simplifies and analytical solution of the GHD can be determined.

Occupation function and effective velocity of fully-polarized states
As anticipated, in the specific case where the local properties of the spin chain in the cell ∆x are described by the ensemble where h(x) is an external magnetic field, the local macrostates n (j) (x, λ) can be analytically determined, see Refs.[8,17,99].The expressions for n (j) (x, λ) are not very instructive and therefore are not reported here.The important information that they reveal is that the occupation functions of strings j = 1, . . ., δ − 2 are exponentially suppressed by the magnetic field while the strings δ − 1, δ are not.In the limit of strong magnetic field where depending on the sign of h(x), the expression for n (j) (λ, x) drammatically simplifies and reads as signaling that only the two largest strings j = δ − 1, δ are responsible for the thermodynamic properties of the fully-polarized cell ∆x.It is then easy to see that the DW initial state for the interacting chain is described by the occupation functions From Eq. ( 49), one can notice that the GHD evolution of the state ( 55) is also determined by the sole behavior of the strings j = δ − 1, δ.Furthermore, a careful inspection of the TBA equation of the model for fully-polarized states reveals that and the time-evolved occupation function We refer to Refs.[8,17,38] for details on the derivation of these results.In Eq. ( 56), we introduced ζ 0 ≡ sin(γ)/ sin(π/P ) (58) and the quantity which, up to the string sign s δ , is equal to the bare physical momentum p δ (λ) of the string, see e.g.Refs.[8,17,38] for details.
Eqs. ( 56)-( 57) allow for the derivation of analytical results for the transport properties of the interacting chain (1) during the DW melting, as discussed in the next paragraph.

Analytical results
As already noticed for the non-interacting case, the solution (57) for the time-evolved occupation function is redundant due to the zero-entropy condition preserved by the GHD equations ( 49) at any time during the melting dynamics.Therefore, in analogy with the non-interacting case, we define a local Fermi contour Γ t (∆) generalizing Eq. ( 36) as [17] Γ with Fermi rapidities λ ± F at space-time position (x, t) obtained by the solution of the equation of motion By noticing that the function k δ (λ) is monotone in the interval [−π/P, π/P ] [8,17,38], one can solve Eq. ( 61) with x 0 ≡ 0 as The other root k δ (λ + F ) comes from a Fermi rapidity λ + F → ∞ initially located at x 0 < 0, corresponding to the momentum Similarly, one can show that the other string, say δ − 1, is related to Fermi points Both strings j = δ, δ − 1 contribute with equal weight δ−1 = δ = P/2 and, for P = 2, they properly reproduce the situation in absence of interactions [38].Notice that the structure of the Fermi contour Γ t is very similar to that of the non-interacting GHD GHD GHD GHD GHD GHD Figure 6 Asymptotic profiles of magnetization (left panels) and of spin current (right panels) as function of the ratio ζ = x /t, plotted for different values of interactions.From top to bottom: γ = 3π /7, π /3, π /4 corresponding to the anisotropy parameter ∆ = 0.222, 0.5, 0.707.The symbols show the numerical data obtained with tensor network simulations for a spin chain of N = 150 sites, while the full lines show the analytical GHD prediction in Eqs. ( 66)- (67).The agreement of the curves is seen extremely good, especially at large times where the convergence towards the GHD is improving.We do observe some small finite-size effects at the edges of the light cones that can be minimized by considering larger spin chains.Though, this task is found non-trivial due to the increasing of entanglement during the dynamics, which significantly slow down tensor-network based simulations, see next section.Finally, we notice that the light cone region shrinks by increasing the interaction (from top to bottom panels), see Eq. (65).In each panel, we marked the light cone positions by a dashed vertical axes.spin chain, although the presence of interactions lead to a shrinking of the light cone region from −1 x /t 1 at ∆ = 0 to [8,17,38] light cone region The charge profiles during the DW melting process follow straightforwardly.In particular the profile of magnetization for |ζ| √ 1 − ∆ 2 is obtained as the weighted integral sum over the available quasiparticles as [8,17] Outside the light cone , the system keeps its initial value of magnetization m = −1/2 (resp.m = 1/2).Similarly, the spin current profile for |ζ| √ 1 − ∆ 2 is obtained as [8,17] (see also Ref. [102][103][104] for a rigorous derivation of (67) in spin chain models) and it vanishes outside the light cone.One can notice that for Q = 1 and P = 2 corresponding to ∆ = 0, the results in Eqs. ( 22)-( 24) are recovered.In Fig. 6, we show the comparison of the analytical GHD prediction in Eqs. ( 66)-( 67) against tensor network numerical simulations for the interacting spin chain, performed with the open-source libray iTensor [105].The agreement of the curves with the numerical data is extremely good.

Quantum description and entanglement dynamics
The goal of this section is to complement the discussion on the transport properties of the DW melting with a study on entanglement.As we shall review in the following paragraphs, an ab-initio characterization of entanglement with exact lattice calculations is very demanding also in absence of interaction, the case ∆ = 0, due to the inhomogeneous and non-equilibrium character of the problem.In addition, even the hydrodynamic evolution established by the local occupation function in Eqs. ( 35)-( 49) is not sufficient for the study of quantum correlations among different cells, hence of entanglement.The reason of this failure is rooted in the assumption that underlies the hydrodynamic limit considered so far.In particular, in deriving the hydrodynamic picture, we assumed that each cell is described by a local density matrix of the form of a GGE, i.e., ˆ t (x) ∝ exp[− i β i (x, t) Qi ] in terms of some local Lagrange multipliers β(x, t) associated with each conserved quantity, see e.g.
Refs. [3,5] for more details.However, by doing this, long-range quantum coherent effects among different coarse grained points are washed out, resulting in the vanishing of equal-time correlation functions and of zero-temperature entanglement.
Nevertheless, a possibility to restore these missing quantum effects is given by the novel framework of Quantum Generalized Hydrodynamics (QGHD) [17,[36][37][38][39][40][41]43].According to this theory, the relevant processes at low energies are in the form of particle-hole excitations generated near the local Fermi points and can be described by a non-standard Luttinger liquid living along the evolving Fermi contour.In the following, we briefly revisit the QGHD framework and detail the solution for the DW melting problem.For the sake of clarity, we shall treat first the case ∆ = 0 and afterwards extend the discussion on entanglement to the interacting spin chain.

Case: ∆ = 0
In absence of interactions, it is well known that the entanglement properties of the lattice model ( 9) can be related to the spectrum of the two-point correlation function (18), see Ref. [106][107][108][109][110][111] for a discussion.In particular, given a bipartition of the spin chain as A ∪ B = [− N /2 + 1, j] ∪ [j + 1, N /2], one has the expression for the Rényi entropies and for the entanglement entropy Here, Λ l are the eigenvalues of the two-point correlation function (18) restricted to the subsystem A, i.e., G A ≡ [G n,m ] n,m∈A .However, such a direct approach is very demanding and often not possible in inhomogeneous quench settings as the one that we are considering.In some special cases, some results have been obtained following this strategy (see e.g.Refs.[112][113][114][115][116]) but these are based on non-trivial lattice calculations and therefore will not be discussed in this context.

Re-quantization of the Fermi contour
In this short review, we rather proceed by considering a re-quantization of the Euler hydrodynamics of Sec. 3 including linear quantum fluctuations at the edges of the local Fermi points as and define the excess density simply as δn t (x) = δ kt (x)/π.Using standard bosonization arguments [117][118][119], one can then relate the fluctuating field δn t (x) to the vertex operators of an effective field theory arising at large scales and at low energies.To this end, we express δn t (x) as where φt (x) is a height field encoding the long-range density fluctuations along the one-dimensional spin chain, see e.g.Ref. [120][121][122][123]. From Eq. ( 70), one can obtain the leading order term of the Haldane harmonic-fluid expansion for the time-dependent lattice fermions as [36,39,43] ĉ † x (t) ∼ C(x, t)e −iϕ(x,t) : e i θt(x) : + h.o.c.(72) with dimensionful non-universal amplitude C(x, t) and semi-classical phase ϕ(x, t) that are unimportant for our scopes, see e.g.Refs.[36,55,121,123] for a discussion.
Here, : • : denotes the normal ordering of fields and θt (x) is a phase-fluctuating field satisfying Higher order terms are obtained from vertices with higher scaling dimension but their contribution is negligible in the low-energy regime.It is then customary to express the fields φt (x) and θt (x) in terms of two chiral bosons φ(±) [118] φt (x) ≡ 1 describing left-and right-moving propagating sound waves along the Fermi contour.The parameter K t (x) is usually called Luttinger parameter and is related to the local compressibility of the off-equilibrium quantum fluid.In absence of interactions, free fermionic limit: irrespective on the space-time position.At this point, by plugging the expansion (74) in ( 9) and retaining only quadratic terms, one arrives to the Luttinger liquid Hamiltonian where the sound velocity v t (x) ≡ sin k F (x, t).In terms of the chiral fields, Finally, from the expression (76) for the Luttinger liquid Hamiltonian, one obtains the action [38,55,56,[121][122][123] S = det(g) dt dx 2πK t (x) where the indices a, b = x, t and we restored the dependence on K t (x) ≡ 1 for future convenience.The action in (78) describes a free massless compact boson with space-time dependent coupling K t (x) and time-dependent non-flat 2-dimensional metric tensor g ab , whose line element reads as It is easy to show that this metric can be mapped to a flat one with a simple change of coordinates where dx → dx = dx/v t (x), see e.g.Ref. [38,39,55,56,[121][122][123] for details.Hence, the action (78) displays conformal invariance in those settings where K t (x) = const, as it is the case for a non-interacting spin chain (cf.Eq. ( 75)).In the following, we shall make use of tools stemming from conformal field theory (CFT) to establish the universal behavior of the entanglement entropy during the DW melting problem.

Universal behavior of entanglement
The universal behavior of Rényi entropy for a bipartition A ∪ B ≡ [−∞, x 0 ] ∪ (x 0 , +∞] with cutting point in a real-space position x 0 can be established using the known relation [124][125][126] obtained from the path-integral representation of the operator [ρ x0 ] p on a p-sheeted Riemann surface.Here, T (p) t (x 0 ) is a boundary field known as twist field connecting local operators among the p-copies of the model across the branch-cut at x 0 , see e.g.Refs.[125,126] for a discussion.Crucially, in our model, the twist field is a primary operator of the CFT along the Fermi contour with scaling dimension and it allows for a chiral decomposition in terms of the chiral components τ (p) t (x 0 , ±), each with scaling dimension d p /2.Notice that these objects are highly non-trivial if compared to their equilibrium counterparts, due to their time dependence and to the non-homogeneous background over which expectation values are evaluated.However, a natural description for these boundary fields is attained by moving from space-time coordinates (x, t) to a coordinate s ≡ s(x, t) along the Fermi contour Γ t at time t.By doing this, Eq. ( 80) becomes where s ± are the boundary points of the subsytem A along Γ t , see the discussion below.Notice that in terms of the coordinate s, sound waves can be described by a unique component circulating along the contour, see e.g.Refs.[39,43].The expectation value in the last equation is then obtained by standard boundary CFT and reads as The coordinate s(x, t) is an isothermal coordinate for the space-time dependent metric (79).Its explicit expression is generically out-of-reach, excluding outstanding cases such as Ref. [36] and Refs.[56].For generic interacting problems, one typically finds an isothermal coordinate s(x, 0) for the initial inhomogeneous configuration as done in Refs.[120][121][122][123] and evolves the initial correlations in time according to the equation derived in Ref. [37].For free systems, the isothermal coordinate s(x, 0) is a comoving coordinate since quantum fluctuations in the initial state do not evolve and are only transported along the contour Γ t during the dynamics, see Refs.[17,[39][40][41]55] for examples.
The DW melting problem falls in the simplest category.An isothermal coordinate for this quench is simply given by [55,56] It follows that the boundary points s ± for a real-space cut at position x 0 are given by and, by elementary algebra, one obtains directly from Eq. ( 84) the result [55] S(p) Despite its simplicity, Eq. ( 85) allows for a non-trivial interpretation that we wish to comment.In fact, the change of coordinates (x, t) → s ≡ k F (ζ) implements a mapping between real-space correlations and momentum-space correlations of the spin chain during the melting dynamics.At t = 0, we prepare our system in a product state, i.e., with exact zero entanglement for any cutting position x 0 .However, the same does not hold true if one considers the entanglement of a bipartition in momentum space, for which standard boundary CFT can be applied.During the time evolution, momentum-space entanglement is gradually transported to real-space through the ballistic propagation of particles, resulting in the logarithmic growth of Eq. (87).In this sense, the isothermal coordinate s(x, t) may be viewed as a bridge between real and momentum space, and therefore it establishes the amount of spatial correlations that are generated via propagation at a given space-time point (x 0 , t).
We illustrate this procedure in Fig. 7.

Short distance regularization and asymptotic results
At this point, we recall that Eq. ( 87) provides only the universal contribution to the Rényi entropies and it has to be complemented with a non-universal regularization at short distances [55].The latter, for homogeneous equilibrium models, typically introduces an additive constant and therefore is often neglected or estimated as fitting parameter.Though, in out-of-equilibrium and non-homogeneous situations, its effect is more pronounced and significantly modify the space-time dependence of Rényi entropies in Eq. ( 87).Hence, by dimensional analysis, we can write the regularized Rényi entropy as where ε x0 (t) is a local UV cutoff.On physical grounds, one expects that the shortdistance regularization is set by the microscopic scale of the model within the cell x 0 , which is the inverse local density, and reads as From exact Fisher-Hartwig calculations made for the microscopic XX model [127,128], one can refine the ansatz in Eq. ( 89) and obtains the correct value of cutoff as [17, 39-41, 43, 55] with C a known non-universal amplitude, see Refs.[127,128] for its expression.Using Eq. ( 90) in Eq. ( 88), after simple algebra one obtains the final result where the additive constant κ p = − (p+1) /12p log[C/2], for instance κ 1 = 0.4785 [128].By setting x 0 = 0 and p = 1 in Eq. ( 91), we see that the half-system entanglement entropy displays a logarithmic growth in time as [17,40,55] This results can be alternatively viewed as a manifestation of Calabrese-Cardy formula [124] (see also Ref. [129] by Holzhey-Larsen-Wilczek), since the size of the correlated region is growing linearly with time in our problem.In Fig. 8, we show the entanglement profiles ( 91) and the half-system entanglement entropy (92) against exact numerical data and we observe a perfect matching of the curves with numerics.

Case: −1 < ∆ < 1
As anticipated, in generic interacting integrable models the Luttinger parameter K t (x) is non-constant and non-homogeneous and, therefore, the conformal invariance of the action ( 78) is generically lost.Its value can be determined (as at equilibrium) from the value of the dressed magnetization evaluated at any of the local Fermi rapidities λ ± F (x, t) in Eq. ( 60), see Ref. [38] for a discussion and e.g.Ref. [100] for details on the calculation.However, the very peculiar features about the Bethe ansatz of the DW state (cf.Refs.[8,17,99] and Sec.3.2.3)that led to the very outstanding possibility of finding an analytical solution for the GHD equation of the model (cf.Refs.[8,17,38] and Sec.3.2.4),permit the analytical calculation of the Luttinger parameter, which for both the strings j = δ − 1, δ equals to Remarkably, the Luttinger parameter for the DW melting problem is constant and depends only on the interaction parameter (actually only on the denominator, recall that γ = arccos(∆) = πQ/P ), displaying a peculiar fractal dependence [38].Consequently, conformal invariance in Eq. ( 78) is not broken even at finite interactions and one can obtain the Rényi entropy of a bipartition with cutting point x 0 employing the same techniques of Sec.4.1.In particular, one can still write the relation and, by scaling arguments, one obtains [38] S (p) x0 (∆) = p + 1 12p log(t) + f p ( x /ζ0t) + κ p (∆) (95) where f p (ζ/ζ 0 ) is a scaling function that has been conjectured in Ref. [38] and κ p (∆) is a non-universal additive constant that we are currently unable to analytically determine.From the numerical analysis of Fig. 9, we observe that the qualitative behavior of the entanglement for the interacting chain |∆| < 1 is not limited to the logarithmic half-system growth (cf.Eq. ( 95) and ( 92)) but instead applies to the whole entanglement profiles, modulo a rescaling by the interactions of the light cone, as pointed out in Eq. ( 65).

Summary and concluding remarks
In this short review, we revisited the out-of-equilibrium physics arising during the unitary evolution of a one-dimensional spin chain prepared in a domain wall configuration.For this problem, we focused on transport properties (i.e., characterization of charge and current profiles) and on the study of entanglement entropy.In Sec. 3, the first aspect has been first addressed for a noninteracting spin chain by exact lattice methods, from which a hydrodynamic approach has been later developed.We then extended such a hydrodynamic approach to the interacting case by means of Bethe ansatz techniques and through the use of generalized hydrodynamics, both of which are briefly commented on in the text.We derived exact results for charge and current profiles and showed how these are found in perfect agreement with the numerical data, obtained by exact methods or by the use of tensor networkbased numerical simulations.
In the second part of the review (Sec.4), we revisited the study of entanglement dynamics by means of quantum generalized hydrodynamics, which consists of the use of CFT tools arising from the requantization of the hydrodynamic background through a Luttinger liquid.Analytical results for Rényi entropy profiles at different times have been derived and compared with numerical data.
As commented in the introduction, the domain wall melting problem has been the subject of numerous studies e.g.[17, 39-42, 47, 53, 55-57, 59, 61, 62, 71, 112, 113], from which numerous analytical results have been collected over years.Besides being one of the very rare cases of non-trivial systems where such solutions can be derived, the study of this setting is also extremely useful as basis for investigating the dynamics of other quench problems with generic integrable models.Hence, the purpose of this review is twofold: first, to concisely provide a summary of available results for domain wall melting, and second, to give an overview accessible to a broad community of the hydrodynamic approach to integrable systems via a thorough study of a prototypical model, which has recently attracted a great deal of attention.

Figure 2
Figure2Magnetization profiles as function of the lattice site j, plotted at different times.The symbols show the numerical data, obtained from the exact diagonalization of the lattice Hamiltonian with N sites and by a further Trotter evolution of the initial two-point function, see e.g.Refs.[39,40] for details.The full lines show the exact result of Eq. (15).

Figure 3
Figure 3 Hydrodynamic description of the DW melting problem for the non-interacting spin chain.Left -The initial state (2) can be described in terms of the local fermionic occupation function(32), here illustrated in the position-momentum plane, i.e., in the phase-space of the free Fermi gas.Right -At times t > 0, the evolution of the occupation function is obtained from the equation of motion of each noninteracting mode and it displays a non-trivial phase-space configuration in the region |x| t (light cone region), where the correlations during the DW melting process are spread from the junction, see main text.

Figure 4
Figure 4 Asymptotic behavior of the charges and current profiles during the DW melting for the noninteracting spin chain.Top panels -Asymptotic magnetization (left) and spin current (right) profiles as function of the rescaled position x/N , for different values of time.The symbols show the numerical data obtained for a chain of N1 via exact diagonalization and subsequent Trotter evolution of the two-point function (see also the caption of Fig.2).The full lines show the scaling function in Eq. (22) (resp.Eq.(24)) for the magnetization (resp.the spin current).Bottom panels -Same plots as function of ζ = x /t.The data collapse at different times signals the presence of ballistic transport in our quench setting.

Figure 7
Figure 7  Illustration of the QGHD framework for the DW melting problem.Left -At t = 0, we include linear quantum fluctuations δn(x) on top of the initial Fermi contour at x = 0 and we map momentumspace correlations into a CFT that lives along a unit circle of length 2π.Right -At t > 0, through the ballistic propagation of particles from the junction, the system develops real-space entanglement for any bipartition with cutting point x 0 inside the light cone.The amount of entanglement in real space is then established by the local Fermi points at (x, t) (green circles) and simply translates into the entanglement of the interval [s − , s + ] for the CFT along the Fermi contour, thanks to the non-interacting nature of the problem.

Figure 8
Figure8Left -Entanglement entropy for the DW melting of a non-interacting spin chain, plotted as function of the cutting point x and for different times.Symbols show the numerical data obtained from Eq. (69) after the exact diagonalization of G A while the full lines show the predictions by QGHD of Eq. (91).Right -Half-system entanglement growth during the DW melting function of time, plotted in semi-logarithmic scale.

Figure 9
Figure 9  Entanglement profiles during the DW melting of the interacting spin chain, plotted at different times and for different values of interaction (from left to right: γ = 3π /7, π /3, π /4 corresponding to ∆ = 0.222, 0.5, 0.707).The data is obtained with tensor network simulations.Dashed vertical axes mark the position of the light cone at each time, according to Eq. (65).