Chaotic dynamical ferromagnetic phase induced by non-equilibrium quantum fluctuations

TL;DR

This paper explores how quantum fluctuations influence the stability of dynamical phase transitions in a fully-connected quantum Ising model, revealing a chaotic phase with sensitive dependence on initial conditions.

Contribution

It demonstrates that increasing quantum fluctuations can induce a chaotic dynamical phase, extending the understanding of non-equilibrium quantum phase behavior beyond mean-field predictions.

Findings

Quantum fluctuations cause the dynamical critical point to fan out into a chaotic phase.

The chaotic phase exhibits strong sensitivity to parameters and initial conditions.

The phenomenon is likely general for systems with broken discrete symmetry.

Abstract

We investigate the robustness of a dynamical phase transition against quantum fluctuations by studying the impact of a ferromagnetic nearest-neighbour spin interaction in one spatial dimension on the non-equilibrium dynamical phase diagram of the fully-connected quantum Ising model. In particular, we focus on the transient dynamics after a quantum quench and study the pre-thermal state via a combination of analytic time-dependent spin-wave theory and numerical methods based on matrix product states. We find that, upon increasing the strength of the quantum fluctuations, the dynamical critical point fans out into a chaotic dynamical phase within which the asymptotic ordering is characterised by strong sensitivity to the parameters and initial conditions. We argue that such a phenomenon is general, as it arises from the impact of quantum fluctuations on the mean-field out of equilibrium dynamics of any system which exhibits a broken discrete symmetry.

Equations

$$H=-\frac{\lambda }{N}\sum _{i,j=1}^N{\sigma }_i^x{\sigma }_j^x-g\sum _{i=1}^N{\sigma }_i^z-J\sum _{i=1}^N{\sigma }_i^x{\sigma }_{i+1}^x,$$

$$\frac{{\sigma }_i^X}{2}\simeq \sqrt{s}{q}_i,\frac{{\sigma }_i^Y}{2}\simeq \sqrt{s}{p}_i,\frac{{\sigma }_i^Z}{2}=s-\frac{q_i^2+p_i^2-1}{2}.$$

