Spin fluctuations after quantum quenches in the S=1 Haldane chain: numerical validation of the semi-semiclassical theory

TL;DR

This paper validates a semi-semiclassical method for describing spin fluctuations after quantum quenches in the S=1 Haldane chain, showing it agrees well with TEBD simulations and improves upon standard semiclassical approaches.

Contribution

The paper introduces and validates a hybrid semi-semiclassical approach that accurately captures quantum quench dynamics in the S=1 Haldane chain, surpassing traditional semiclassical methods.

Findings

Semi-semiclassical method matches TEBD results.

Standard semiclassical approach fails to describe dynamics.

Full scattering matrix improves predictive accuracy.

Abstract

We study quantum quenches in the $S=1$ Heisenberg spin chain and show that the dynamics can be described by the recently developed semi-semiclassical method based on particles propagating along classical trajectories but scattering quantum mechanically. We analyze the non-equilibrium time evolution of the distribution of the total spin in half of the system and compare the predictions of the semi-semiclassical theory with those of a non-Abelian time evolving block decimation (TEBD) algorithm which exploits the SU(2) symmetry. We show that while the standard semiclassical approach using the universal low energy scattering matrix cannot describe the dynamics, the hybrid semiclassical method based on the full scattering matrix gives excellent agreement with the first principles TEBD simulation.

Tables

Table 1. Table 1: Left: Ground state spin structure of a semi-infinite spin-1 Haldane chain. Right: Slopes of the probability distributions at short times, r S = d ​ P ​ ( S , t ) / d ​ t | t = 0 subscript 𝑟 𝑆 evaluated-at d 𝑃 𝑆 𝑡 d 𝑡 𝑡 0 r_{S}={\rm d}P(S,t)/{\rm d}t|_{t=0} , as predicted by the semiclassical theory, and as extracted from the TEBD simulations for quench duration t Q = 1.6 / J 1 subscript 𝑡 𝑄 1.6 subscript 𝐽 1 t_{Q}=1.6/J_{1} .

aa $S_0$ aa	a $P_0\left(S_0\right)$ a
0	0.2426
1	0.7388
2	0.0186

Equations

$$\hat{H}(t)=\sum _i{J}_1{\vec{S}}_i\cdot {\vec{S}}_{i+1}+J_2(t){\left({\vec{S}}_i\cdot {\vec{S}}_{i+1}\right)}^2.$$

$$i{\partial }_t|\Psi (t)⟩=\hat{H}(t)|\Psi (t)⟩$$

$$|\Psi ⟩=\sum _{\{{\sigma }_l\}}(\dots M_{{\sigma }_{l-1}}^{[l-1]}{M}_{{\sigma }_l}^{[l]}\dots )|\{{\sigma }_l\}⟩,$$

$$|\Psi ⟩=\sum _{t_l}\frac{{\Lambda }_{t_l}^{[l]}}{\sqrt{2J_{t_l}+1}}\sum _{m_l=-J_{t_l}}^{J_{t_l}}{|t_l,m_l⟩}_{\mathrm{left}}{|t_l,{\overline{m}}_l⟩}_{\mathrm{right}},$$

$$\begin{gathered} {|t_l,m_l⟩}_{\mathrm{left}}=\sum _{t_{l-1},{\eta }_l}{A}_{t_{l-1}{\eta }_l}^{[l]t_l}\sum _{m_{l-1},{\sigma }_l}{C}_{m_{l-1}{\sigma }_l}^{m_l{\eta }_l} \\ {|t_{l-1},m_{l-1}⟩}_{\mathrm{left}}|{\sigma }_l⟩. \end{gathered}$$

$$(M_{{\sigma }_l}^{[l]}{)}_{t_{l-1},m_{l-1}}^{t_l,m_l}=\sum _{{\eta }_l}(A^{[l]}{)}_{t_{l-1}{\eta }_l}^{t_l}{C}_{m_{l-1}{\sigma }_l}^{m_l{\eta }_l},$$

$$P(S,t)=\sum _{t_l|J_{t_l}=S}|\Lambda (t{)}_{t_l}{|}^2.$$

$${\omega }_q=\sqrt{{\Delta }^2+c^2{q}^2},$$

$$\Delta \approx 0.4105J_1,c\approx 2.472J_{1}a.$$

$${\rho }^{-1}>c/\Delta .$$

$$P(S,t)=\sum _{S_0}\sum _{S_{\mathrm{qp}}}{\delta }_{S\in S_0\otimes S_{\mathrm{qp}}}\left(2S+1\right)\frac{P_0(S_0)}{2S_0+1}\frac{P_{\mathrm{qp}}(S_{\mathrm{qp}};t)}{2S_{\mathrm{qp}}+1}.$$

$$Q(t)=\rho {\int }_0^{c}vtf(v)dv\equiv \frac{t}{2\tau },$$

$$P_{\text{qp}}(S_{\text{qp}}=1;t)$$

$$P_{\text{qp}}(S_{\text{qp}}=0;t)$$

$$P(0,t)=P_0(0)-\left(P_0(0)-\frac{P_0(1)}{9}\right)\frac{t}{2\tau }+\dots ,$$

$$P(S,t)=P_0(S)$$

$$+\frac{2S+1}{6}\left(\frac{P_0(S-1)}{2S-1}-\frac{2P_0(S)}{2S+1}+\frac{P_0(S+1)}{2S+3}\right)\frac{t}{\tau }+\dots$$

$$p(n_{+},n_{-})=\frac{1}{n_{+}!n_{-}!}q(t{)}^{n_{+}+n_{-}}{e}^{-2q(t)},$$

$$\begin{gathered} P_{\text{qp}}^{\text{refl}}(S_{\text{qp}},t)=\sum _{n_{+},n_{-}=0}^{\infty }p(n_{+},n_{-})\frac{1+(-1{)}^{S_{\text{qp}}}(-1{)}^{n_{+}+n_{-}}}{2} \\ =\frac{1}{2}(1+(-1{)}^{S_{\text{qp}}}{e}^{-t/\tau }). \end{gathered}$$

