Domain wall melting across a defect

A single localised impurity or defect can alter the global structure of a many-body quantum system, as well known from the textbook examples of Anderson orthogonality catastrophe [1] and the Kane-Fisher model [2,3]. In the latter, it has been shown that for repulsive interactions, the electrons are completely reflected by even the smallest scatterer, leading to a truly insulating weak link disconnecting the two halves. Conversely for attractive bulk interactions, the weak link is irrelevant, i.e., it is washed away at large scales. As a consequence free fermions represent the most interesting system in which the defect is marginal and there is a line of fixed points characterised by the defect strength [4,5].

In recent years, the physics of impurities in one-dimensional (1D) free-fermionic systems has been investigated a lot through the lens of entanglement. The marginality of the defect is reflected into a logarithmic scaling of the entanglement entropy with a prefactor that depends continuously on the defect strength [6–20]. Overall, thanks to all these studies nowadays we have a rather complete understanding of the physics of defects in equilibrium free fermionic systems. The same is definitively not true when the free fermionic chain is driven out of equilibrium; in fact, in spite of several works about the non-equilibrium behaviour across one defect (see, e.g., refs. [21–31]), a complete understanding is still far because of the many different ways of driving a system away from equilibrium.

In this manuscript, we prepare an initial state with a domain wall localised at the defect and we let it melt. Without the defect, this is a protocol that has been studied intensively in the free fermionic literature [32–45] and it represents a case study for the application of hydrodynamics to non-equilibrium quantum systems (recently adapted also to interacting integrable models [46,47] with the generalised hydrodynamics formalism [48,49]). Anyhow, the presence of the defect is expected to alter dramatically the evolution at a qualitative level. The density, the currents, and other local quantities have been characterised in ref. [50] where the emergence of a local non-equilibrium stationary state (NESS) has been rigorously established. However, little is known for the entanglement entropy, whose behaviour is affected, as any other non-local observable, by non-local correlations generated by the defect. Some lattice results have been derived for the domain wall melting with defect in ref. [21], and another important step forward has been done by Fraenkel and Goldstein [51,52] in a slightly different context, but a general scheme to describe non-local correlations is still missing. In this regard, it is very natural to wonder whether even a small defect alters substantially the domain wall melting. As we shall see, this is the case.

We consider a 1D chain of free spinless fermions with 2N sites and with nearest-neighbour hopping with a defect of strength λ located at the center of the system. The Hamiltonian is

$$\begin{equation} \hat{H}=\sum_{i,j=-N+1}^{N} h_{i,j} \ \hat{c}^\dagger_i\hat{c}_j, \end{equation} \tag{ 1 }$$

with

$$\begin{equation} h_{i,j}=-\frac{1}{2}(\delta_{i,j+1}+\delta_{i+1,j}), \quad \forall i,j\neq 0,1 \end{equation} \tag{ 2 }$$

and defect

$$\begin{equation} h_{0,1}= h_{1,0}=-\frac{\lambda}{2}, \quad h_{0,0}= -h_{1,1}=\frac{1}{2}\sqrt{1-\lambda^2}. \end{equation} \tag{ 3 }$$

Here $\hat{c}^\dagger_j$, $\hat{c}_j$ are the creation and annihilation operators of spinless fermions at site j, satisfying $\{\hat{c}^\dagger_j,\hat{c}_i\} = \delta_{ij}$. For $\lambda = 1$ the Hamiltonian (1) reduces to a standard hopping model. The system is initially prepared in the product state

$$\begin{equation} |{\psi_0}\rangle = \bigotimes_{j=-N+1}^{0} |{1}\rangle_j \bigotimes_{j=1}^{N} |{0}\rangle_j, \end{equation} \tag{ 4 }$$

