How order melts after quantum quenches

TL;DR

This paper investigates how long-range order in a quantum many-body system, specifically the spin-1/2 Heisenberg XXZ chain, melts over time after a quantum quench, revealing detailed dynamics through the distribution of the order parameter.

Contribution

It provides the first detailed quantum mechanical distribution function of the order parameter during the melting process after a quench from the Néel state.

Findings

Order melts over time after a quantum quench.

Intermediate times show a broad probability distribution of the order parameter.

Detailed quantum distribution functions characterize the melting process.

Abstract

Injecting a sufficiently large energy density into an isolated many-particle system prepared in a state with long-range order will lead to the melting of the order over time. Detailed information about this process can be derived from the quantum mechanical probability distribution of the order parameter. We study this process for the paradigmatic case of the spin-1/2 Heisenberg XXZ chain. We determine the full quantum mechanical distribution function of the staggered subsystem magnetization as a function of time after a quantum quench from the classical N\'eel state. This provides a detailed picture of how the N\'eel order melts and reveals the existence of an interesting regime at intermediate times that is characterized by a very broad probability distribution.

Figures

Equations

$$P_{A}(m,t)={\rm Tr}\big{[}\rho(t)\delta({\cal O}_{A}-m)\big{]}.$$

$$H_{\Delta }=\sum _j{S}_j^x{S}_{j+1}^x+S_j^y{S}_{j+1}^y+\Delta S_j^z{S}_{j+1}^z.$$

$$P_{\ell }(m,t)$$

$$P_{\ell }(\mu ,t)=\frac{1}{\sqrt{2\pi \sigma (t)}}\exp \left\{-\frac{[\mu -\bar{m}(t){]}^2}{2\sigma (t)}\right\}.$$

$$H_{\Delta }=\sum _j\frac{1}{2}\left[c_j^{†}{c}_{j+1}^{†}+h.c.\right]+\Delta n_j(1-n_{j+1}),$$

$$n_j{n}_{j+1}$$

$$|{\Psi }_t⟩=T\exp \left[-i{\int }_0^{t}d{t}^{\prime }{H}_{\mathrm{MF}}(t^{\prime })\right]|0⟩.$$

$$\rho_{\rm SS}=\sum_{j\geq 0}\rho_{\rm SS}^{(j)}\ ,\ \rho_{\rm SS}^{(j)}={\cal O}\big{(}e^{-\beta_{\rm eff}j\Delta/2}\big{)}\ ,$$

$$P_{2\ell }^{(0)}(m)$$

$$P_{2\ell }^{(1)}(m)$$

$$|{\Psi }_0⟩$$

$$S_{\ell }$$

$$H_{\Delta }^{\prime }$$

$$P_{\ell }(m,t)=⟨{\Psi }_t|\delta (S_{\ell }-m)|{\Psi }_t⟩,|{\Psi }_t⟩=e^{-i{H}_{\Delta }^{\prime }t}|{\Psi }_0⟩.$$

$$S_j^{+}=\prod _i^{j-1}(1-2n_i)c_i,S_j^z=\frac{1}{2}-n_j,$$

$$H_{\Delta }^{\prime }=\sum _j\frac{1}{2}\left[c_j^{†}{c}_{j+1}^{†}-c_j{c}_{j+1}\right]+\Delta n_j\left[1-n_{j+1}\right],$$

$$|{\Psi }_0⟩=|0⟩.$$

$$c_j^{†}{c}_{j+1}^{†}{c}_{j+1}{c}_j$$

$$H_{\mathrm{MF}}(t)$$

$$A_j(t)$$

$$B_j(t)$$

$$h_j(t)$$

$$⟨\mathcal{O}{⟩}_t$$

$$U(t)$$

$$-i\frac{d}{dt}{c}_j(t)$$

$$\frac{d}{dt}\mathbf{d}(t)=i{\mathbb{T}}_{\Delta }[\mathbb{C}(t)]\cdot \mathbf{d}(t),$$

$${\left[\mathbb{C}(t)\right]}_{nm}=⟨0|d_n(t)d_m^{†}(t)|0⟩.$$

$$\frac{d\mathbb{C}(t)}{dt}=i{\mathbb{T}}_{\Delta }[\mathbb{C}(t)]\cdot \mathbb{C}(t)-i\mathbb{C}(t)\cdot {\mathbb{T}}_{\Delta }^{†}[\mathbb{C}(t)].$$