$$P_{\mathrm{qp}}^{\mathrm{trans}}(S_{\mathrm{qp}};t)=\sum _N{P}_{\mathrm{spin}}(S_{\mathrm{qp}}|N)P_{\mathrm{cross}}(N;t),$$

$$M_{S|N}=M_{S+1|N-1}+(1-{\delta }_{S,0})\left(M_{S|N-1}+M_{S-1|N-1}\right)$$

$$n(q)\propto \frac{q^2}{{\sinh }^2({\omega }_q/2E_0)},$$

$${\hat{H}}_0=\sum _{q,\sigma }{\omega }_q{b}_{q,\sigma }^{†}{b}_{q,\sigma }.$$

$$\begin{gathered} \hat{H}=\sum _{q,\sigma }\{{\omega }_q(t)b_{q,\sigma }^{†}{b}_{q,\sigma }+ \\ \frac{1}{2}{g}_q(t)(-1{)}^{\sigma }\left(b_{q,\sigma }^{†}{b}_{-q,-\sigma }^{†}+b_{q,\sigma }{b}_{-q,-\sigma }\right)\}, \end{gathered}$$

$${\omega }_q(t)={\omega }_q+J_2\lambda (t){\Delta }_q,g_q(t)=J_2\lambda (t)g_q.$$

$$n(q)\propto|g_{q}|^{2}\;\frac{J_{2}^{2}}{\omega_{q}^{2}}\;\;\biggl{|}\widetilde{\frac{{\rm d}\lambda}{{\rm d}t}}(2\omega_{q})\biggr{|}^{2}\,,$$

$$\lambda (t)=\frac{1}{2}\left(1-\tanh (t/t_Q)\right).$$

$$n(q)\propto J_2^2\frac{q^2}{{\sinh }^2(\pi t_Q{\omega }_q)}.$$

$$\rho \lesssim {\rho }_e/\Delta .$$

$$c{\rho }_e<{\Delta }^2$$

$${\left[U_{\mathrm{sym}}\right]}_{{\eta }_l^{\prime }{\eta }_{l+1}^{\prime }}^{{\eta }_l{\eta }_{l+1}}=\sum _{\begin{array}{c}{m}_{l-1}\\ m_l\\ m_l^{\prime }\\ m_{l+1}\end{array}}\sum _{\begin{array}{c}{\sigma }_l\\ {\sigma }_{l+1}\\ {\sigma }_l^{\prime }\\ {\sigma }_{l+1}^{\prime }\end{array}}\frac{1}{2J_{l+1}+1}{U}_{{\sigma }_l^{\prime }{\sigma }_{l+1}^{\prime }}^{{\sigma }_l{\sigma }_{l+1}}$$

$$C_{m_{l-1}{\sigma }_l}^{m_l{\eta }_l}{C}_{m_l{\sigma }_{l+1}}^{m_{l+1}{\eta }_{l+1}}{\left(C_{m_{l-1}{\sigma }_l^{\prime }}^{m_l^{\prime }{\eta }_l^{\prime }}\right)}^{\ast }{\left(C_{m_l^{\prime }{\sigma }_{l+1}^{\prime }}^{m_{l+1}{\eta }_{l+1}^{\prime }}\right)}^{\ast }.$$

Full text

I Introduction

Understanding non-equilibrium dynamics in interacting quantum many-body systems is one of the major challenges in today’s statistical physics Cazalilla2010 ; Polkovnikov2011 ; Eisert2015 ; Calabrese2016 . The formation of current-carrying steady states Prosen2015 ; Bertini2016 ; Bernard2016a ; Collura2018 , the details of the thermalization process Moeckel2008 ; Bertini2015a , entropy production Calabrese2005 ; Schuch2008 ; Alba2017 ; Nahum2017 ; Collura2018a or the interplay with disorder Znidaric2008 ; Bardarson2012 ; Nanduri2014 ; Serbyn2014 ; Vasseur2016 and topology Halasz2012 ; Mazza2014 ; DAlessio2015 ; McGinley2018 are just a few examples of the intriguing open questions which, due to recent breakthroughs in quantum simulation, can now be investigated experimentally.

While experimental results are abounding Kinoshita2006 ; Trotzky2012 ; Gring2012 ; Chenau2012 ; Langen2015 ; Kaufman2016 , available theoretical tools are quite limited in number and power. One dimensional systems represent in this regard an exception and a theoretical testing ground for all methods and investigations, especially since in one dimension powerful analytical and numerical methods such as Bethe AnsatzKorepin_book , bosonization Haldane_Bosonization ; Giamarchi_book , Density Matrix Renormalization Group (DMRG) White1992 ; Schollwock2011 exist to study equilibrium properties in detail. These methods can be extended to non-equilibrium steady statesProsen2009 ; Cui2015 , while dynamics under non-equilibrium conditions can be efficiently simulated by the tDMRG White2004 and the Time-Evolving Block Decimation (TEBD) algorithms Vidal2007 . However, in dynamical simulations “exactness” or precision is typically lost after relatively short times due to the rapidly increasing entanglement of the state.

Recently, a semi-semiclassical approach (SSA) has been proposed to study generic, gapped one dimensional quantum systems at longer times, and demonstrated on the sine–Gordon model PascuMarciGergo2017 . As a generalization of the original semiclassical approach of Sachdev, Young, and Damle SachdevYoung1997 ; SachdevDamle1997 , the semi-semiclassical method treats trajectories of quasiparticles classically, but it accounts for the precise quantum evolution of the internal degrees of freedom fully quantum-mechanically, and allows one to capture the associated quantum entanglement generation as well as simultaneous probabilistic processes. Though approximate, the SSA method is versatile, conceptually simple, and has also been extended to study dynamics in non-equilibrium steady states (NESS) within the nonlinear sigma model Marci_Pascu_Gergo_NESS .