where $|{\zeta=0,1}\rangle_j$ are the eigenstates of the fermionic number operator $\hat{c}^\dagger_j\hat{c}_j$ with eigenvalues $\zeta = 0,1$. In the spin language, the initial state (4) corresponds to the domain wall $|{\psi_0}\rangle = |{\dots\uparrow\uparrow\downarrow\downarrow\dots}\rangle$. Thereafter, our discussions can be extended to the XX spin chain with defect [8], obtained from the model (1) through a Jordan-Wigner transformation. For t > 0, the state (4) is unitarily evolved with Hamiltonian (1), $|{\psi_t}\rangle = e^{-{\rm i} t\hat{H}}|{\psi_0}\rangle$.

The structure of the defect (3) does not spoil the exact solvability of the free fermionic model, which has spectrum $\omega_q = -\cos(k_q)$, $k_q = \pi q/(2N)\,(q = 1,\dots,2N)$, for any value of $\lambda\in(0,1]$ [10,26]. Moreover, the eigenstates of $\hat{H}$ can be related to the eigenstates $\Psi_{q,2N}(j) = \sin(k_q j)/\sqrt{N}$ of (1) in the absence of defect $(\lambda = 1)$ as [26]

$$\begin{equation} \Psi^\text{def}_{q,2N(j)}=\Theta(-j)\alpha^+_q(\lambda)\Psi_{q,2N} + \Theta(j)\alpha^-_q(\lambda) \Psi_{q,2N}, \end{equation} \tag{ 5 }$$

with $\Theta(j)$ the Heaviside step function and coefficients $\alpha^\pm_q(\lambda)=[1\pm(-1)^q\sqrt{1-\lambda^2}]^{1/2}$. For N → ∞, this eigenproblem reduces to a scattering of plane waves across a localised defect, i.e.,

$$\begin{equation} \Psi^\text{def}_{k,\infty}(j)\propto\Theta(j) \lambda \ e^{{\rm i} k j} + \Theta(-j)(\sqrt{1-\lambda^2} \ e^{-{\rm i} k j} + e^{{\rm i} kj}), \end{equation} \tag{ 6 }$$

with transmission probability ${T}(\lambda)\equiv \lambda^2$ and reflection probability ${R}(\lambda)\equiv 1-\lambda^2$. These parameters do not depend on the momentum k of the scattered particle and so the defect (3) is also known as conformal defect [8,10,21].

Exact asymptotic results for the charges profiles can be obtained in the hydrodynamic limit $N \rightarrow \infty$, $j \rightarrow \infty$, $t \rightarrow \infty$ at fixed j/t. The lattice index j and the quantised momenta k_q are replaced by continuous variables for the position $x = j a\in\mathbb{R}$ (a is the lattice spacing) and for the momenta $-\pi\leq k \leq \pi$. In such scaling limit, the essential information on the initial state (4) is retained by the local fermionic occupation

$$\begin{equation} n_0(x,k)=\begin{cases} 1, \qquad \text{if $x\leq 0$\enskip and\enskip $-\pi\leq k \leq \pi$};\\ 0, \qquad\text{otherwise}. \end{cases} \end{equation} \tag{ 7 }$$

In the absence of defect $(\lambda = 1)$, the evolution of the occupation function is given by the Euler equation

$$\begin{equation} (\partial_t + \sin k \partial_x) n_t(x,k)=0, \end{equation} \tag{ 8 }$$

which can be simply solved as

$$\begin{equation} n_t(x,k) = n_0(x-t\sin k, k), \end{equation} \tag{ 9 }$$

 i.e., it stays equal to 0 and 1 and only the Fermi contour separating the two values evolves. Intuitively, this solution encodes the fact that each non-interacting particle moves along the ballistic trajectory with constant velocity $v(k) = \sin k$. Instead, for $\lambda\neq 1$, a particle of momentum k > 0 traveling from x < 0 is scattered by the defect in such a way that it is reflected with probability $R(\lambda)$ and transmitted with probability $T(\lambda)$. Accordingly, the time-evolved occupation function in the presence of the defect takes the form [50]