$${\left[\mathbb{C}(0)\right]}_{nm}=\{\begin{array}{ll}{\delta }_{n,m}& \text{if }n>L\\ 0& \text{else}.\end{array}$$

$$P_{\ell }(m,t)$$

$$F_{\ell }(\theta ,t)$$

$$a_j^x\equiv c_j^{†}+c_j,a_j^y\equiv i(c_j^{†}-c_j).$$

$$S_{\ell }=\sum _{j=1}^{\ell }\frac{i}{2}{a}_j^y{a}_j^x.$$

Full text

Introduction. —

A fundamental objective of quantum theory is to determine probability distribution functions of observables in given quantum states. In few-particle systems the time evolution of such probability distributions provides a lot of useful information beyond what is contained in the corresponding expectation values. Recent advances in cold-atom experiments have made possible not only the study of non-equilibrium time evolution of (almost) isolated many-particle systems Langen2015 ; Altman2015 ; kww-06 ; GM:col_rev02 ; hacker10 ; tetal-11 ; shr-12 ; cetal-12 ; langen13 ; MM:Ising13 ; ronzheimer13 ; zoran1 ; lev2017 , but given access to the full quantum mechanical probability distributions of certain observables HLSI08 ; KPIS10 ; KISD11 ; GKLK12 ; Greiner . This provides an opportunity to gain new insights about the coherent dynamics of many-particle quantum systems. One intriguing question one may ask is how order melts, or forms, when an isolated many-particle system is driven across a phase transition. Related questions have been studied in solids, but there one essentially deals with open quantum systems and has access to very different observables, see e.g. Hill ; Cavalieri . The basic setup we have in mind is as follows. Let us consider a system of quantum spins with Hamiltonian $H$ that is initially prepared in a state with density matrix $\rho(0)$. In this state there is long-range order characterized by a local order parameter ${\cal O}=\sum_{j=1}^{L}{\cal O}_{j}$, where $j$ runs over the sites of the lattice and ${\cal O}_{j}$ is a local operator. We are interested in the probability distribution $P_{A}$ of the order parameter ${\cal O}_{A}$ in a contiguous subsystem of linear size $|A|$

[TABLE]

Here $\rho(t)$ is the density matrix of the system at time $t$ and $P_{A}(m,t)$ is the probability that the subsystem order parameter ${\cal O}_{A}$ takes the value $m$ in the state $\rho(t)$. We are interested in cases where the system is initially well ordered at all length scales and $P_{A}(m,t)$ is therefore narrowly peaked around the average ${\cal O}_{A}$. Under time evolution the order melts and at late times and large subsystem sizes $P_{A}(m,t)$ is believed to approach a Gaussian distribution centred around zero Deutsch1991 ; Srednicki1994 ; Rigol2008 ; Rigol2007 ; Cassidy2011 ; Caux2013 ; Ilievski2015 . The question of interest is how $P_{A}(m,t)$ evolves as a function of time and subsystem size $|A|$. Varying the latter provides information about how well the system is ordered at length scale $|A|$.

We find that when the initially ordered system is quenched well into the unbroken symmetry phase of the Hamiltonian, the (local) order quickly disappears and the PDF acquires a simple Gaussian shape. In contrast, when the quantum quench is to an energy density where Hamiltonian eigenstates retain short-range order, the PDF exhibits a complex structure both at finite times and in the steationay state. In the following we focus on the example of the spin-$1/2$ Heisenberg XXZ chain, but note that the picture we put forward is general and has a wide range of applicability.

Model and Setup. —

We investigate the time evolution of antiferromagnetic (short-ranged) order after a quantum quench in the spin-$1/2$ XXZ chain

[TABLE]

Here $S^{\alpha}_{j}$ are spin-$1/2$ operators acting on the site $j$ and we restrict our analysis to the range $\Delta>0$. The phase diagram of (2) is well established: at $T=0$ there is a BKT phase transition at $\Delta=1$ that separates a quantum critical phase at $\Delta<1$ and an antiferromagnetically ordered phase at $\Delta>1$. At any finite temperature the antiferromagnetic order melts. The Hamiltonian (2) is invariant under rotations by an arbitrary angle around the z-axis, translations by one site, and rotations around the x-axis by 180 degrees. In the thermodynamic limit at $\Delta>1$ and zero temperature the last symmetry gets broken spontaneously and one of the two degenerate ground states $|{\rm GS}^{\pm}_{\Delta}\rangle$, characterized by equal but opposite expectation values of the staggered magnetization per site, gets selected. In the Ising limit $\Delta\to\infty$ the ground states become the classical Néel states, i.e. $|{\rm GS}^{+}_{\infty}\rangle=|\cdots\uparrow\downarrow\uparrow\downarrow\cdots\rangle$ and $|{\rm GS}^{-}_{\infty}\rangle=|\cdots\downarrow\uparrow\downarrow\uparrow\cdots\rangle$. In order to investigate the melting of antiferromagnetic order we consider the following quantum quench protocol: (i) we prepare the system in the classical Néel state $|\Phi_{0}\rangle=|{\rm GS}^{+}_{\infty}\rangle$, which exhibits saturated antiferromagnetic long-ranged order; (ii) We consider unitary time evolution with Hamiltonian $H_{\Delta}$. The state of the system at time $t$ is thus $|\Phi_{t}\rangle=\exp(-iH_{\Delta}t)|\Phi_{0}\rangle$. This quench is integrable Pozsgay:13a ; FE_13b ; FCEC_14 and exact results on the stationary state are available neel_overlap ; QANeel ; XXZung ; XXZunglong ; Ilievski2015 ; IQNB . We employ the infinite Time-Evolving Block-Decimation (iTEBD) algorithm iTEBD1 ; iTEBD2 to obtain a very accurate description of $|\Phi_{t}\rangle$ in the thermodynamic limit. However, the growth of the bipartite entanglement entropy limits the time window accessible by this method. Retaining up to $\chi_{max}=1024$ auxiliary states, we are able to reach a time $t_{max}\simeq 12$ without significant error ($\lesssim 10^{-3}$).

PDF dynamics. —

Detailed information on how the antiferromagnetic order melts as the system evolves in time is provided by the probability distribution of the staggered subsystem magnetisation $M_{\ell}\equiv\sum_{j=1}^{\ell}(-1)^{j}S^{z}_{j}$

[TABLE]

where second line follows from the fact that the eigenvalues of $M_{\ell}$ are half-integer numbers. We note that the probabilities satisfy the normalisation condition $\sum_{m=-\ell/2}^{\ell/2}\widetilde{P}_{\ell}(m,t)=1$. The initial Néel state is an eigenstate of the staggered subsystem magnetization $M_{\ell}$ and concomitantly the probability distribution is a delta function $P_{\ell}(m,0)=\delta(m-\ell/2)$. This reflects the long-range magnetic order in the initial state. In Fig. 1 we show the evolution of $P_{\ell}(m,t)$ in time obtained by iTEBD for subsystem size $\ell=15$ and several values of the interaction strength $\Delta$ in the “post-quench” Hamiltonian (2).

We observe that the probability distribution depends strongly on $\Delta$: for small values of $\Delta$ the antiferromagnetic short-ranged order melts quickly and $P_{\ell}(m,t)$ is narrowly peaked around its average, which exhibits a damped oscillatory behaviour around zero BPGDA09 . The behaviour for $\Delta\gtrsim 2$ is very different: short-ranged order persists for some time while the probability distribution broadens and becomes more symmetric in $m$. This nicely chimes with the expectation (see below) that in the stationary state reached at late times the probability distribution to become symmetric in $m$. In Fig. 2 we plot the weights of $P_{\ell}(m,t)$ at several times and compare them to a Gaussian approximation based on the first two moments $\bar{m}(t)=\langle\Phi_{t}|M_{\ell}|\Psi_{t}\rangle$, $\sigma^{2}(t)=\langle\Phi_{t}|M_{\ell}^{2}|\Psi_{t}\rangle-\bar{m}^{2}(t)$

[TABLE]

We see that at $\Delta=0$ the probability distribution is approximately Gaussian at all times, while for $\Delta=2,3$ it exhibits a pronounced even/odd structure at short times and even at the latest times shown is strongly non-Gaussian.

“Small”- $\Delta$ regime. —

At small values of $\Delta$ and short and intermediate times we can use a time-dependent self-consistent mean-field approximation to determine the evolution of $P_{\ell}(m,t)$. We first map the Hamiltonian (2) to a model of spinless fermions by means of a Jordan-Wigner transformation, where we use the positive (negative) z-direction in spin space as quantization axis for even (odd) sites (see Supplementary material). This results in a spinless fermion Hamiltonian

[TABLE]

where $n_{j}=c^{\dagger}_{j}c_{j}$ and $\{c_{j},c^{\dagger}_{k}\}=\delta_{j,k}$. The staggered subsystem magnetization maps to $M_{\ell}=\sum_{j=1}^{\ell}(1/2-n_{j})$, while the initial Néel state maps to the fermion vacuum $|\Psi_{0}\rangle=|0\rangle$. Our self-consistent approximation corresponds to the replacement

[TABLE]

which leads to an explicitly time-dependent Hamiltonian $H_{\rm MF}(t)$, cf. Ref. SotiriadisCardy, . The expectation values in (6) are calculated self-consistently $\langle.\rangle_{t}=\langle\Psi_{t}|.|\Psi_{t}\rangle$, where

[TABLE]

Following Ref. Groha, we can express the characteristic function of $P_{\ell}(m,t)$ as a determinant of a $2\ell\times 2\ell$ matrix (see Supplementary material), which is easily evaluated numerically. This provides us with exact results at $\Delta=0$ for all times, cf. Fig. 3, and a highly accurate short-time approximation even for $\Delta=3$ as is shown in Fig. 4.

Late times and large $\Delta$. —

We now turn to the behaviour at late times after the quench. The stationary state is characterized by a finite correlation length $\xi(\Delta)$. On length scales $\ell\lesssim\xi(\Delta)$ we expect short-ranged antiferromagnetic order to remain, while it will have melted at scales $\ell>\xi(\Delta)$. We also expect the spin-rotational symmetry by 180 degrees around the x-axis to be restored in the stationary state as we are dealing with a one dimensional system with short-range interactions. The situation is completely analogous to that at finite temperatures – in fact adding a very small integrability-breaking term to the Hamiltonian would result in a steady state that is very close to the thermal state of the XXZ chain QANeel . In contrast to the steady state after our quench, the probability distribution of the staggered subsystem magentization at finite temperature $P_{\ell}(m,\beta)$ can be computed by matrix product state methods and for the aforementioned reasons it is instructive to consider it. Results for two values of $\Delta$ are shown in Fig. 5. We see that the probability distributions are symmetric in $m$, reflecting the unbroken symmetry of rotations around the x-axis by 180 degrees. At $\Delta=4$ we further observe that when the subsystem size exceeds the thermal correlation length $\xi_{\Delta}(\beta)$ antiferromagnetic short-ranged order melts and we obtain a Gaussian probability distribution centred around $m=0$. On the other hand, for $\ell\lesssim\xi_{\Delta}(\beta)$ the probability distribution is very broad and peaked at the maximal values $\pm\ell/2$, signalling the presence of both kinds of antiferromagnetic short-ranged order. For $\Delta=1$ the thermal correlation length is smaller than one lattice site in the temperature regime shown, which is why no traces of short-range order are visible and the probability distribution is a Gaussian centred around $m=0$.

The large-$\Delta$ regime is characterized by a low density of excitations and it is therefore possible to understand the behaviour observed above by combining a $1/\Delta$-expansion with a linked-cluster expansion along the lines of Refs EK:finiteT, ; JEK, ; GKE, ; PT10, ; ST12, ; CEF1, ; CEF2, ; SE:Ising, ; sine-Gordon, . As the physics we wish to describe is not tied to integrability, and the non-integrable case is easier to discuss, we focus on the latter. We consider the regime $\Delta\gg 1$ and break integrability by adding a small perturbation to the Heisenberg Hamiltonian, e.g. consider time evolution under $H=H_{\Delta}+\Delta^{-n}V$, where $n$ is a positive integer and $V$ some perturbation involving short-ranged spin-spin interactions that has the same symmetries as $H_{\Delta}$. We define linked clusters following the general formalism of Ref. EK:finiteT, and then implement a $1/\Delta$-expansion through a unitary transformation $\widetilde{H}=e^{iS}He^{-iS}$ SchriefferWolf , see the Supplementary Material for details. The result is an expansion of the stationary state density matrix of the form

[TABLE]

where $\rho_{\rm SS}^{(j)}$ are given as power series in $1/\Delta$. The leading term in the expansion is $\rho_{\rm SS}^{(0)}=\frac{1}{2}\sum_{\sigma=\pm}|{\rm GS}_{\Delta}^{\sigma}\rangle\langle{\rm GS}_{\Delta}^{\sigma}|$, where $|{\rm GS}_{\Delta}^{\sigma}\rangle$ are the two ground states of the model at anisotropy $\Delta$. The small parameter $e^{-\beta_{\rm eff}\Delta/2}$ is proportional to the density of domain-wall excitations over the ground states at large $\Delta$. The expansion (8) of the steady-state density matrix leads to a corresponding expansion of the probability distribution of the staggered subsystem magnetization $P_{\ell}(m,\infty)=\sum_{j}P^{(j)}_{\ell}(m)$ (see Supplemental Material)

[TABLE]

[TABLE]

where the dots denote subleading terms in $1/\Delta$. The expansions (10) hold as long as the subsystem size $2\ell$ is small compared to the correlation length in $\rho_{\rm SS}$ and establish that for large anisotropies $\Delta$ the probability distribution in the steady state is symmetric in $m$ and close to the average over the two ground states. In addition there is an exponentially suppressed ”background” contribution arising from a dilute gas of domain walls.

“Symmetrization” of the PDF in time. —

A characteristic feature of the time-evolution of $P_{\ell}(m,t)$ is that it becomes increasingly symmetric in $m$. In order to ascertain the associated time scale in the most interesting large-$\Delta$ regime it is useful to compare the probabilities for $M_{\ell}$ to be maximal ($\ell/2$) or minimal ($-\ell/2$) respectively. Results for $\Delta=6$ are shown in Fig. 6.

We see that $P_{\ell}(-\ell/2,t)$ grows linearly in time, while $P_{\ell}(\ell/2,t)$ shows a corresponding linear decrease. For the integrable XXZ chain the associated velocity is expected to be the maximal group velocity of elementary excitations over the stationary state BEL . In presence of weak integrability-breaking interactions in the large-$\Delta$ regime we expect qualitatively similar prethermal behaviour prethermal .

Conclusions. —

We have considered the full quantum mechanical probability distribution of the staggered subsystem magnetization $P_{\ell}(m,t)$ after a quantum quench from a classical Néel state to the spin-1/2 Heisenberg XXZ chain. We have shown that $P_{\ell}(m,t)$ provides detailed information how short-range antiferromagnetic order melts and have shown how to understand our numerical findings by analytical approaches valid in certain limits. Our setup should be realizable in cold-atom experiments like the ones by the Harvard group Greiner .

Acknowledgments. —

We are grateful to P. Calabrese, A. De Luca and M. Fagotti for stimulating discussions and to the Erwin Schrödinger International Institute for Mathematics and Physics for hospitality and support during the programme on Quantum Paths. M.C. thanks the Galileo Galilei Institute in Florence for hospitality during the workshop “Entanglement in quantum systems”. This work was supported by the the European Union’s Horizon 2020 under the Marie Sklodowska-Curie grant agreement No. 701221 NET4IQ (M.C and F.H.L.E.), the EPSRC under grant EP/N01930X (F.H.L.E.), the National Science Foundation under Grant No. NSF PHY-1748958 (F.H.L.E.) and by the BMBF and EU-Quantera via QTFLAG (M.C.).

Appendix A SUPPLEMENTARY MATERIAL

A.1 Self-consistent time dependent mean-field approximation

In the small-$\Delta$ regime the quench dynamics of the Heisenberg XXZ chain can be analyzed by means of a self-consistent time-dependent mean-field theory SotiriadisCardy as we will now demonstrate. For convenience we rotate the spin-quantization axis on all odd sites, which maps the classical Néel state on the saturated ferromagnetic state with all spins up. This transformation is induced by $P=\prod_{j\;{\rm odd}}\sigma^{x}_{j}$ and we have

[TABLE]

while the Hamiltonian maps onto

[TABLE]

The probability distribution of the staggered subsystem magnetization in the XXZ chain is then expressed as

[TABLE]

Performing the Jordan-Wigner transformation to spinless fermions $\{c_{i},c^{{\dagger}}_{j}\}=\delta_{ij}$

[TABLE]

the Hamiltonian becomes (up an unimportant constant)

[TABLE]

where $n_{i}=c^{{\dagger}}_{i}c_{i}$. The initial state maps onto the fermion vacuum, i.e.

[TABLE]

We now decouple the four-fermion interaction in a self-consistent fashion

[TABLE]

which leads to a time-dependent mean-field Hamiltonian

[TABLE]

Here we have defined

[TABLE]

The expectation values at time $t$ in (LABEL:SCMF) and (9) are calculated self-consistently

[TABLE]

The self-consistent mean-field approximation obtained in this way is expected to work well for small values of $\Delta$ and sufficiently short times. As we will see it works quite well even for intermediate values of $\Delta$.

The Heisenberg equations of motion in our self-consistent mean-field approximation are

[TABLE]

It is convenient to cast them in matrix form

[TABLE]

where ${\bf d}^{{\dagger}}\equiv(c^{{\dagger}}_{1},\dots,c^{{\dagger}}_{L},c_{1},\dots,c_{L})$ and we have defined a $2L\times 2L$ correlation matrix $\mathbb{C}(t)$ by

[TABLE]

The time evolution of the correlation matrix is governed by a first-order differential equation

[TABLE]

Using that our initial state is the fermion vacuum we have

[TABLE]

The correlation matrix at time $t$ can now be straightforwardly obtained by numerically integrating (14). By comparing the results to iTEBD computations we observe very good agreement for small values of $\Delta$. For large values $\Delta>1$ the approximation still gives a fair description of the dynamics at short times $t\lesssim 1$.

A.1.1 Full counting statistics

A nice feature of the self-consistent mean-field approximation described above is that it allows us to determine the probability distribution $P_{\ell}(m,t)$. To that end we define the associated characteristic function by

[TABLE]

Following Ref. Groha, we can derive a determinant representation for $F_{\ell}(\theta,t)$ that can be efficiently evaluated. We start by introducing Majorana fermion operators by

[TABLE]

These have anticommutation relations $\{a^{\alpha}_{i},a^{\beta}_{j}\}=2\delta_{\alpha\beta}\delta_{ij}$ and the observable of interest is bilinear in them

[TABLE]

Next we introduce an auxiliary reduced density matrix acting on sites $1$,…,$\ell$ by

[TABLE]

where ${\bf a}^{{\dagger}}=(a^{x}_{1},\dots,a^{x}_{\ell},a^{y}_{1},\dots,a^{y}_{\ell})$, $\mathbb{W}=i\sigma^{y}\otimes\mathbb{1}_{\ell}$ and

[TABLE]

The characteristic function can then be expressed as

[TABLE]

where $\rho_{\ell}(t)$ is the reduced density matrix of the subsystem $[1,\ell]$ at time $t$

[TABLE]

Finally we use an identity for the trace of the product of two Gaussian density matrices proved in Ref. Fagotti,

[TABLE]

where $\Gamma_{j}$ are correlation matrices defined by

[TABLE]

In our case the two correlation matrices are

[TABLE]

where $\mathbb{C}(t)$ is the correlation matrix obtained from (14) and

[TABLE]

Our final result for the characteristic function of the probability distribution $P_{\ell}(m,t)$ is thus

[TABLE]

Combining (27) and (14) we can obtain numerically exact results in the framework of our time-dependent self-consistent mean-field approximation for $F_{\ell}(\theta,t)$ in the thermodynamic limit and large subsystem sizes.

Appendix B Combined linked-cluster and $1/\Delta$ expansions

In the large-$\Delta$ regime we can gain insights about the late time behaviour after the quench by means of combined linked-cluster and $1/\Delta$ expansions. We again carry out a rotation of the spin-quantization axis, which leads us to a Hamiltonian of the form

[TABLE]

Here $\tau^{\alpha}_{j}$ are Pauli matrices and we have introduced an energy scale $J$, which should be set to $1/4$ in order to arrive at the conventions used in the main text. In the rotated basis the staggered subsystem magnetization maps onto

[TABLE]

and our goal is to calculate

[TABLE]

where $|\uparrow\rangle$ is the saturated ferromagnetic state. As the physics we are interested in does not rely on the integrability of the spin-1/2 XXZ chain and the non-integrable case is easier to analyze, we proceed by adding a short-ranged integrability-breaking interaction $V$ that does not break any of the discrete symmetries of $H^{\prime}(\Delta)$ and consider

[TABLE]

where $n$ is a positive integer. In the regime $\Delta\gg 1$ we may carry out a $1/\Delta$-expansion following Ref. SchriefferWolf, . We consider a basis transformation of the form

[TABLE]

where $S^{(j)}$ are to be chosen in such a way that order by order in the expansion $\widetilde{H}$ commutes with the domain-wall number operator

[TABLE]

The first order term in (32) is

[TABLE]

In the transformed basis the Hamiltonian has the following expansion

[TABLE]

and domain wall number is a good quantum number. The ground states of $\widetilde{H}$ are the saturated ferromagnetic states $|\uparrow\rangle=|\uparrow\uparrow\dots\rangle$ and $|\downarrow\rangle=|\downarrow\downarrow\dots\rangle$ with energy $E_{0}=-J\Delta L$ (in the original spin basis these correspond to the two Néel states). Low-lying excited states involve ferromagnetic domain walls. In order to analyze properties at finite but low energy densities by means of the method first introduced in Ref. EK:finiteT, we need to impose open boundary conditions on (LABEL:Htilde) so that excited states with an odd number of domain walls are allowed. To deal with this complication we need to shift the subsystem for which we determine the probability distribution to the centre of the open chain, i.e. for even $\ell$ we take

[TABLE]

As we are interested in $\ell$ fixed and $L\to\infty$ boundary terms will not contribute and we will recover a translationally invariant result. Excited states involving a single domain wall can be expanded in a basis formed by the states

[TABLE]

and their spin-reversed analogues $|\bar{n}\rangle=C|n\rangle$ where $C=\prod_{j=1}^{L}\tau^{x}_{j}$. A basis of eigenstates of (LABEL:Htilde) in the sector with a single domain wall is given by $\{|Q^{\pm}_{\alpha}\rangle,C|Q^{\pm}_{\alpha}\rangle\}$, where

[TABLE]

and

[TABLE]

The corresponding excitation energies are

[TABLE]

Energy eigenstates in the original spin basis are obtained by undoing the basis transformation (32). In particular we have (up to boundary terms)

[TABLE]

The density matrix describing the steady state after a quantum quench from the classical Néel state is defined by the requirement

[TABLE]

where ${\cal O}$ is any local operator. In our case $\rho_{\rm SS}$ is a thermal density matrix (as we have broken integrability though the $V$-term) and in the large-$\Delta$ regime we have (in a large, finite volume) EK:finiteT ; JEK ; GKE ; PT10 ; ST12 ; CEF1 ; CEF2 ; SE:Ising ; sine-Gordon

[TABLE]

where

[TABLE]

The first two terms of the low-density expansion involve only the ground states and single domain wall excitations

[TABLE]

By construction of the linked cluster expansion the last term precisely subtracts the contribution that diverges in the infinite volume limit. Using the explicit form (B) of the $1/\Delta$ expansion for the ground states we obtain an explicit expression for the leading term in the characteristic function

[TABLE]

where we have used that the second order correction to the ground state has zero overlap with the ferromagnetic states. We see that in this contribution the probabilities of $S^{z}_{\ell}$ taking values less than the maximal possible ones $\pm\ell/2$ are suppressed by powers of $1/\Delta$. The subleading term in the low-density (linked cluster) expansion is obtained from the matrix elements

[TABLE]

Evaluating the matrix elements gives

[TABLE]

where for even $\ell$ we have

[TABLE]

Combining (48) with (45) and turning the momentum sums into integrals we obtain the following explicit expressions for the first subleading contribution to the characteristic function

[TABLE]

In Fig. 7 we compare the PDF in the thermal state with inverse temperature fixed by the condition $\langle{\rm GS}_{\infty}^{+}|H_{\Delta}|{\rm GS}_{\infty}^{+}\rangle={\rm Tr}\left[\rho_{\rm SS}H_{\Delta}\right]$ with the results of the large-$\Delta$/low density expansion. As expected, the agreement is increasingly better for larger values of the anisotropy.