Here we investigate quasiparticle creation and spin propagation after a quantum quench within the Haldane-gapped phase of the the antiferromagnetic spin-1 Heisenberg chain and compare the predictions of the semi-semiclassical method with TEBD simulations in detail. In particular, we consider the $S=1$ Haldane chain with a time-dependent biquadratic interaction $J_{2}(t)$

[TABLE]

Here $J_{2}(t)=J_{2}\gamma(t)$, and the function $\gamma(t)$ changes from $1$ to [math] around $t=0$. We start from the ground state of $H(-\infty)$ and perform infinite volume TEBD (iTEBD) simulations Vidal2007 to obtain the full time dependent wave function of the chain, $|\psi(t)\rangle$. To reach sufficiently long times, comparable with the quasiparticle collision times, we need to exploit non-Abelian symmetries in our simulations.

The quench protocol described above generates a gas of (entangled pairs of) quasiparticles in the final state. Cutting the infinite chain into two at time $t$ and then measuring the spin distribution on one side, $P(S,t)$, we can explore the propagation and collision of these quasiparticles (Fig. 1). We can, in particular, compute $P(S,t)$ in terms of the semi-semiclassical approach, and compare it to the results of TEBD simulations to find an astonishing agreement. As we shall also demonstrate, a theory based on completely reflective quasiparticle collisions SachdevDamle1997 ; Damle2005 ; RappZarand ; Evangelisti2013 ; Kormos_Zarand_PRE_2016 is not able to account for the observed behavior, and the full semi-semiclassical approach is needed to get agreement with TEBD computations.

The paper is structured as follows. In Sec. II we overview the basics of the microscopic non-Abelian iTEBD simulations and the SSA method. In Sec. III the short time ballistic behavior after the quench is analyzed. Then in Sec. IV we discuss the collision dominated regime of the dynamics. In Sec. V we develop a perturbative quench theory and test its scope of validity. Our conclusions are summarized in Sec. VI.

II Basic concepts and numerical methods

Before presenting the main results, let us shortly review the two methods used to investigate the quench dynamics.

II.1 Non-abelian TEBD lattice simulations

Here we discuss our non-abelian TEBD algorithm only briefly, since our flexible and general way of treating non-Abelian symmetries in MPS simulations will be presented in a separate publication Werner_future .

The post-quench dynamics of the system is described by the microscopic many-body Schrödinger equation,

[TABLE]

which we solve numerically by means of the infinite chain time-evolving block decimation (iTEBD) algorithm Vidal2007 . The TEBD algorithm describes the real time dynamics of the system based on the matrix product state (MPS) description of the quantum state $\left|\Psi(t)\right\rangle$ Schollwock2011 ; Orus2014 ,

[TABLE]

where the $|\sigma_{l}\rangle$ with $\sigma_{l}=\{0,+,-\}$ refer to the three quantum states of spin $S_{l}$ at site $l$.

The MPS factorization of the quantum state $\left|\Psi\right\rangle$ relies on Schmidt decomposition Schmidt1907 . Let us now focus on the case where $\left|\Psi\right\rangle$ is a spin-singlet, and cut the chain into two halves between sites $l$ and $l+1$, i.e. treat the Hilbert space of the full chain as the product of its left and right halves. The Schmidt decomposition of $\left|\Psi\right\rangle$ then reads asMcCulloch ; Singh_Vidal2010

[TABLE]

where $\left|t_{l},m_{l}\right\rangle_{\mathrm{left}}$ and $\left|t_{l},\overline{m}_{l}\right\rangle_{\mathrm{right}}$ are the so-called Schmidt-pairs footnote_right_state . Here $t_{l}$ labels multiplets, while $m_{l}$ refers to internal states within this multiplet having a total spin $J_{t_{l}}$. The values $\Lambda^{[l]}_{t_{l}}/\sqrt{2J_{t_{l}}+1}>0$ are the so called Schmidt values and are independent of the internal label $m_{l}$.

The non-Abelian MPS (NA-MPS) representation of $|\Psi\rangle$ can be constructed by relating neighboring left Schmidt states while moving the cut position forward by one site,

[TABLE]

Here the tensor $C^{m_{l}\,\eta_{l}}_{m_{l-1}\,\sigma_{l}}$ contains the Clebsch–Gordan coefficients $\left\langle J_{l-1},m_{l-1};S_{l},\sigma_{l}\right.\left|J_{l},m_{l}\right\rangle$, while the superindex $\eta_{l}\to\{J_{l-1},S_{l},J_{l}\}$ runs over allowed values of the three spins footnote_on_outer_multiplicity . An iterative application of Eq. (5) yields the non-Abelian MPS (NA-MPS) representation (3) with

[TABLE]

corresponding to the two-layer NA-MPS structure sketched in Fig. 3. The tensors $A^{[l]}$ describe the transformation at the level of multiplets, and the contraction of the index $\eta_{l}$ between $A^{[l]}$ and $C$ ensures that the values of the three representation indices in the two tensors match.

The upper layer of the tensors $A^{[l]}$ contains all important physical information on the state $|\Psi\rangle$, and our TEBD time-evolver operators act only on this upper NA-MPS layer (see Appendix A). In other words, all expensive contractions within the Clebsch–Gordan layer are eliminated.

Here we use the infinite-chain TEBD algorithm that applies for translation invariant statesVidal2007 . In this case, both the tensors $A^{[l]}$ and the Schmidt values $\Lambda^{[l]}_{t_{l}}$ are independent of the site index $l$. Numerically, this one site translation invariance is, however, lowered to a two-site translation invariance due to the Suzuki–Trotter time evolution scheme Trotter ; Suzuki , i.e. the tensors on the even and odd sublattices are slightly different.