$$\begin{align} n^{(\lambda)}_t(x,k) &= \lambda^2 \Theta(x) n_t(x,k) \nonumber\\ &\quad+ \Theta(-x)\left[ (1-\lambda^2)n_t(-x,-k) + n_t(x,k)\right], \end{align} \tag{ 10 }$$

as illustrated in fig. 1. In our notations $n_t^{(\lambda\equiv1)}= n_t$ in eq. (9). The occupation function (10) gives us access to the asymptotic profiles of conserved charges as elementary integrals over the modes k, properly weighted with the single-particle eigenvalue of the associated charge [32,33]. For instance, the particle density profile for $0 < x\leq t$ is

$$\begin{equation} n^{(\lambda)}_t(x)=\int_{-\pi}^\pi \frac{\mathrm{d} k \ n^{(\lambda)}_t(0< x\leq t,k)}{2\pi}=\lambda^2 \frac{\arccos(x/t)}{\pi}. \end{equation} \tag{ 11 }$$

For $-t\leq x < 0$, the profile is obtained via particle-hole symmetry and reads $n^{(\lambda)}_t(x)=1-\lambda^2 \arccos(|x|/t)/\pi$. Outside the correlated region, i.e., for $|x| > t$, the systems keeps its initial configuration with constant density $n^{(\lambda)}_t = 1$ ($n^{(\lambda)}_t = 0$) on its left (right) part. In fig. 2, we compare numerical results obtained from exact lattice calculations with the hydrodynamic result (11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The crucial observation now is that the local occupation function (10) for $\lambda\neq1$ assumes values which are different from 0 and 1. Consequently, the local entropy is non-vanishing, resulting in an extensive entanglement. Physically this phenomenology roots back to the correlation between reflected and transmitted modes generated at the defect.

We now move to our main goal which is characterising the entanglement dynamics. Specifically, we focus on a bipartition of the system $A\cup B$ with a reduced density matrix $\hat\rho_t(A) = {\rm tr}_B|{\psi_t}\rangle \langle{\psi_t}|$. The n-Rényi entropy is

$$\begin{equation} S_n(A,t)=\frac{1}{1-n}\log {\rm tr}[\hat\rho_t(A)^n] \end{equation} \tag{ 12 }$$

that provides the entanglement entropy in the limit $n \rightarrow 1$, i.e., $S_1(A,t) = -{\rm tr}[\hat\rho_t(A)\log\hat\rho_t(A)]$.

An ab initio description of the entanglement dynamics is very demanding even in the absence of defect, due to the non-equilibrium and non-homogeneous character of the quench problem under analysis. However, the asymptotic behaviour of entanglement can be determined with hydrodynamic arguments as follows. First of all, we recall the definition of the local Yang-Yang Rényi entropy [52–55]

$$\begin{align} s_n(x,t) &\equiv \frac{1}{1-n} \int_{-\pi}^\pi \frac{\mathrm{d} k}{2\pi}\log[n_t^{(\lambda)}(x,k)^n \nonumber\\ &\quad + (1-n_t^{(\lambda)}(x,k))^n], \end{align} \tag{ 13 }$$

which is non-zero because $n_t^{(\lambda)}\neq 0,1$. From an entanglement perspective, this entropy measures the correlations between the particles at position x with the ones at $-x$ generated by the scattering at the defect, according to the quasiparticle picture [56,57]. The total entanglement entropy of a region A which is entirely to the right or to the left of the defect (say $A = [-\infty,x_0]$) is then given by (see also [27,48])

$$\begin{equation} S_n(x_0,t)=\int_{A} \mathrm{d} x s_n(x,t) = \frac{{\cal N}_t(A)}{1-n}\log\left[\lambda^{2n}+(1-\lambda^2)^n\right], \end{equation} \tag{ 14 }$$

with ${\cal N}_t(A)$ being the total number of entangled particles in the region A at time t. For example for $A = [-\infty,x_0]$ with x₀ < 0, we have

$$\begin{equation} {\cal N}_t([-\infty,x_0])=\frac{t}{\pi}\left(\sqrt{1-\frac{x_0^2}{t^2}}-\frac{x_0}{t}\arccos\frac{x_0}{t}\right). \end{equation} \tag{ 15 }$$

By setting $x_0 = 0$, eq. (14) predicts a linear growth of entanglement (see also ref. [21])

$$\begin{equation} S_n(0,t)= \frac{t}{\pi(1-n)} \log\left[\lambda^{2n}+(1-\lambda^2)^n\right]. \end{equation} \tag{ 16 }$$

In the absence of defect, previous studies highlighted a half-system entanglement growth $S_n(0,t)\sim (n+1)/12n \log(t)$, arising from subleading contributions [35,45,58]. Interestingly, the entanglement transition from logarithmic to linear law is observed even for values of λ very close to unit, see fig. 3(a) for a comparison with exact lattice calculations. We refer to, e.g., refs. [34,58] for details on the numerical implementation.

Zoom In Zoom Out Reset image size

Equation (14) fails to capture the behaviour of entanglement for a subsystem straddling the defect because it counts also for the pairs of entangled particles which are both in A, but on different sides of the defect. Such over-counting is however easily cured within the quasiparticle picture [56,57]. First, for the case $A = [-\infty,x_0]$ with x₀ > 0, using particle-hole symmetry and $S_n(A,t) = S_n(\bar{A},t)$, we have $S_n(x_0,t) = S_n(-x_0,t)$ where the r.h.s. is in eq. (14). The validity of eq. (14) is tested against exact lattice calculations in fig. 3(b).

For subsystems $A = [x_0',x_0]$ consisting of an interval straddling the defect (i.e., $x'_0 < 0$ and x₀ > 0), we can simply correct the over-counting of eq. (14) by subtracting the doubly counted particles (see fig. 4 for an illustration), resulting finally in

$$\begin{equation} S_n([x_0',x_0],t)=\left|S_n(|x_0|,t)-S_n(|x_0'|,t)\right|. \end{equation} \tag{ 17 }$$

In fig. 5, numerical results for the lattice model are compared with the hydrodynamic prediction in eq. (17), showing an excellent agreement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our results also show the appearance of long-range entanglement in the sense of ref. [52]. Indeed, if we compute the mutual information $I_{A_1:A_2}\equiv S_{A_1}+S_{A_2}-S_{A_1\cup A_2}$ between $A_1 = [-\infty,-x_0]$ and $A_2 = [x_0,\infty]$ we have

$$\begin{equation} I_{A_1:A_2}= 2 S_n([x_0,\infty],t) \end{equation} \tag{ 18 }$$

because $S_{A_1\cup A_2} = 0$ (up to subleading terms, see also the next section). The same remains true if A₁ and A₂ are two symmetric finite intervals. Such large (actually extensive for large t) mutual information is due to the constant presence of shared pairs between symmetric intervals, exactly as in ref. [52].

When the subsystem A is placed symmetrically across the defect, i.e., $A = [-x_0,x_0]$, the entanglement resulting from correlated pairs of particles in eq. (17) vanishes. This is clearly due to the fact that entangled pairs have symmetric positions and so they are either both in A or in the complement. As a consequence, the behaviour of entanglement is entirely due to subleading contributions associated with quantum fluctuations, which cannot be determined with a semiclassical approach. For the homogeneous Hamiltonian $(\lambda = 1)$, a useful way to incorporate quantum fluctuations in our description is established by quantum generalised hydrodynamics [34,45,47,58–62]. According to this theory, the relevant contribution to the entanglement in zero-entropic states is given by linear quantum fluctuations $\delta\hat{n}_t(x)$ at the edges of n_t (x, k), corresponding to the formation of particle-hole pairs near the local Fermi points.

The time evolution of the entanglement entropy is a straightforward but tedious adaptation of the calculations reported, e.g., in ref. [45] for a slightly different situation. We only report here the final result (which indeed coincides with the one in refs. [45,63])

$$\begin{equation} S_n([-x_0,x_0],t)=\frac{n+1}{12n}\log\left[x_0^2(1-{x_0^2}/{t^2})^3\right]+2\kappa_n, \end{equation} \tag{ 19 }$$

with $\kappa_n$ a known non-universal amplitude [64,65], for instance $\kappa_1\simeq 0.4785$. For $t\gg x_0$, eq. (19) predicts a saturation of the half-system entanglement to the value $S_1([-x_0,x_0],t\gg x_0)\approx 1/3\log(x_0)+2\kappa_1$.

Numerical exact calculations for the lattice model reveal a similar behaviour for the half-system entanglement even in the presence of the defect $\lambda\in(0,1)$, see fig. 6(a). We expect that the large-time plateaus in the figure scale like $\log(x_0)$ for large x₀, i.e., we expect

$$\begin{equation} S_1([-x_0,x_0],\infty)\sim \frac{c_\text{eff}(\lambda)}{3}\log(x_0)+\gamma(\lambda). \end{equation} \tag{ 20 }$$

To have an unbiased estimate of $c_\text{eff}(\lambda)$ we proceed as following. We first consider the difference of plateaus reached at fixed λ for different sizes of A, i.e.,

$$\begin{equation} \Delta S^{(\lambda)}_1(x_0,x_0')\equiv S_1([-x_0,x_0],\infty)-S_1([-x_0',x_0'],\infty), \end{equation} \tag{ 21 }$$

(where by "infinite time" we just mean confidently within the plateau); then we take the ratio

$$\begin{equation} r(\lambda)=\Delta S^{(\lambda)}_1(x_0,x_0')/\Delta S^{(1)}_1(x_0,x_0'), \end{equation} \tag{ 22 }$$

that for large $x_0, x_0'$ converges to $c_\text{eff}(\lambda)$ by construction. The perfect collapse in fig. 6(b) of the ratios $r(\lambda)$ for different pairs $x_0, x_0'$ confirms the conjectured behaviour of eq. (20). Moreover, the resulting factor $c_\text{eff}(\lambda)$ is numerically consistent with the effective central charge appearing in the ground-state entanglement of free fermions with defects [8] given by

$$\begin{align} &c_\text{eff}(\lambda)= -\frac{6}{\pi^2} \{(1+\lambda)\mathrm{Li_2}(-\lambda) + (1-\lambda) \mathrm{Li_2}(\lambda) \nonumber\\ &+[ (1+\lambda) \log (1+\lambda) + (1-\lambda) \log (1-\lambda)] \log \lambda \}\,, \end{align} \tag{ 23 }$$

that satisfies $c_\text{eff}(0) = 0$ and $c_\text{eff}(1) = 1$. We believe that, being these logarithmic contributions related to zero-point fluctuations, it should be possible to map explicitly the equilibrium entanglement to the non-equilibrium one. However, this goes beyond the scope of this work.

Zoom In Zoom Out Reset image size

We studied the time evolution of the entanglement entropy in a domain wall melting across a conformal defect. We showed that the pure-system logarithmic growth in time of the entanglement entropy is turned, by the smallest defect, into a linear one with an extensive stationary value corresponding to a non-vanishing thermodynamic Yang-Yang entropy. Furthermore we showed that there are extensive long-range correlations between sites which are mirror images of each other with respect to the defect. None of these effects has an equilibrium counterpart. We also characterised the subleading logarithmic contributions, which are not captured by the quasiparticle picture.

A natural extension of this work could be the study of the dynamics in the presence of multiple defects. In that case, we expect a richer pattern of long-range correlations arising from multiple scattering across the defects, e.g., along the lines of ref. [52]. However, we still do not know how to deal systematically with those effects and how to incorporate them in a quasiparticle picture for the entanglement.

LC is grateful to Viktor Eisler for useful discussions. The authors acknowledge support from ERC under Consolidator grant No. 771536 (NEMO).

Data availability statement: All data that support the findings of this study are included within the article (and any supplementary files).