$$\{\begin{array}{rl}\frac{d\theta }{dt}=& 4\left[\bar{\lambda }\rho (t)-J{\delta }^{pp}(t)\right]\sin \theta \cos \varphi \sin \varphi \\ & +4J{\delta }^{qp}(t)\cos \theta \sin \theta {\cos }^2\varphi ,\\ \frac{d\varphi }{dt}=& -2g+4\left[\bar{\lambda }\rho (t)-J{\delta }^{qq}(t)\right]\cos \theta {\cos }^2\varphi \\ & +4J{\delta }^{qp}(t)\sin \varphi \cos \varphi ,\end{array}$$

$$\begin{array}{rl}{\Delta }_k^{qq}(t)& \equiv ⟨{\tilde{q}}_k(t){\tilde{q}}_{-k}(t)⟩,{\Delta }_k^{pp}(t)\equiv ⟨{\tilde{p}}_k(t){\tilde{p}}_{-k}(t)⟩,\\ {\Delta }_k^{qp}(t)& \equiv ⟨{\tilde{q}}_k(t){\tilde{p}}_{-k}(t)+{\tilde{p}}_k(t){\tilde{q}}_{-k}(t)⟩/2.\end{array}$$

$$ϵ(t)\equiv \frac{1}{N/2}\sum _{k\ne 0}({\Delta }_k^{qq}+{\Delta }_k^{pp}-1)/2,$$

$$H=-g{\tilde{\sigma }}_0^z-\frac{\bar{\lambda }}{N}{\left({\tilde{\sigma }}_0^x\right)}^2-\frac{J}{N}\sum _{k\ne 0}\cos k{\tilde{\sigma }}_k^x{\tilde{\sigma }}_{-k}^x,$$

$$\hat{Z}\equiv \left(\begin{array}{c}\sin \theta \cos \varphi \\ \sin \theta \sin \varphi \\ \cos \theta \end{array}\right);\hat{X}\equiv \left(\begin{array}{c}\cos \theta \cos \varphi \\ \cos \theta \sin \varphi \\ -\sin \theta \end{array}\right);\hat{Y}\equiv \left(\begin{array}{c}-\sin \varphi \\ \cos \varphi \\ 0\end{array}\right).$$

$$V{\sigma }_j^x{V}^{†}=\hat{X}\cdot {\vec{\sigma }}_j\equiv {\sigma }_j^X,V{\sigma }_j^y{V}^{†}=\hat{Y}\cdot {\vec{\sigma }}_j\equiv {\sigma }_j^Y,V{\sigma }_j^z{V}^{†}=\hat{Z}\cdot {\vec{\sigma }}_j\equiv {\sigma }_j^Z,$$

$$\frac{d}{dt}\sigma_{j}^{\alpha}=\frac{1}{i}\big{[}\sigma_{j}^{\alpha},\widetilde{H}\big{]},$$

$$\tilde{H}\equiv H+iV{\dot{V}}^{†},$$

$$iV{\dot{V}}^{†}=-s\vec{\omega }(t)\cdot \sum _j{\vec{\sigma }}_j,$$

$${\vec{\sigma }}_j=\hat{X}{\sigma }_j^X+\hat{Y}{\sigma }_j^Y+\hat{Z}{\sigma }_j^Z.$$

$$\begin{array}{rl}\frac{\tilde{H}}{N}=& -g\left[\left(\hat{X}\cdot \hat{z}\right)\frac{{\tilde{\sigma }}_0^X}{N}+\left(\hat{Y}\cdot \hat{z}\right)\frac{{\tilde{\sigma }}_0^Y}{N}+\left(\hat{Z}\cdot \hat{z}\right)\frac{{\tilde{\sigma }}_0^Z}{N}\right]-\bar{\lambda }{\left[\left(\hat{X}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_0^X}{N}+\left(\hat{Y}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_0^Y}{N}+\left(\hat{Z}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_0^Z}{N}\right]}^2\\ & -J\sum _{k\ne 0}\cos k\left[\left(\hat{X}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_k^X}{N}+\left(\hat{Y}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_k^Y}{N}+\left(\hat{Z}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_k^Z}{N}\right]\left[\left(\hat{X}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_{-k}^X}{N}+\left(\hat{Y}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_{-k}^Y}{N}+\left(\hat{Z}\cdot \hat{x}\right)\frac{{\tilde{\sigma }}_{-k}^Z}{N}\right]\\ & +\sin \theta s\dot{\varphi }\frac{{\tilde{\sigma }}_0^X}{N}-s\dot{\theta }\frac{{\tilde{\sigma }}_0^Y}{N}-\cos \theta s\dot{\varphi }\frac{{\tilde{\sigma }}_0^Z}{N},\end{array}$$

$$\{\begin{array}{rl}{\sigma }_i^X& =\frac{q_i}{\sqrt{s}}+\mathcal{O}\left((q_i,p_i/\sqrt{s}{)}^3\right),\\ {\sigma }_i^Y& =\frac{p_i}{\sqrt{s}}+\mathcal{O}\left((q_i,p_i/\sqrt{s}{)}^3\right),\\ {\sigma }_i^Z& =1-\frac{b_i^{†}{b}_i}{s}\equiv 1-\frac{q_i^2+p_i^2-1}{2s},\end{array}$$

$$\{\begin{array}{rl}\frac{{\tilde{\sigma }}_k^X}{N}=& \frac{{\tilde{q}}_k}{\sqrt{Ns}}+\mathcal{O}\left((\tilde{q},\tilde{p}/\sqrt{Ns}{)}^3\right),\\ \frac{{\tilde{\sigma }}_k^Y}{N}=& \frac{{\tilde{p}}_k}{\sqrt{Ns}}+\mathcal{O}\left((\tilde{q},\tilde{p}/\sqrt{Ns}{)}^3\right),\\ \frac{{\tilde{\sigma }}_k^Z}{N}=& {\delta }_{k,0}-\sum _{k^{\prime }}\frac{{\tilde{q}}_{k^{\prime }}{\tilde{q}}_{k-k^{\prime }}+{\tilde{p}}_{k^{\prime }}{\tilde{p}}_{k-k^{\prime }}-{\delta }_{k,0}}{2(\sqrt{Ns}{)}^2}.\end{array}$$

$${\vec{\Sigma }}_{\text{tot}}\equiv s\sum _j{\vec{\sigma }}_j,$$

$$\big{\lvert}\vec{\Sigma}_{\text{tot}}\big{\rvert}^{2}=\left(Ns-\sum_{k\neq 0}\tilde{b}_{k}^{\dagger}\tilde{b}_{k}\right)\left(Ns-\sum_{k\neq 0}\tilde{b}_{k}^{\dagger}\tilde{b}_{k}+1\right),$$

$$\big{\langle}\tilde{\sigma}_{0}^{X}(t)\big{\rangle}=\big{\langle}\tilde{\sigma}_{0}^{Y}(t)\big{\rangle}=0\quad\text{or, equivalently,}\quad S_{X}(t)=S_{Y}(t)=0.$$

$$\sum_{k\neq 0}\langle n_{k}(t)\rangle=\sum_{k\neq 0}\Big{\langle}\frac{\tilde{q}_{k}(t)\tilde{q}_{-k}(t)+\tilde{p}_{k}(t)\tilde{p}_{-k}(t)-1}{2}\Big{\rangle}\ll Ns,$$

$$n_k\equiv {\tilde{b}}_k^{†}{\tilde{b}}_k.$$

$$\frac{d}{dt}\big{\langle}\tilde{\sigma}_{0}^{X}(t)\big{\rangle}=\frac{d}{dt}\big{\langle}\tilde{\sigma}_{0}^{Y}(t)\big{\rangle}=0,\quad\mbox{with}\quad\big{\langle}\tilde{\sigma}_{0}^{X}(t=0)\big{\rangle}=\big{\langle}\tilde{\sigma}_{0}^{Y}(t=0)\big{\rangle}=0.$$

$$\epsilon(t)\equiv\frac{1}{Ns}\sum_{k\neq 0}\big{\langle}n_{k}(t)\big{\rangle}=\frac{1}{Ns}\sum_{k\neq 0}\bigg{\langle}\frac{\tilde{q}_{k}(t)\tilde{q}_{-k}(t)+\tilde{p}_{k}(t)\tilde{p}_{-k}(t)-1}{2}\bigg{\rangle},$$

$$\frac{{\tilde{\sigma }}_0^Z}{N}=1-\frac{1}{Ns}\sum _{k\ne 0}\frac{{\tilde{q}}_k(t){\tilde{q}}_{-k}(t)+{\tilde{p}}_k(t){\tilde{p}}_{-k}(t)-1}{2},$$

$$\frac{d}{dt}ϵ=0,$$

$$U_1$$

$$U_2$$

$$U_3$$

$$\left\{\begin{split}\frac{d}{dt}\theta=&+4\bar{\lambda}(1-\epsilon)\sin\theta\cos\phi\sin\phi\\ &-4J\bigg{(}\frac{1}{Ns}\sum_{k\neq 0}\cos k\;\left\langle\tilde{p}_{k}\tilde{p}_{-k}\right\rangle\bigg{)}\sin\theta\cos\phi\sin\phi\\ &+4J\Bigg{(}\frac{1}{Ns}\sum_{k\neq 0}\cos k\frac{\left\langle\tilde{q}_{k}\tilde{p}_{-k}+\tilde{p}_{k}\tilde{q}_{-k}\right\rangle}{2}\Bigg{)}\cos\theta\sin\theta\cos^{2}\phi,\\ \frac{d}{dt}\phi=&-2g+4\bar{\lambda}(1-\epsilon)\cos\theta\cos^{2}\phi\\ &-4J\bigg{(}\frac{1}{Ns}\sum_{k\neq 0}\cos k\;\left\langle\tilde{q}_{k}\tilde{q}_{-k}\right\rangle\bigg{)}\cos\theta\cos^{2}\phi\\ &+4J\Bigg{(}\frac{1}{Ns}\sum_{k\neq 0}\cos k\frac{\left\langle\tilde{q}_{k}\tilde{p}_{-k}+\tilde{p}_{k}\tilde{q}_{-k}\right\rangle}{2}\Bigg{)}\sin\phi\cos\phi,\end{split}\right.$$

$$\begin{array}{rl}{\Delta }_k^{qq}(t)& \equiv ⟨{\tilde{q}}_k(t){\tilde{q}}_{-k}(t)⟩,\\ {\Delta }_k^{pp}(t)& \equiv ⟨{\tilde{p}}_k(t){\tilde{p}}_{-k}(t)⟩,\\ {\Delta }_k^{qp}(t)& \equiv \frac{1}{2}⟨{\tilde{q}}_k(t){\tilde{p}}_{-k}(t)+{\tilde{p}}_k(t){\tilde{q}}_{-k}(t)⟩.\end{array}$$

$$\{\begin{array}{rl}\frac{d}{dt}{\tilde{q}}_k=& +4\bar{\lambda }{\cos }^2\varphi {\tilde{p}}_k-4J\cos k{\sin }^2\varphi {\tilde{p}}_k+4J\cos k\cos \theta \cos \varphi \sin \varphi {\tilde{q}}_k,\\ \frac{d}{dt}{\tilde{p}}_k=& -4\bar{\lambda }{\cos }^2\varphi {\tilde{q}}_k+4J\cos k{\cos }^2\theta {\cos }^2\varphi {\tilde{q}}_k-4J\cos k\cos \theta \cos \varphi \sin \varphi {\tilde{p}}_k,\end{array}$$

$$\{\begin{array}{rl}\frac{d}{dt}{\Delta }_k^{qq}=& 8J\cos k\cos \theta \cos \varphi \sin \varphi {\Delta }_k^{qq}+8\left(\bar{\lambda }{\cos }^2\varphi -J\cos k{\sin }^2\varphi \right){\Delta }_k^{qp},\\ \frac{d}{dt}{\Delta }_k^{qp}=& -4\left(\bar{\lambda }{\cos }^2\varphi -J\cos k{\cos }^2\theta {\cos }^2\varphi \right){\Delta }_k^{qq}+4\left(\bar{\lambda }{\cos }^2\varphi -J\cos k{\sin }^2\varphi \right){\Delta }_k^{pp},\\ \frac{d}{dt}{\Delta }_k^{pp}=& -8\left(\bar{\lambda }{\cos }^2\varphi -J\cos k{\cos }^2\theta {\cos }^2\varphi \right){\Delta }_k^{qp}-8J\cos k\cos \theta \cos \varphi \sin \varphi {\Delta }_k^{pp}.\end{array}$$

Full text

Introduction — Throwing dice on the table floor is a prototypical random (or pseudorandom) process Strzalko et al. (2009). Its aleatory nature is a consequence of a few ingredients: the die, initially out of equilibrium, dissipates its energy rolling on the table, and hence relaxes onto one of few possible equilibrium configurations. In this work we show that the same ingredients play an important role in the physics of quantum many-body systems undergoing a Dynamical Quantum Phase Transition (DQPT), leading to the emergence of an intriguing chaotic dynamical phase once the collective dynamics of these systems gets damped by quantum fluctuations.

DQPTs are among the most interesting phenomena occurring in quantum many-body systems after a sudden change of the system parameters (quantum quench) Polkovnikov et al. (2011); *Lamacraft2012; *Eisert2015a, a type of process which can be realized both with ultra-cold gases Greiner et al. (2002a); *Greiner2002b; *Bloch2008; *Trotzky2008; *Cheneau2012; *Kauf16; *Jorg15; *rauer17 and trapped ions Hess et al. (2017). Such DQPTs Sciolla and Biroli (2011); Smacchia et al. (2015); Maraga et al. (2016); Halimeh et al. (2017); *Halimeh2017; *Gambassi2011 are characterised by the vanishing of a non-equilibrium order parameter (accompanied by critical scaling behaviour Sciolla and Biroli (2013); Chiocchetta et al. (2015)) and to be distinguished from those signalled by non-analyticities in the temporal evolution of the Loschmidt echo Heyl et al. (2013) (see Ref. Zunkovic et al. (2016) for connections between the two notions). They not only provide a genuine instance of classical and quantum criticality out of equilibrium Sciolla and Biroli (2013); Chandran et al. (2013); *Chiocchetta2015; Chiocchetta et al. (2017), but demonstrate also the emergence of intermediate stages of relaxation with nontrivial time-dependent fluctuations and dynamics Smacchia et al. (2015). A DQPT separates “phases” characterised by qualitatively different quasi-stationary states Sciolla and Biroli (2011); Gambassi and Calabrese (2011), anomalous coarsening Maraga et al. (2015), aging Janssen et al. (1989); *Calabrese2005; Maraga et al. (2015); Chiocchetta et al. (2017), as well as by a non-trivial dynamical evolution of observables and their fluctuations Chandran et al. (2013); Smacchia et al. (2015); Maraga et al. (2016); Zunkovic et al. (2016); Halimeh et al. (2017); *Halimeh2017. Due to the lack of spatial and temporal collective scales upon approaching a DQPT, they display features reminiscent of equilibrium critical points.

DQPTs are expected to be strongly affected by quantum fluctuations: recent investigations beyond mean-field approximations Schiró and Fabrizio (2011); Sandri et al. (2012); Sciolla and Biroli (2013); Peronaci et al. (2015) showed that these fluctuations influence, e.g., the early stages of the evolution Maraga et al. (2015); Chiocchetta et al. (2015, 2017). In this work we demonstrate a more dramatic effect of fluctuations on the dynamics of the order parameter, which induces a qualitative modification of the dynamical phase diagram, in particular close to the dynamical critical point. We study the non-equilibrium dynamics of an infinite-range (mean-field) ferromagnetic system perturbed by additional short-range interaction terms, which rule the strength of quantum fluctuations. We show that the dynamical phases are robust, whereas the impact of non-equilibrium quantum fluctuations makes the dynamical critical point open up in a novel chaotic dynamical phase where the dynamics are reminiscent of that of a coin toss: The asymptotic stationary state displays a finite magnetization whose positive or negative sign is highly sensitive to initial conditions and system parameters, as we show in Fig.1.

The model — In this work, for the sake of definiteness, we focus on a fully-connected quantum Ising ferromagnet in a transverse magnetic field $g$, in the presence of additional nearest-neighbor couplings in one spatial dimension, governed by the Hamiltonian

[TABLE]

where $\sigma^{\alpha}_{i}$ are the standard Pauli matrices at lattice site $i$. In the limit $J\to 0$, $H$ maps to the exactly solvable Lipkin-Meshkov-Glick (LMG) model Lipkin et al. and displays both a quantum critical point in equilibrium Dutta and Bhattacharjee (2001); *knap13 at $g=2\lambda$ and a DQPT after a quench Sciolla and Biroli (2011); Zunkovic et al. (2016), with the longitudinal global magnetization $S_{x}(t)$ — with $S_{\alpha}\equiv\langle\sum_{i}\sigma_{i}^{\alpha}\rangle/N$ — being the dynamical order parameter of the DPT. For example, for quenches starting from the ferromagnetic ground state at $g_{0}=0$, $S_{x}(t)$ evolves periodically with a period set by the post-quench values $g$ and $\lambda$ of the couplings. In particular, the time average $\bar{S}_{x}=\lim_{T\to\infty}\int^{T}_{0}{\rm d}t\,S_{x}(t)/T$ vanishes for $g>\lambda$ (because the oscillations of $S_{x}(t)$ are symmetric around zero) while it does not for $g<\lambda$ (because the oscillations do not change the sign of $S_{x}$), corresponding to the dynamically paramagnetic and ferromagnetic “phases”, respectively. At the dynamical critical point, $g=\lambda$, the order parameter decays exponentially to zero with $S_{x}(t)\sim{\rm e}^{-gt}$ for $t\gg g^{-1}$. Within mean-field theory, which is an exact treatment of the LMG model in the thermodynamic limit Zunkovic et al. (2016), the DQPT can be rationalized Gambassi and Calabrese (2011); Sciolla and Biroli (2013); Zunkovic et al. (2016) in terms of the motion of a classical particle with position $S_{x}(t)$ in an effective, double-well even potential ${\cal U}(S_{x})$. If $S_{x}(t=0)$ is such that ${\cal U}(S_{x}(0))>{\cal U}(0)$, then $S_{x}(t)$ explores both wells and $\bar{S}_{x}=0$; otherwise the motion is localized within one well and $\bar{S}_{x}\neq 0$.

In this work we study how this mean-field non-equilibrium phase diagram is affected by quantum fluctuations. While for $J=0$ all spins perform a coherent collective motion, turning on a short-range perturbation is expected to damp the persistent classical oscillations of $S_{x}(t)$, altering the features of the mean-field evolution and inducing relaxation towards a stationary and eventual thermal state.

In order to address these questions we develop a spin-wave theory in the reference frame aligned with the instantaneous average total spin, with the spin-coherent state in this direction representing the instantaneous spin-wave vacuum. While for $J=0$ the length of the total spin is constantly maximal, i.e., $\lvert S(t)\rvert\equiv 1$, in the presence of a small short-range perturbation $J\neq 0$ a finite density $\epsilon(t)$ of spin-wave excitations is generated by the precessing collective spin, yielding $\lvert S(t)\rvert=1-\epsilon(t)$. As long as $\epsilon(t)\ll 1$, the non-linear, inelastic scattering among spin-waves is negligible and thermalization is expected to occur at longer times. Accordingly, the temporal regime with $\epsilon(t)\ll 1$, within which the mean-field motion receives correction from having $J\neq 0$ while keeping its non-equilibrium features (e.g., the DQPT), can be qualified as being pre-thermal, in analogy with similar cases Langen et al. (2016).

Outline of results — In the presence of quantum fluctuations, one would expect the collective motion of $S_{x}(t)$ to be damped by the generation of spin-wave excitations with a finite rate, leading to the breakdown of the approximation $\epsilon(t)\ll 1$. (Throughout the paper we fix energy units such that $\bar{\lambda}\equiv\lambda+J=1$.) We find, instead, that for small $J\lesssim 0.25$, $\epsilon(t)$ always saturates, implying that the dynamical paramagnetic and ferromagnetic phases indicated by A and B, respectively, in Fig. 1 are stable. In particular, $S_{x}(t)$ approximately oscillates with a period which is perturbatively close to the mean-field one. A numerical analysis based on the time-dependent variational principle (MPS-TDVP) Lubich et al. (2015); Haegeman et al. (2016) indicates that this stability extends to larger values of $J$ where the spin-waves density is no longer small. This implies that $H$ inherits the dynamical phase diagram of the classical case with $J=0$. However, the presence of spin-wave excitations makes the dynamical critical point at $\lambda=g$ for $J=0$, fan out in a chaotic dynamical ferromagnetic phase denoted by C in Fig. 1. Within C, the non-equilibrium quantum fluctuations in the form of spin-waves act effectively as a self-generated bath responsible for the localization of the system, initially with $S_{x}>0$, into one of the two wells of ${\cal U}$, with either sign of $S_{x}$ and for the possible hopping of the collective spin $S_{x}$ between the two of them; these processes are sketched as (a) and (b), respectively, in the inset of Fig. 1. The strong sensitivity of the long-time ferromagnetic ordering to the values of $g$, $\lambda$, and $J$ and of the initial data can be regarded as a signature of a collective chaotic behavior. The stability of this picture upon increasing $N$ in both the analytic and the numerical approaches leads us to conclude that such behavior carries over to the thermodynamic limit.

Time-dependent spin-wave theory — We now briefly outline the non-equilibrium spin-wave theory at the core of this work SM . We first introduce a time-dependent reference frame $\mathcal{R}=(\hat{X},\hat{Y},\hat{Z})$ in the spin space, with its $\hat{Z}$-axis following the collective motion of $\vec{S}(t)$. The change of frame is implemented by the time-dependent global rotation operator V\big{(}\theta(t),\phi(t)\big{)}=\exp\left(-i\phi(t)\sum_{i}\sigma_{i}^{z}/2\right)\exp\left(-i\theta(t)\sum_{i}\sigma_{i}^{y}/2\right) parameterized by the angles $\theta(t)$ and $\phi(t)$, which are eventually determined in such a way that $S_{X}(t)\equiv S_{Y}(t)\equiv 0$. For $J=0$, when $H$ is a function of the total spin only, this requirement translates into a closed pair of (classical) ordinary differential equations, the solution of which determines the evolution of the order parameter $S_{x}(t)=\sin\theta(t)\cos\phi(t)$ Zunkovic et al. (2016). For $J\neq 0$, the additional short-range interaction renders $H$ a function of not only the total spin, i.e., the $k=0$ Fourier mode of the spins, but also of all the $k$-modes of the spins, which now contribute to the dynamics. In order to make the equations of motion tractable and to set up a systematic expansion, we introduce the canonically conjugated spin-wave variables $q_{i}$ and $p_{i}$ at site $i$ with respect to the instantaneous $\hat{Z}$-axis via the Holstein–Primakoff (HP) transformation

[TABLE]

We then express all the spin operators in $H$, see Eq. (1), in terms of the spin-wave coordinates in Fourier space $\tilde{q}_{k}$, $\tilde{p}_{k}$, and retain up to quadratic terms in the spatial fluctuations modes $(\tilde{q}_{k},\tilde{p}_{k})$ with $k\neq 0$ (i.e., we neglect collisions among spin-waves). After averaging the Heisenberg equations of motion of the spins over the nonequilibrium state SM , we find that $\theta$ and $\phi$ evolve according to (recall $\bar{\lambda}\equiv\lambda+J$)

[TABLE]

where $\delta^{\alpha\beta}(t)\equiv\sum_{k\neq 0}\Delta^{\alpha\beta}_{k}\cos k/(Ns)$ with $\alpha,\beta\in\{p,q\}$ is the quantum “feedback” given by the correlation functions of the spin-waves,

[TABLE]

The relevance of these spin-wave excitations is controlled by the quantity

[TABLE]

i.e., by the total number of spin-waves divided by $N/2$. In Eq. (34), $\rho(t)=1-\epsilon(t)$ is the ratio between the expectation value of the modulus of the total spin of the system and its maximal value $N/2$, which is conserved by the dynamics only when $J=0$ Zunkovic et al. (2016). The evolution of $\Delta^{\alpha\beta}_{k}$ in Eq. (4) is ruled by a system of linear differential equations involving $\theta(t)$ and $\phi(t)$ SM . The quadratic approximation is justified as long as the density of excited spin-waves is small, i.e., $\epsilon(t)\ll 1$. For a quench starting from the spin coherent state fully polarized in the $\hat{x}$ direction, considered here (i.e., from $g=0$) the initial data of Eqs. (34) are $\theta(0)=\pi/2$, $\phi(0)=0$ with $\Delta^{qq}_{k}(0)=\Delta^{pp}_{k}(0)=1/2$, and $\Delta^{qp}_{k}(0)=0$ for $k\neq 0$; in particular, $\epsilon(0)=0$ (note that at time $t=0$ the mobile $\hat{Z}$-axis is aligned with the fixed $\hat{x}$ direction). Equation (34) includes the feedback terms $\delta^{\alpha\beta}(t)$ from quantum fluctuations, which both “dress” the value of $\bar{\lambda}$ and generate new terms of pure quantum origin in addition to the classical mean-field dynamics corresponding to $J=0$ in Eq. (34).

Nonequilibrium quantum phase diagram — Via a joint numerical integration of Eq. (34) and of the evolution equations (see Eqs. (26) in Ref. SM ) of $\Delta^{\alpha\beta}_{k}$ in Eq. (4), for a range of post-quench values of $g$ and $J$, we obtained the dynamical “phase” diagram portrayed in Fig. 1. In particular, for each integration, we compute the direction $(\theta(t),\phi(t))$ of the total spin $\vec{S}$ and the density $\epsilon(t)$ of spin-waves, verifying that the latter always settles around a small value at long times within the range of parameters considered here. Then, we compute the long-time average $\bar{S}_{x}$ of $S_{x}(t)$ and color the corresponding point in light yellow if $\bar{S}_{x}>0$, in orange if $\bar{S}_{x}=0$, and in blue if $\bar{S}_{x}<0$. The results of this procedure, in Fig. 1, shows that the two dynamical ferromagnetic and paramagnetic phases present for $J=0$ survive at $J>0$: for $g\ll\bar{\lambda}$ the order parameter has non-zero time-average (its value being perturbatively close to the mean-field one), while for $g\gg\bar{\lambda}$ it vanishes for all the values of $J$ within the range considered here. In particular, the persistent oscillations of $S_{x}$, characteristic of the mean-field solution, are not wiped out by the spin-wave bath, which does not produce a significant noise during the pre-thermal stage of dynamics and leaves the overall motion perturbatively close to perfect coherence for all observation times.

Near the dynamical transition point $g\simeq\bar{\lambda}$, the system becomes extremely sensitive to the non-equilibrium quantum fluctuations, realizing a peculiar intermediate phase, the existence of which is intimately related to the conservation of the energy. In a typical point of this region, the dynamics of $S_{x}$ is characterized by the processes illustrated by the inset of Fig. 1 and by Fig. 2: the decay from a transient paramagnetic behavior to one of the two ferromagnetic sectors, possibly followed by one or more hops between them. Typically, after an initial transient with the energy of the macroscopic total spin slightly above the barrier ${\mathcal{U}}(0)$ separating the two ferromagnetic wells of the effective potential ${\mathcal{U}}$, the production of spin-waves causes the dynamics to get trapped within one of the two. The system thus shows ferromagnetic order at long times, though it might occasionally hop to the opposite well, assisted by the absorption of energy from the spin-waves. The asymptotic sign of $S_{x}(t)$, and hence of $\bar{S}_{x}$, sensitively depends on the specific values of the parameters in a large part of this novel ferromagnetic region (indicated by C in Fig. 1), implying a collective chaotic character of the dynamics within it, as illustrated in Fig. 2. Unlike the quantum critical cone emanating from equilibrium quantum critical points at finite temperature Sachdev (2011), the boundaries of region C are expected to be sharp both towards the ferromagnetic and the paramagnetic phases, signalling two transitions expected to be characterized by diverging time scales.

The dynamical phases discussed here turn out to be robust against the perturbation of $J$ even for values $J\simeq 0.67$ at which the low-density spin-wave expansion need not be accurate. To show this, we performed numerical simulations by using a time-dependent variational principle on the matrix product state manifold (MPS-TDVP) Lubich et al. (2015); Haegeman et al. (2016), resulting in the evolution reported in Fig. 3; this approach allows us to explore the dynamics of $S_{x}(t)$ up to times of the order $\sim 60\bar{\lambda}^{-1}$. For $g\lesssim\bar{\lambda}$ we find a ferromagnetic region with $\bar{S}_{x}\neq 0$ of the same sign as the initial magnetisation $S_{x}(t=0)$. For large values of $g\gtrsim\bar{\lambda}$, instead, we find a paramagnetic phase with $\bar{S}_{x}=0$ SM , while for intermediate values, $\bar{S}_{x}$ does not vanish but it may have a sign opposite to that of $S_{x}(0)$; this observation is consistent with what observed at smaller values of $J$, see Fig. 3. In addition, in this regime, the final value of $\bar{S}_{x}$ may sensibly depend on the system size $N$: For $N\approx 100$ (see caption of Fig. 3) we observe $\bar{S}_{x}\neq 0$ of the same sign as $S_{x}(t=0)$, while for a slightly larger system $N=125$, $S_{x}(t)$ at long times has the opposite sign, which is eventually observed also in a system with $N=400$, as shown by the inset of Fig. 3. This is consistent with the sensitivity to the parameters predicted by the spin-wave approach, see Fig. 1. These numerical simulations of the exact quantum evolution fully confirm — and even extend to a larger region of the phase diagram — the scenario outlined by the time-dependent spin-wave theory, i.e., the robustness of the two dynamical phases and the emergence of a chaotic region in between.

Perspectives — In summary, the non-equilibrium quantum fluctuations due to spin-wave excitations modify qualitatively the mean-field phase diagram, turning the $J=0$ quantum critical point into a phase with unusual dynamical properties. The non-equilibrium spin-wave theory at the core of this work can be straightforwardly extended to a wide variety of spin systems, in higher dimensions, with other types of integrability breaking terms (of short or long-range character) or non-equilibrium protocols: a chaotic dynamical phase always arises whenever a mean-field system undergoing a ferromagnetic transition is subject to the impact of out-of-equilibrium quantum fluctuations Lerose et al. (2017). In addition, the phenomena discussed here could be within experimental reach, considering recent progress in realising spin models Simon and et al. (2011); *Meinert13; *Labuhn16 as well as in highlighting universal scaling behaviour close to dynamical critical points using cold gases Nicklas et al. (2015).

Acknowledgements — We thank M. Fabrizio for interesting discussions and A. Rosch for insightful comments on the manuscript. A. L. acknowledges hospitality from the University of Cologne. J. M. acknowledges support from the Alexander von Humboldt foundation. B. Ž. is supported by the ERC project OMNES.

Appendix A Supplemental Material:

Chaotic dynamical ferromagnetic phase

induced by non-equilibrium quantum fluctuations

In Section 1, a comprehensive exposition of the time-dependent spin-wave theory is provided. In Section 2, perturbative estimates for the quantum feedback terms near the mean-field dynamical critical point are presented. Finally, in Section 3, details about the numerical study of the full many-body non-equilibrium quantum evolution via TDVP based on matrix product states are discussed.

Appendix B 1. Time-dependent spin-wave theory

As outlined in the Letter, our strategy to study the effect of quantum fluctuations on the dynamics of the spin chain under consideration is based on expressing the quantum dynamical evolution in a time-dependent rotating reference frame ${\mathcal{R}}$ adapted to the instantaneous direction of the average total spin (see [A. Ruckriegel, et al., Phys. Rev. B 85, 054422 (2012)] for a similar approach). In ${\mathcal{R}}$ we can introduce the canonical coordinates corresponding to the spatial fluctuations of the spin field, i.e., to the spin-waves, and write the coupled equations of motion for both the total spin and the spin-waves. The evolution of the latter enters as a feedback into the dynamics of the former while the motion of the reference frame is determined self-consistently under the condition that it remains constantly aligned with the instantaneous average value of the total spin. Below we describe the main technical steps of this approach.

B.1 Time-dependent coordinates

As usual when dealing with spin-wave theory, it is convenient to consider a spin $\vec{\Sigma}$ of generic value $s$, instead of focussing directly on $s=1/2$, in order to keep track of the classical limit $s\to\infty$. Accordingly, we consider the original Hamiltonian (see Eq. (1) of the Letter), in which $\vec{\sigma}_{i}\equiv\vec{\Sigma}_{i}/s$ are now normalized spin-$s$ operators: they reduce to the standard Pauli matrices for $s=1/2$, while their $1/s$ normalization makes the classical limit $s\to\infty$ well-defined, assuming that time is measured in units of $2s$.

Assuming periodic boundary conditions along the chain, we introduce Fourier modes $\tilde{\vec{\sigma}}_{k}\equiv\sum_{j}e^{-ikj}\vec{\sigma}_{j}$, with wavevectors $k=2\pi n/N$ and $n=0,1,\dots,N-1$; accordingly, the Hamiltonian in Eq. (1) of the Letter becomes

[TABLE]

where $\bar{\lambda}\equiv\lambda+J$ is the effective value of $\lambda$ which has been renormalized by the term of the Hamiltonian corresponding to the Fourier mode $k=0$. Note that $\langle\tilde{\vec{\sigma}}_{k=0}\rangle=\langle\sum_{j}\vec{\sigma}_{j}\rangle=N\vec{S}$ is the normalized average total spin.

We now perform a time-dependent change of reference frame, introducing the right–handed Cartesian triple $(\hat{X},\hat{Y},\hat{Z})$ of unit vectors in the mobile frame ${\mathcal{R}}$, whose components in the original fixed frame $(\hat{x},\hat{y},\hat{z})$ are parameterized by the Euler angles $\theta=\theta(t)$ and $\phi=\phi(t)$ as

[TABLE]

Such a change of frame is implemented by the unitary operator $V=V(\theta,\phi)$, corresponding to a (time-dependent) global rotation defined by its action:

[TABLE]

which can be written as $V=\exp\left(-i\phi\,s\sum_{j}\sigma_{j}^{z}\right)\,\exp\left(-i\theta\,s\sum_{j}\sigma_{j}^{y}\right)$. The equations of motion of the components $\sigma_{j}^{\alpha}$ of the spins, with $\alpha\in\{X,Y,Z\}$, in the mobile frame $\mathcal{R}$ read

[TABLE]

with the modified Hamiltonian

[TABLE]

which includes an additional $V$-dependent term, in analogy with the emergence of apparent forces in classical mechanics when rotating coordinates are introduced. In our case, a simple calculation shows that

[TABLE]

where $\omega^{X}\equiv\hat{X}\cdot\vec{\omega}=-\sin\theta\;\dot{\phi}$, $\omega^{Y}\equiv\hat{Y}\cdot\vec{\omega}=\dot{\theta}$, and $\omega^{Z}\equiv\hat{Z}\cdot\vec{\omega}=\cos\theta\;\dot{\phi}$. Again, in analogy with classical mechanics, this latter result can be seen as a generalization of Larmor’s theorem: the effect of the rotation of a reference frame is equivalent to the presence of a time-dependent external magnetic field.

The spins on the lattice can now be decomposed on the basis of $\mathcal{R}$ as

[TABLE]

Accordingly, the Hamiltonian Eq. (6), in the mobile frame $\mathcal{R}$, can be written as

[TABLE]

where the last line takes into account the additional term in Eq. (10), while the components of the unit vectors in $\mathcal{R}$ with respect to those of the fixed frame are given in Eq. (7) — for instance, $\hat{X}\cdot\hat{z}=-\sin\theta$, $\hat{Y}\cdot\hat{z}=0$, $\hat{Z}\cdot\hat{z}=\cos\theta$, and analogously for all the other components.

B.2 Time-dependent Holstein-Primakoff transformation

In the rotating frame $\mathcal{R}$ with unit vectors $(\hat{X},\hat{Y},\hat{Z})$, we introduce the spin-wave canonical variables via the Holstein–Primakoff transformation

[TABLE]

where $q_{i}$ and $p_{i}$ are the conjugate canonical variables representing small deviations of the spin away from the $\hat{Z}$-axis, and along the directions $\hat{X}$ and $\hat{Y}$, respectively, while $b_{j}=(q_{j}+ip_{j})/\sqrt{2}$ in our notation. The formal expansion of the operators $\sigma_{i}^{X,Y,Z}$ is in powers of $q_{j}/\sqrt{s}$ and $p_{j}/\sqrt{s}$ and in Eqs. (14) we have retained the leading orders. Accordingly, defining the Fourier space coordinates $\tilde{q}_{k}=N^{-1/2}\sum_{j}e^{-ikj}q_{j}$ and $\tilde{p}_{k}=N^{-1/2}\sum_{j}e^{-ikj}p_{j}$, we get

[TABLE]

At the lowest non-trivial order in the density of spin-waves — controlled by $\epsilon(t)$ introduced below, see Eq. (22) — a straightforward calculation shows that the modulus $|\vec{\Sigma}_{\text{tot}}|$ of the total spin

[TABLE]

is given by

[TABLE]

where $\tilde{b}_{k}\equiv N^{-1/2}\sum_{j}e^{-ikj}b_{j}$. Note that all the excitations with momenta $k$ decrease the total spin projection along the instantaneous direction $\hat{Z}$ of the “vacuum” $\langle\vec{\Sigma}\rangle$ (from Eq. (15), one gets $\Sigma_{\text{tot}}^{Z}\equiv s\tilde{\sigma}_{0}^{Z}=Ns-\sum_{k}\tilde{b}_{k}^{\dagger}\tilde{b}_{k}$), but only excitations with $k\neq 0$ decrease the modulus \big{\lvert}\vec{\Sigma}_{\text{tot}}\big{\rvert} of the total spin, see Eq. (17). In other words, the spin-wave operators with $k=0$ dictate the motion of the spin-wave vacuum.

In order to derive the equations of motion for the spins, one should substitute the expansions (15) into the Hamiltonian (13) and calculate its commutator with the canonical variables $\tilde{q}_{k}$, $\tilde{p}_{k}$. Truncating these expansions at the lowest orders is justified as long as:

The time-dependent angles $\theta(t)$ and $\phi(t)$ which control the rotating frame $\mathcal{R}$ are chosen such that the total spin $\vec{\Sigma}$ remains constantly aligned with the rotating $\hat{Z}$ axis, i.e., $\forall t>0$,

[TABLE] 2. 2.

The spin-waves population remains small, i.e.,

[TABLE]

where we defined

[TABLE]

The first condition is fulfilled by requiring the equations of motion for $S_{X}(t)$ and $S_{Y}(t)$ to be trivially

[TABLE]

These two equations determine the motion of the rotating frame, i.e., the time-evolution of the Euler angles $\theta(t)$ and $\phi(t)$. The validity of the second condition, i.e., of Eq. (19), can be checked by monitoring the time evolution of the total spin-wave density

[TABLE]

which also quantifies to the deviation of the total spin \big{\lvert}\vec{\Sigma}_{\text{tot}}\big{\rvert} from its maximal value $Ns$ according to Eq. (17). The approximations introduced above are consistent as long as $\epsilon(t)\ll 1$; if the total spin-wave density happens to become larger, then higher-order terms in the canonical spin-wave coordinates are expected to contribute to the dynamics at longer times.

B.3 Equations of motion within the Gaussian approximation

The simplest non-trivial approximation consists in treating quantum fluctuations as being harmonic, i.e., within the Gaussian approximation; this means that the expansion in Eq. (15) is substituted into the Hamiltonian (13) and only the linear terms in the vacuum coordinates $(\tilde{q}_{0},\tilde{p}_{0}$) and the quadratic terms in the spatial fluctuation coordinates $(\tilde{q}_{k},\tilde{p}_{k})$ with $k\neq 0$ are kept.

Let us first discuss the mean-field case $J=0$ (discussed in Ref. [21] of the Letter) in which the Hamiltonian $H_{J=0}$ is a function of the total spin $\tilde{\vec{\sigma}}_{0}$ only. Accordingly, the modes $\tilde{q}_{k}$ and $\tilde{p}_{k}$ with $k\neq 0$ enter the Hamiltonian within the Gaussian approximation only via

[TABLE]

i.e., $\langle\tilde{\sigma}^{Z}_{0}\rangle=N(1-\epsilon)$, while each spin-wave number $n_{k}$ (see Eq. (20)) is a constant of motion. In fact, $[n_{k}\,,\,\tilde{\vec{\sigma}}_{0}]=0$ implies $[n_{k}\,,\,H_{J=0}]=0$ and therefore, from Eq. (22),

[TABLE]

which corresponds to the conservation of the total spin (see Eq. (17)). Imposing the condition (18), we get the classical evolution equation of the total spin in the present mean-field case $J=0$, which coincides with the equations of motion found in Ref. [21] of the Letter.

Consider, now, the case $J\neq 0$ which introduces the short-range interaction term $U$ (corresponding to the second line of Eq. (13)) in the infinite-range Hamiltonian discussed above; in particular, our main goal consists in understanding the influence of $U$ on the mean-field dynamics discussed above.In this respect it is convenient to write $U$ as $U=U_{1}+U_{2}+U_{3}$ where

[TABLE]

Expanding $U_{1,2,3}$ by means of Eqs. (15), it is easy to realise that $U_{1}$ gives rise to quadratic terms in $(\tilde{q}_{k},\tilde{p}_{k})$; $U_{3}$ gives rise to a quartic (i.e., two-body) interaction among spin-waves, which is therefore negligible in the low-density limit $\epsilon\ll 1$; $U_{2}$, instead, generates contributions which are simultaneously linear in the vacuum coordinates $(\tilde{q}_{0},\tilde{p}_{0})$ and quadratic in the spin-waves modes $(\tilde{q}_{k},\tilde{p}_{k})$. These terms therefore couple the motion of the vacuum with the spin-waves motion at the lowest non-trivial order: accounting for them is crucial in order to understand the modifications (which we refer to as “feedback”) to the mean-field motion caused by the quantum fluctuations, which is the goal of our work.

The equations of motion of the mobile frame (i.e., by construction, of the collective spin), including the feedback of quantum fluctuations due to $U$, are found by imposing Eq. (21), and read

[TABLE]

where $\epsilon$ is defined in Eq. (22). The equal-time correlation functions appearing in these equations are

[TABLE]

Using now the equations of motion for the spin-waves coordinates, found by computing their commutators with the Hamiltonian $\tilde{H}$,

[TABLE]

one obtains

[TABLE]

Note that the spin-waves — and therefore their correlators $\Delta_{k}^{\alpha\beta}$ with $\alpha$, $\beta\in\{p,q\}$ — have a dynamics even for $J=0$ and $\lambda\neq 0$, though trivial, as it amounts at conserving the spin-wave occupation numbers $n_{k}$ in Eq. (20). These, numbers, instead, are not conserved for $J\neq 0$. In addition, the equations of motion written above are not actually independent because a Gaussian wavefunction such as the one of the spin-waves within the present harmonic approximation is completely specified by two parameters rather than the three $\Delta^{qq}$, $\Delta^{qp}$, and $\Delta^{pp}$. In fact, the latter quantities are actually related by the condition

[TABLE]

which is satisfied at all times and for all values of $k$.

The “feedback” terms $\propto J$ appearing in the equations of motion (28) of the vacuum are of the form

[TABLE]

hence Eq.s (28) can be written as

[TABLE]

where, from Eqs. (22) and (29),

[TABLE]

(cf. Eq. (3) in the Letter). Equations (34) and (31) provide the final system of $2N$ coupled ordinary differential equations which yield the post-quench dynamics at linear order in the spin-wave density $\epsilon$, where we recall that $N$ is the number of spins on the lattice, i.e., the number of possible discrete momenta $k$. These equations are expected not to provide accurate results whenever $\epsilon(t)$ increases and approaches values of order $1$. Note that, as $J\to 0$, the motion of $\theta$ and $\phi$ decouples from the quantum fluctuations and we retrieve the mean-field limit; the same happens of course in the formal classical limit $s\to\infty$.

In order to solve simultaneously the evolution equations (31) and (34), we need to prescribe the initial conditions, which depend on the specific quench under consideration. For quenches of the Hamiltonian $H$ in Eq. (6) originating from the ground state corresponding to $g=g_{0}$ and $J=0$, the initial state is the perfectly coherent state with all the spins pointing in the direction given by the minimum of the classical pre-quench Hamiltonian. Accordingly, the initial conditions turn out to be

[TABLE]

for all $k\neq 0$; in particular, $\epsilon(t=0)=0$. In the Letter we always consider quenches with $g_{0}=0$, for which the pre-quench value of $J$ is actually inconsequential and it may be taken equal to the post-quench value, so that the only parameter affected by the quench is $g$. As a final remark, we note that our approach can be easily generalized in order to deal with a more general class of spin models on an arbitrary graph structure (including higher spatial dimensions) and with arbitrary spin-spin couplings.

Appendix C 2. Perturbation theory at the dynamical critical point: evaluation of $\delta^{qq}(t)$

In this section we discuss the form of the quantum feedback terms $\delta^{\alpha\beta}$ (see Eqs. (33) and (29)) and their long-time behavior in a regime which is analytically tractable. In particular, we consider the unperturbed dynamics with $J=0$ and at the dynamical critical point, as a reference for introducing a leading-order perturbation theory in $J\ll g$, $\lambda$.

The corresponding mean-field dynamics of the angles $\theta(t)$ and $\phi(t)$ reads $\tan\phi(t)=-\tanh(\bar{\lambda}t)$ and $\cos\theta(t)=\tanh^{2}(\bar{\lambda}t)$ (see, for instance, Ref. [21] in the Letter). Inserting these expressions into the system (31) and considering the leading-order contributions in an expansion at long times with $\bar{\lambda}t\gg 1$, the evolution equations take the simpler form

[TABLE]

where, for convenience, we have chosen $\Delta^{qq}_{k}\pm\Delta^{pp}_{k}$ instead of $\Delta^{qq}_{k}$ and $\Delta^{pp}_{k}$, as new active variables. Assuming generic initial conditions of the form $(\Delta^{qq}_{k}+\Delta^{pp}_{k})|_{t=0}=a_{1}$, $(\Delta^{qq}_{k}-\Delta^{pp}_{k})|_{t=0}=a_{2}$, $(\Delta^{pq}_{k})|_{t=0}=a_{3}$, we solve the dynamics prescribed by Eq. (37) at the lowest non-trivial order in $J$. Such a solution allows us to calculate explicitly, for instance, the quantum feedback $\delta^{qq}(t)$ (see Eqs. (33) and (29)) in the thermodynamic limit ($\sum_{k\neq 0}...\to N\int^{2\pi}_{0}\frac{dk}{2\pi}...$)

[TABLE]

where $I_{n}(x)$ is the Bessel function of first kind, with index $n$ and argument $x$. In the long-time limit $Jt\gg 1$, one can employ the asymptotic expansion of $I_{n}(x)$ for $x\gg 1$, specifically for $n=0,1$,

[TABLE]

in order to find that, in addition to a constant term, $\delta^{qq}(t)$ decays as $(Jt)^{-1/2}$ modulated by oscillatory terms of the form $\cos(\pi/4\mp 2(\bar{\lambda}\pm J)t)$ and $\sin(\pi/4\pm 2(\bar{\lambda}\mp J)t)$, which result from beats of the two frequencies $2J$ and $2\bar{\lambda}$. Analogous qualitative results hold for $\delta^{pp}(t)$ and $\delta^{qp}(t)$.

Appendix D 3. Convergence of the MPS-TDVP

As discussed in the Letter, we investigated the dynamical behavior of the system under study also for values of the parameters at which the spin-wave approximation discussed in the previous sections is not expected to be accurate. In this case, we used the matrix product state time-dependent variational principle (MPS-TDVP [23, 24]) and we assessed its viability for investigating the dynamics of the average longitudinal magnetization $S_{x}(t)$ at long times by studying its finite-size scaling, i.e., how it changes upon increasing the systems size $N$. In addition, for each value of $N$ investigated here, we also studied the dependence of the numerical result on the bond dimension $D$. For all simulations, we used a fourth-order integrator with time step 0.02 and the MPS-TDVP formulation developed in  [23, 24].

In Fig. 4 we compare the time evolution of the order parameter $S_{x}(t)$, obtained by using two different bond dimensions $D=80$ and $D=128$, three values of the system size $N=100,~{}200,~{}400$, and three post-quench values of $g=0.5$, 0.83 (corresponding to the ferromagnetic phase), and 1.33 (paramagnetic phase), while keeping fixed $J=2/3$ and $\lambda=1/3$. The systematic error of this approach can be estimated as the difference between curves which differ only for the value of $D$ and it turns out to be of the order of 10% at intermediate times; however, we consistently observe faster convergence in the time averaged order parameter $\bar{S}_{x}$, for which the estimated error is few percents. Interestingly, we observe faster convergence with the bond dimension for larger systems.

The behavior of the dynamics of the order parameter $S_{x}$ in the limit $N\to\infty$ is inferred here from observing how it changes upon increasing $N$, see Fig. 4. While at intermediate times the finite-size effects are still large, the long-time averaged observables such as $\bar{S}_{x}$ converge much faster to the thermodynamic limit, as no changes are observed upon further increasing $N$. Moreover, in the ferromagnetic regime we observe that the order parameter increases upon increasing the system size $N$, indicating a non-vanishing order parameter as $N\rightarrow\infty$. This last observation is valid also in the case where the final magnetization is reversed with respect to the initial state.

Similar fast convergence of $\bar{S}_{x}$ is observed upon increasing the bond dimension $D$ and the system size $N$, also in the chaotic region with initial and final magnetizations of opposite sign, as shown by the curves in Fig. 5. In fact, while they display a significant dependence on $N$ as a consequence of the chaotic behavior consistently observed in that region of the parameter space, the time averaged order parameter $\bar{S}_{x}$ is eventually independent of the bond dimension $D$.