The TEBD simulation provides the time-dependent tensors $A(t)^{\,t_{l}}_{t_{l-1}\,\eta_{l}}$ and the Schmidt values $\Lambda(t)_{t_{l}}$ (for details see Appendix A). The distribution of the total spin of the half chain is then simply related to the Schmidt values,

[TABLE]

As shown in later sections, this distribution contains the essential information on the spin structure and spin entanglement of excited quasiparticles in the post-quench state.

II.2 The semi-semiclassical approach

The numerical resources required for TEBD grow exponentially in time, and make the microscopic simulations tractable only for short times. However, the single-site resolved knowledge of the quantum state is not necessary for answering many questions. The model (1) is known to be gapped Haldane1 ; AffleckWhite2008 , and its low-energy excitations are $S=1$ triplet quasiparticles AffleckWhite2008 described by the $O(3)$ nonlinear sigma model, an integrable relativistic field theory. Low energy quasiparticles therefore have a (close to) relativistic dispersion

[TABLE]

where the gap and the “speed of light” are given in terms of the microscopic parameters of the Hamiltonian asAffleckWhite2008

[TABLE]

The finite gap ensures that, in the case of small quenches, the post-quench state is a dilute gas of quasiparticles.

From the locality and the translation invariance of the Hamiltonian we can also conclude that shortly after the quench the state consists of spatially localized uncorrelated quasiparticle pairs. Moreover, since the quench protocol conserves SU(2), these quasiparticle pairs must form SU(2) singlets. These assumptions form the basis of the semi-semiclassical approximation PascuMarciGergo2017 , where the quasiparticles’ spatial degrees of freedom are treated classically in a Monte Carlo sampling of the possible world-line configurations, while their internal spin states are followed at the quantum level. As sketched in Fig. 2, the evolution of the quasiparticle spin states can be described by consecutive application of the two-particle S-matrix at collision events – if the gas is dilute enough.

The condition for the applicability of semiclassical approach is that the mean interparticle distance must be larger than the Compton wavelength of the particlesfootnote_debroglie

[TABLE]

As we shall see, this condition is satisfied even for relatively large quenches in our model.

This semi-semiclassical formalism allows us to compute the time dependent spin distribution of a half-chain, $P(S,t)$. Half-chain spin fluctuations have two sources in the semi-semiclassical approach: (1) cutting the vacuum state of an infinite Haldane chain into two gives rise to a non-trivial spin distribution, $P_{0}(S_{0}),$ and (2) quench generated quasiparticles carry spins (and entanglement) across the cut.

To determine the first contribution, we first constructed the post-quench ground state with $J_{2}=0$ using TEBD, cut it into two, and determined $P_{0}(S_{0})$ from the Schmidt values. As shown in Table 1, we clearly observe the presence of two topologically protected spin 1/2 end spins after the cut AKLT ; Kennedy1990 : these yield a triplet state with almost 74% probability and a singlet state with more than 24% likelihood. If the infinitely long half chains were fully separated from each other already before the cut, these probabilities would be just 75% and 25% demonstrating that the states of the topological end spins of the left half chain are independent of each other leading to a triplet-singlet degeneracy in the ground state of the half chain. However, if the exchange coupling between the two half chains is finite before the cut, there is a small probability that virtual quasiparticles generated by vacuum-fluctuations cross the cut and alter the total spin of the half chain, although the probability of these processes is less than about 2%.

Semi-semiclassics can be used to determine the second contribution, that of pairs of quasiparticles created by the $J_{2}$ quench. This gives rise to the quasiparticles’ spin distribution, $P_{\mathrm{qp}}(S_{\mathrm{qp}})$. Assuming that the spin orientation of the quasiparticles is independent of those of the vacuum fluctuations, we obtain the total half chain spin distribution,

[TABLE]

The distribution $P_{\mathrm{qp}}(S_{\mathrm{qp}})$ can, in general, be determined only numerically. However, as described in the next two sections, we have simple analytical expressions at very short times as well as in the limits of completely reflective and completely transmissive scatterings, valid for very cold and very hot gases of quasiparticles, respectively.

III Short time ballistic behavior

First we analyze the initial change of $P(S)$ shortly after the quench, where the quasiparticle picture gives simple predictions. As discussed before, the initial state is a superposition of states containing randomly localized spin singlet quasiparticle pairs with random velocities $\pm v$, and a quasiparticles density $\rho$. The distribution of the magnitude of the velocity, $f(v),$ depends on the details of the quench. The total spin of the left half-chain changes when the first $S=1$ quasiparticle crosses the position of the cut and carries spin from one half to the other. The probability of this happening within time $t$ is simply

[TABLE]

because the quasiparticle of velocity $v$ can come from an interval of length $|v|t$ on either side, touching the cut, but it must move to the right direction. The collision time $\tau$ above is defined as the ratio of the mean inter-particle spacing and the mean velocity.

For short times, we can neglect multiple crossings, so in Eq. (11)

[TABLE]

and we find

[TABLE]

Using the vacuum spin probabilities in Table 1, we can compute, independently of $f(v)$, the ratios of the slopes $r_{S}=\left.\frac{dP(S,t)}{dt}\right|_{t=0}$ of the initial linear time dependences with the result shown in Table 1.

We have confronted these predictions with microscopic TEBD simulations, shown in Fig. 4. The short time $P(S,t)$ functions are plotted for $S=0,1,2$ in panel (b) as functions of time, for a representative sudden quench. The relative rates $r_{S}/r_{0}$ extracted from these and similar curves obtained with various quench sizes are displayed in the main panel (a). In accordance with the quasiparticle picture, the relative rates are independent of the quench size. The agreement with the quasiparticle prediction is excellent for $S=0,1,2,$ while there is a small deviation for sudden quenches for $S=3.$ The latter can be attributed to the difficulty of extracting the universal initial rate due to transient oscillations (c.f. upper curve in panel (c), comparing the results for a sudden and a finite time quench for $S=3$). The oscillations are not present for smooth finite time quenches, and for these we get excellent agreement with the semiclassical prediction. This agreement between the prediction of the quasiparticle picture and the numerics gives strong evidence that the quasiparticle picture is correct. In the following sections we test the validity of the semiclassical description at longer times.

IV Collision dominated regime

At later times after the quench, several particles can cross from one half-chain to the other from both directions. Moreover, one has to take into account the effect of collisions. In this section, after considering two analytically tractable limiting cases, we apply the semi-semiclassical method and compare its results with the TEBD numerics.

IV.1 Simple limits

Let us first consider two limits, those of completely reflective and completely transmissive collisions, in which we can compute the spin distribution function analytically.

IV.1.1 Completely reflective limit

In the universal low-energy limit, the two-particle scattering matrix of gapped models with short range interaction is a permutation matrix, corresponding to perfect reflection of the incoming quasiparticles. This limiting S-matrix was used in the early works SachdevDamle1997 ; Damle2005 ; RappZarand ; Evangelisti2013 ; Kormos_Zarand_PRE_2016 on the semiclassical method to describe the dynamics at low temperatures. The S-matrix for the $O(3)$ nonlinear sigma model, describing our spin Hamiltonian, is exactly known, and also describes perfectly reflective processes in the limit of small relative rapidities (see Appendix B).

The initial state consists of pairs of quasiparticles that form spin singlets, and in this reflective limit, neighboring quasiparticle pairs remain singlets even after many collisions (see Fig. 5). If we cut the system into two half-chains, the quasiparticle contribution to the total spin of the left part at a given time is $S_{\text{qp}}=0$ if the cut lies between pairs, while $S_{\text{qp}}=1$ if the cut breaks a pair. It is easy to see that the first situation is realized if the number of quasiparticles crossing from one half to the other up to the given time is even, and the second if this number is odd. The crossing number can be computed using the straight lines in Fig. 5 (the would-be trajectories of non-interacting particles). The total number of crosses from the left and right, $n_{+}$ and $n_{-}$, follow independent Poisson distribution,

[TABLE]

where $q(t)=\frac{1}{2}\rho\,t\int_{0}^{c}dvvf(v)=Q(t)/2=1/(4\tau)$ is the probability that a trajectory crosses the cut from the left (or from the right). The probabilities of having even ($S_{\text{qp}}=0$) or odd ($S_{\text{qp}}=1$) number of crossings are then

[TABLE]

Note that expanding for short time we recover the expressions in Eqs. (13). Substituting the above result into Eq. (11) the total spin distribution can be computed in the totally reflective limit.

IV.1.2 Completely transmissive (ultrarelativistic) limit

Very high energy quasiparticles do not interact with each other. This can also be verified on the exact S-matrix of the nonlinear sigma model in Appendix B, in the limit of very large rapidity differences, i.e. of ultrarelativistic quasiparticles. In this limit, the original spin singlets remain singlets, but now the members of a pair keep moving away from each other following the light cone (see Fig. 6). Since quasiparticles crossing the cut are independent of each other, the quasiparticle spin distribution is given by

[TABLE]

where $P_{\mathrm{cross}}(N;t)=\frac{Q(t)^{N}}{N!}e^{-Q(t)}$ is the Poisson probability distribution of the number of worldlines that cross the cut from any side, while $P_{\mathrm{spin}}(S|N)=(2S+1)\,M_{S|N}/3^{N}$ is the distribution of the total spin of $N$ particles with random spin orientations. Here $M_{S|N}$ counts the multiplets of spin $S$ in the $N$-particle space, and it can be calculated iteratively from the recursion relation

[TABLE]

with initial condition $M_{S|0}=\delta_{S,0}$. Solving these equations iteratively, we can quickly compute $P_{\mathrm{qp}}^{\mathrm{trans}}$ at any time.

IV.2 Hybrid semiclassical dynamics

In general, the scattering is neither fully reflective nor fully transmissive. The result of each collision is instead a superposition of possible outgoing states with respective amplitudes given by the scattering matrix. In the O(3) nonlinear sigma model the total spin and also its $z$-component are conserved in the scattering of the $S=1$ quasiparticles, but transmissive and reflective processes as well as quasiparticle spin flips all occur with finite scattering amplitude. In the $S_{T}^{z}=0$ scattering channel, for example, a superposition of transmissions and reflections occur for incoming particles $(+,-)$, but even the process $(+,-)\longrightarrow(0,0)$ is allowed.

The two-body scattering matrix is exactly known for the O(3) nonlinear sigma model. In our hybrid semiclassical method we use this S-matrix: whenever there is a collision of quasiparticles, we act on the two colliding quasiparticle spins by the corresponding O(3) S-matrix. This goes beyond standard semiclassical treatments not only by allowing nontrivial scattering processes, but also by treating the spin part of the many-body wave function fully quantum mechanically PascuMarciGergo2017 .

Apart from the S-matrix, the other main input for the method is the momentum distribution of the quasiparticles, which is not easy to measure or calculate. In our simulation we used the distribution

[TABLE]

with $\omega_{q}=\sqrt{\Delta^{2}+q^{2}c^{2}}$ the quasiparticle dispersion relation and $E_{0}$ a tunable cutoff parameter. This particular functional form is motivated by the perturbative calculation outlined in Sec. V.

We compare the results of the semi-semiclassical method with the microscopic TEBD numerics for sudden quenches in Fig. 7. Different symbols correspond to quenches of different magnitude. In the semi-semiclassical approach, the distribution of scattering matrices depends exclusively on the velocity distribution of quasiparticles, which is, in turn, determined exclusively by the microscopic quench protocol. In the perturbative limit, we expect this distribution to be independent of the amplitude of the quench. The amplitude of the quench is only supposed to influences the density of the quasiparticles created, and thus the collision time. Therefore, we expect that the functions $P(S,t)$ can be scaled on the top of each other by rescaling time.

This is indeed what we find by studying sudden quenches of different sizes, $J_{2}/J_{1}\in[0.12,0.36]$, for which the curves $P(S,t)$ collapse when plotted against $t/\tau$ (see Fig. 7). For smaller quenches, the collision time $\tau$ is determined from the early time slope of the functions in Eq. (14), while for larger quenches, where transient oscillations are superposed at early times, we just rescaled the curves to achieve the best collapse. Remarkably, the same rescaling factor worked for all different spin values in each case.

The resulting curves can be compared to semiclassical predictions. The first observation is that neither of the special limits can reproduce the results of the microscopic numerical simulation (c.f. green dotted and red dash-dotted curves in panels (b)-(f)). In particular, the standard semiclassical method based on totally reflective collisions cannot account for the post-quench dynamics, and the numerically determined time evolution lies between the predictions of completely reflective and completely transmissive computations.

Remarkably, however, the semi-semiclassical results give an excellent agreement with the TEBD data. We recall that there is a single free parameter in the semiclassical calculation, the cutoff parameter $E_{0}$. The same choice, $E_{0}=4J_{1},$ gave the best agreement for all spin values $S$. This cutoff is close to the quasiparticles’ bandwidthAffleckWhite2008 .

The effect of the cutoff on the semi-semiclassical spin distributions for $S=1,\dots,4$ is illustrated in Fig. 8. A very small cutoff leads to slow quasiparticles with predominantly reflective scattering, and dynamics close to the standard semiclassical prediction of Subsection IV.1. For large cutoffs, on the other hand, the majority of the quasiparticles is ultrarelativistic and collide mostly by transmissive scatterings. Varying $E_{0},$ the semi-semiclassical approach interpolates between these two limiting cases.

V Perturbative quench theory

In this section, we develop a perturbative description of the quench in terms of an effective field theory, and compare its implications with the results of our TEBD simulations and the semi-semiclassical interpretation thereof. Certain details of the derivations are given in Appendix C.

Our first step is to replace the original Hamiltonian (1) by a phenomenological quasiparticle Hamiltonian, and to express the post-quench Hamiltonian as

[TABLE]

Here the operators $b^{\dagger}_{q,\sigma}$ create quasiparticles of momentum $q$, energy $\omega_{q}$, and spin $\sigma=\pm 1,0$. As discussed earlier, switching off $J_{2}$ creates singlet pairs of quasiparticles with opposite momenta. This is implemented within our effective field theory by the termfootnote_effective_Hamiltonian

[TABLE]

with the time dependent dispersion and pair creation amplitudes

[TABLE]

Here $\lambda(t)$ describes the time-dependent profile of the quench with $\lambda(t\rightarrow-\infty)=1$ and $\lambda(t\rightarrow\infty)=0.$ The precise momentum dependence of the couplings $\Delta_{q}$ and $g_{q}$ depends on microscopic details of the Hamiltonian $\hat{H}_{0}$ and the perturbing quench operator. At leading order in perturbation theory, the momentum distribution of quasiparticles is determined by the quasiparticle dispersion $\omega_{q}$, the specific form of $\lambda(t)$, and by the pair creation amplitude, $g_{q}$. This latter typically vanishes in the limit $q\to 0$ linearly, reflecting the fact that low energy quasiparticles behave as hard-core bosons.

Being quadratic, the simple model outlined above can be treated analytically, and can be used to determine the momentum distribution of the quasiparticles analytically. For small and smooth quenches we find

[TABLE]

where tilde denotes the Fourier transform. Note that $\lambda(t)$ does not depend on the amplitude of the quench, only on its duration and shape. Quite naturally, the quasiparticle density is proportional to the squared magnitude of the quench. Consequently, the collision rate $1/\tau$ defined in Eq. (12) is also just proportional to $J_{2}^{2}$.

This prediction is checked in Fig. 9.a, where we plot the inverse collision times extracted from TEBD as functions of $J_{2}/J_{1}$ for three different quench speeds. The quadratic dependence is valid up to $J_{2}/J_{1}\approx 1/3$ which is just the so-called AKLT point AKLT , suggesting that perturbation theory holds for quenches even up to this size. Also, slower, more adiabatic quenches create less quasiparticles, reflected in a lower collision rate.

As the quench magnitude $J_{2}$ only appears in the prefactor in Eq. (23), the normalized velocity distribution is predicted to be independent of it in this perturbative regime. The energy density should therefore be proportional to the quasiparticle density with a proportionality factor depending on the shape of $\lambda(t)$. As a consequence, the product of the collision time and the energy density $\rho_{\text{e}}$ should be independent of the size of the quench if perturbation theory holds. This prediction is tested in Fig. 9(b), where we plot the product $\rho_{\text{e}}\tau$ against $J_{2}/J_{1}$ for various quench durations. Here the quasiparticle energy density was measured in the TEBD simulation by taking the difference of the post-quench and the vacuum energy densities. It can be seen that the product is indeed approximately constant up to $J_{2}/J_{1}\approx 1/3,$ i.e. in the domain of perturbation theory identified above.

In order to make a more quantitative comparison between TEBD and our semiclassical method, we specify the shape of the quench as

[TABLE]

We furthermore assume a simple linear dependence of $g_{q}$ around $q=0$ (see Appendix C) and use a relativistic dispersion relation $\omega_{q}$ with the microscopic parameters given in Eq. (9). With these approximations, we obtain the quasiparticle momentum distribution used in our semi-semiclassical simulations, Eq. (19) with $E_{0}=t_{Q}$,

[TABLE]

With this distribution, we can compute the product of the energy density and the collision time for different quench durations $t_{Q}.$ The results are shown in Fig. 10 together with the TEBD data for quenches with $J_{2}/J_{1}=0.24.$ There is a reasonable agreement for not very fast quenches. This demonstrates convincingly that the semiclassical picture based on relativistic quasiparticles with perturbatively computed momentum distribution provides a valid qualitative and quantitative description of the post-quench dynamics.

Let us close this Section by verifying the consistency condition (10) of the semi-semiclassical approach. Even though we cannot measure the particle density directly, we can set an upper bound on it using the energy density:

[TABLE]

This implies that the semiclassical condition (10) is certainly satisfied if

[TABLE]

holds. We checked that this inequality was always satisfied in our simulations, except for the largest sudden quench with $J_{2}/J_{1}=0.36$, where $c\rho_{E}/\Delta^{2}\approx 1.2$ was found. However, for sudden quenches (26) strongly overestimates the particle density, i.e. (10) is expected to be valid even in this largest sudden quench.

VI Summary

In this work, we have studied numerically the statistics of spin transfer after a quantum quench in the Haldane chain, and compared it with the predictions of the semi-semiclassical approach, developed in Ref. PascuMarciGergo2017, . As we demonstrate, the probability distribution $P(S,t)$ of transferred spins, extracted from non-Abelian TEBD quench simulations, contains a surprising amount of information: in addition to detecting the presence of topologically protected end states AKLT ; Kennedy1990 , it gives insight to the internal spin structure and velocity distribution of quasiparticles as well as to the nature of their collisions.

Application of SU(2) symmetries allows us to reach sufficiently large bond dimensions in the range of 2000 multiplets corresponding to 10,000 states Werner_future , and to reach long enough simulation times to overlap with the range of validity of semi-semiclassics, and to observe effects related to quasiparticle collisions. We find that spin distributions, as computed through our full TEBD simulations, are in perfect agreement with those determined from the semi-semiclassical approach: the short time behavior of $P(S,t)$ reveals the presence of spin $S=1$ quasiparticles traveling ballistically, while at longer times we enter a collision-dominated regime, where collisions are neither reflective SachdevDamle1997 ; RappZarand ; Evangelisti2013 ; Kormos_Zarand_PRE_2016 nor transmissive Altshuler2006 . In this regime, it appears to be absolutely necessary to incorporate the coherent time evolution of the spin wave function, and to follow its quantum mechanical evolution, as performed by the hybrid semi-semiclassical approach.

For a sudden quench, the functions $P(S,t)$ are found to be universal functions of $t/\tau$, with $\tau$ the collision time. This universal scaling can also be simply explained within the semi-semiclassical theory: the size of the quantum quench has a direct influence on the quasiparticle density, $\varrho\sim J_{2}^{2}$, and thus the scattering rate, $1/\tau\sim\rho$, however, it does not change the velocity distribution of the quasiparticles and the structure and distribution of scattering matrices. Therefore, the statistics of collisions and the statistical time evolution of the wave function is independent of the size of the quench if time is measured in units of collision time.

The precise velocity distribution does depend, however, on details of the quench protocol and the quench time $t_{Q}.$ Correspondingly, for finite time quenches, we observe an explicit dependence of the functions $P(S,t)$ on the quench time (quench protocol). This explicit dependence is very well captured by simply changing the energy cut-off of the quasiparticles’ velocity distribution, $E_{0}\sim 1/t_{Q}$. This latter correspondence, and the velocity distributions used in our semi-semiclassical simulations, have both been motivated by a perturbative field theoretical quantum quench theory, presented in Section V.

The extension of perturbative quench regime as well as the range of validity of our semi-semiclassical approach appear to be astonishingly large. First of all, $1/\tau$, directly proportional to the density of quasiparticles is found to scale as $\sim J_{2}^{2}$ up to quench sizes of $J_{2}/J_{1}\approx 0.3,$ a value very close to the AKLT point at $J_{2}/J_{1}=1/3.$ Second, our assumption of a quench size independent velocity distribution is also verified by direct measurements of the quasiparticles’ energy density on the Haldane chain. Also quite remarkably, our perturbative field theoretical calculation assuming a simple relativistic quasiparticle spectrum yields a quantitative estimate for the dimensionless energy density $\rho_{e}\tau\sim\rho_{e}/\rho$, in very good agreement with our full TEBD spin chain simulations for quench times $t_{Q}\gtrsim 0.5/J_{1}$. This is, again, quite astonishing, since sudden quenches generate many quasiparticles presumably outside the regime of our effective field theory. Apparently, however, quasiparticles close to the gap appear to play the dominant role in the spin transport studied here.

The calculations presented here have thus two general conclusions: on the one hand, they show that semi-semiclassical calculations are simple methods that are able to capture important long-time features of non-equilibrium quantum dynamics at time scales much above the microscopic time scales. On the other hand, our studies reveal the power and richness of full counting statistics and full distributions, such as $P(S,t)$, containing precious information on correlations, equilibration and dynamics Gustavsson2006 ; Lamacraft2008 ; Hofferberth2008 ; Batalhao2014 ; Izabella2017 ; Collura2017 .

Acknowledgements. This research has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) through Grant Nos. SNN118028, K120569, and the Hungarian Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017- 00001). C.P.M. aknowledges support from the Romanian National Authority for Scientific Research and Innovation (UEFISCDI) through Grant No. PN-III-P4-ID-PCE-2016-0032. M.K. was also supported by a “Prémium” postdoctoral grant of the Hungarian Academy of Sciences. Ö. L. also acknowledges financial support from the Alexander von Humboldt foundation.

Appendix A NA-TEBD

The MPS based calculations of the time-evolution – both the microscopic simulations of the $S=1$ chain and the hybrid semi-semiclassical simulations – were performed by a variant of the standard TEBD algorithm Vidal2007 that exploits the presence of the non-Abelian SU(2) symmetry; we call our variant of the algorithm NA-TEBD, and our code is similar to the approaches described in Refs. Singh_Vidal2010, ; Wilson_NRG, ; Pascu_SU3_NRG, ; Weichselbaum2012, , but by the introduction of the $\eta$-legs the generalization to other non-Abelian symmetries becomes easier Werner_future ; footnote_on_outer_multiplicity .

In the TEBD algorithm, the time-evolution of the MPS wave function is simulated by consecutive application of two-site unitary evolvers Vidal2007 , based on the Trotter–Suzuki decomposition Trotter ; Suzuki . Considering now the SU(2) symmetric case, the evolvers like $U_{\sigma_{l}^{\prime}\sigma_{l+1}^{\prime}}^{\sigma_{l}\sigma_{l+1}}$, acting on the neighboring sites $l$ and $l+1$, preserve the SU(2) symmetry, and one naturally would like to exploit this property. Formally, one should contract the legs ${\sigma_{l},\sigma_{l+1}}$ of the tensor $U$ with the ones of the Clebsch–Gordan tensors in the two-layer NA-MPS (See Fig. 3 and Fig. 11). However, if the SU(2) symmetry is preserved by the two-site evolver, then after the application of it, the Clebsch-layer of the NA-MPS should remain the same, i.e. only the tensors $A^{[l]}$ and $A^{[l+1]}$ change. Using the standard orthogonality relations of the Clebsch–Gordan coefficients, we can determine the “symmetric” version of the two-site evolver that can be directly applied on the legs $\eta_{l}$ and $\eta_{l+1}$ of the tensors $A^{[l]}$ and $A^{[l+1]}$. This “symmetric” evolver is defined as (see also Fig.11a)

[TABLE]

The resulting tensor can be determined and stored before the simulation procedure, and later can be directly applied on the upper layer of the NA-MPS. It is important to remark that the blocks of the resulting tensor $\left[U_{\mathrm{sym}}\right]_{\eta_{l}^{\prime}\eta_{l+1}^{\prime}}^{\eta_{l}\eta_{l+1}}$ are specified by eight irreducible representation (“total-spin”) indices that are $\{J_{l-1},S_{l},J_{l},S_{l+1},J_{l+1},S_{l}^{\prime},J_{l}^{\prime},S_{l+1}^{\prime}\}$.

After the determination of $U_{\mathrm{sym}}$ the standard TEBD algorithm as described in Ref. Vidal2007, can be used to evolve the upper layer of the NA-MPS without any relevant modification as shown in Fig.11b. In other words, the layer of Clebsch–Gordan coefficients does not appear during the simulations, making our NA-TEBD code fast and efficient.

Numerical details. The rapidly growing truncation error (the weight of discarded Schmidt states) that spoils the accuracy of TEBD simulations leads to a short cutoff time beyond which results get unreliable. In our simulations we defined this cutoff time where the truncation error reached $10^{-7}$. To reach long enough times ($t\gtrsim\tau$), the MPS bond-dimension was set up to $M_{SU(2)}=2000$, i.e. 2000 multiplets (around 10,000 states) with the largest Schmidt values were kept.

Appendix B The relativistic S-matrix

In the hybrid semi-semiclassical simulations the relativistic S-matrix of the low-energy field theory of the $S=1$ Heisenberg model (namely, the $O(3)$ non-linear sigma model) was used to describe the collision events. We used the same S-matrix as in Ref. Marci_Pascu_Gergo_NESS, , but we transformed it in the more convenient form

[TABLE]

Here the (incoming and outgoing) basis states are $\{\left|++\right\rangle,\left|+0\right\rangle,\left|+-\right\rangle,\left|0+\right\rangle,\left|00\right\rangle,\left|0-\right\rangle,\left|-+\right\rangle,\left|-0\right\rangle,\left|--\right\rangle\}$, where $\{+,0,-\}$ stand for the eigenvalues $S^{z}=\{1,0,-1\}$ of the particles. In the two-particle basis states we used the convention that the left and right spin-state always stand for the particles in the left and right position, both in the incoming and outgoing states. The values $\sigma_{i}$ depend on the relative rapidity $\theta_{\mathrm{rel}}=\theta_{1}-\theta_{2}$ of the colliding particles as

[TABLE]

Appendix C Perturbative description of the quench

In this appendix we provide some details on the quench dynamics within our effective field theory, presented in Section V.

We shall assume in what follows that our initial quench state is a squeezed state Fioretto2010 ; Horvath2016 . Using that our time dependent Hamiltonian (21) is quadratic, one can show that the wave function remains a squeezed state upon time evolution,

[TABLE]

The exponential form corresponds to independent pairs of quasiparticles.

Our task is to determine the full time evolution of the functions $K_{q}(t)$ throughout a quantum quench. First we observe that different momentum modes decouple as $\left|\Psi(t)\right\rangle=\prod_{q>0}\left|\Psi_{q}(t)\right\rangle$, with

[TABLE]

The time-dependent Schrödinger equation leads to an equation for the pair creation amplitude $K_{q},$

[TABLE]

Within first order perturbation theory we can neglect the nonlinear $K_{q}(t)^{2}$ term in this equation. Using the initial values of the parameters, $\omega_{q}(t\rightarrow-\infty)=\omega_{q}+J_{2}\Delta_{q}$ and $g_{q}(t\rightarrow-\infty)=J_{2}g_{q}$, and exploiting that in the limit $t\rightarrow-\infty$ we initialize the system in its stationary vacuum state, we find the following initial conditions,

[TABLE]

In general, one needs to solve Eq. (33) numerically. However, we can easily obtain its approximate analytical solution up to leading order in $J_{2}$,

[TABLE]

Using this expression we find that the density of quasiparticles created by the quench is just given by Eq. (23).

For the particular quench profile

[TABLE]

the Fourier transformation can be performed analytically with the result

[TABLE]

For slow enough quenches the denominator cuts off the momentum distribution at low momenta, so we can substitute $g_{q}$ with its small momentum behavior, $g_{q}\approx g\;|q|,$ and we obtain the momentum distribution (25).