Exact relaxation dynamics and quantum information scrambling in multiply quenched harmonic chains

TL;DR

This paper provides exact analytical results on the relaxation and information scrambling in multiply quenched harmonic chains, revealing local thermalization despite integrability.

Contribution

It presents the first exact analysis of relaxation dynamics and quantum information scrambling in multiply quenched harmonic chains, a trivially integrable system.

Findings

Out of time correlator indicates scrambling at large times.

Entanglement entropy saturates to a steady state.

Reduced density matrix becomes diagonal long after quenches.

Abstract

The quantum dynamics of isolated systems under quench condition exhibits a variety of interesting features depending on the integrable/chaotic nature of system. We study the exact dynamics of trivially integrable system of harmonic chains under a multiple quench protocol. Out of time ordered correlator of two Hermitian operators at large time displays scrambling in the thermodynamic limit. In this limit, the entanglement entropy and the central component of momentum distribution both saturate to a steady state value. We also show that reduced density matrix assumes the diagonal form long after multiple quenches for large system size. These exact results involving infinite dimensional Hilbert space are indicative of local thermal behaviour for a trivially integrable harmonic chain.

Figures

Equations

$$H(t)$$

$$=\frac{1}{2}\left[\sum _{j=1}^N{p}_j^2+X^T.\Sigma (t).X\right],$$

$$i\frac{\partial }{\partial t}\mid \psi ⟩=H(t)\mid \psi ⟩$$

$$\psi (x_1,....,x_N,t)=$$

$$\times \exp \left[-\frac{X^T\Omega X}{2}\right],$$

$${\ddot{b}}_j+{\lambda }_j(t)b_j=\frac{{\lambda }_j(0)}{b_j^3},$$

$${\lambda }_j(t)={\omega }^2(t)+2k(t)-2k(t)\cos (2\pi j/N)$$

$${\eta }_j(t)\equiv b_j^2(t),$$

$${\ddot{\eta }}_j{\eta }_j-\frac{1}{2}{\dot{\eta }}_j^2+2{\lambda }_j(t){\eta }_j^2=2{\lambda }_j(0).$$

$${\eta }_{j,1}(t)={\alpha }_{j,1}\cos (2\sqrt{{\lambda }_j(t)}t)+{\beta }_{j,1}\sin (2\sqrt{{\lambda }_j(t)}t)+{\gamma }_{j,1}$$

$$0<t<T$$

$${\eta }_{j,2}(t)={\alpha }_{j,2}\cos (2\sqrt{{\lambda }_j(t)}t)+{\beta }_{j,2}\sin (2\sqrt{{\lambda }_j(t)}t)+{\gamma }_{j,2}$$

$$T<t<2T$$

$${\eta }_{j,i}(t)={\alpha }_{j,i}\cos 2(\sqrt{{\lambda }_j(t)}t)+{\beta }_{j,i}\sin (2\sqrt{{\lambda }_j(t)}t)+{\gamma }_{j,i}$$

$$(i-1)T<t<iT$$

$$⋮$$

$${\eta }_{j,n}(t)={\alpha }_{j,n}\cos 2(\sqrt{{\lambda }_j(t)}t)+{\beta }_{j,n}\sin (2\sqrt{{\lambda }_j(t)}t)+{\gamma }_{j,n}$$

$$(n-1)T<t<\infty$$

$${\vec{\eta }}_j(t)=\left(\begin{array}{c}{\eta }_j(t)\\ {\dot{\eta }}_j(t)\\ {\ddot{\eta }}_j(t)\end{array}\right)=A_j(t)\left(\begin{array}{c}{\alpha }_{j,i}\\ {\beta }_{j,i}\\ {\gamma }_{j,i}\end{array}\right),$$

$$A_j(t)=$$

$$\left(\begin{array}{ccc}\cos (2\sqrt{{\lambda }_j(t)}t)& \sin (2\sqrt{{\lambda }_j(t)}t)& 1\\ -2\sqrt{{\lambda }_j(t)}\sin (2\sqrt{{\lambda }_j(t)}t)& 2\sqrt{{\lambda }_j(t)}\cos (2\sqrt{{\lambda }_j(t)}t)& 0\\ -4{\lambda }_j(t)\cos (2\sqrt{{\lambda }_j(t)}t)& -4{\lambda }_j(t)\sin (2\sqrt{{\lambda }_j(t)}t)& 0\end{array}\right).$$

$${\vec{\eta }}_j(t^{+})=B_j(t^{+}\leftarrow t^{-}){\vec{\eta }}_j(t^{-}),$$

$$B_j=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 2({\lambda }_j(t^{-})-{\lambda }_j(t^{+}))& 0& 1\end{array}\right).$$

$$\left(\begin{array}{c}{\alpha }_{j,i+1}\\ {\beta }_{j,i+1}\\ {\gamma }_{j,i+1}\end{array}\right)=A^{-1}(t^{+})B(t^{+}\leftarrow t^{-})A(t^{-})\left(\begin{array}{c}{\alpha }_{j,i}\\ {\beta }_{j,i}\\ {\gamma }_{j,i}\end{array}\right),$$

$$\begin{array}{r}\left(\begin{array}{c}{\alpha }_{j,1}\\ {\beta }_{j,1}\\ {\gamma }_{j,1}\end{array}\right)=\left(\begin{array}{c}\frac{{\lambda }_{j,f}-{\lambda }_{j,in}}{2{\lambda }_{j,f}}\\ 0\\ \frac{{\lambda }_{j,f}+{\lambda }_{j,in}}{2{\lambda }_{j,f}}\end{array}\right)\end{array}.$$

$${\rho }_{red}(x_i,x_i^{\prime },t)=\int \prod _{a=1\ne i}^{N}d{X}^{a}d{X}^{\prime a}{\psi }^{\ast }(X^a,t)\psi (X^{\prime a},t).$$

$${\rho }_{red}(x_i,x_i^{\prime },t)=$$

$$\times \exp \left[-\frac{1}{2}\gamma \left(x_i^2+x_i^{\prime 2}\right)+x_i\beta x_i^{\prime }\right],$$

$$Z(t)={\tilde{b}}_{1\times 1}-{\tilde{b}}_{N-1\times 1}^{†}{\Omega }_{N-1\times N-1}^{-1}{\Omega }_{N-1\times 1},$$

$$\gamma (t)={\Omega }_{1\times 1}-\frac{1}{2}{\Omega }_{N-1\times 1}^{†}{\Omega }_{N-1\times N-1}^{-1}{\Omega }_{N-1\times 1}$$

$$+2{\tilde{b}}_{N-1\times 1}^{†}{\Omega }_{N-1\times N-1}^{-1}{\tilde{b}}_{N-1\times 1},$$

$$\beta (t)=\frac{1}{2}{\Omega }_{N-1\times 1}^{†}{\Omega }_{N-1\times N-1}^{-1}{\Omega }_{N-1\times 1}$$

$$+2{\tilde{b}}_{N-1\times 1}^{†}{\Omega }_{N-1\times N-1}^{-1}{\tilde{b}}_{N-1\times 1}.$$

$${\rho }_{red}(x_i,x_i+\Delta ,t)={\left(\frac{\gamma -\beta }{\pi }\right)}^{\frac{1}{2}}\exp \left[-iZ\left(\Delta x_i+{\Delta }^2\right)\right]$$

$$\times \exp \left[-(\gamma -\beta ){\left(x_i-\frac{\Delta }{2}\right)}^2-{\Delta }^2\left(\gamma -\frac{\gamma -\beta }{4}\right)\right],$$

$$r$$

$$=\exp \left[-iZ\left(\Delta x_i+{\Delta }^2\right)\right]\exp \left[-(\gamma -\beta )x_i\Delta -\gamma {\Delta }^2\right].$$

$$n(p,t)$$

$$=\int d{x}_{i}d{x}_i^{\prime }{\left(\frac{\gamma -\beta }{4{\pi }^3}\right)}^{\frac{1}{2}}{e}^{[iZ(x_i^2-x_i^{\prime 2})-ip(x_i-x_i^{\prime })]}$$

$$\times e^{[-\frac{\gamma }{2}(x_i^2+x_i^{\prime 2})+\beta x_i{x}_i^{\prime }]}.$$

Full text

I Introduction

The behaviour of isolated quantum systems under non-equilibrium conditions is of great interest from both theoretical Rigol et al. (2007); Marcos et al. (2008); Polkovnikov et al. (2011); Srednicki (1994); D’Alessio et al. (2016); Deutsch (2018) as well as experimental point of view Islam et al. (2015); Vidmar et al. (2015); Kaufman et al. (2016). A generic method to reach a non-equilibrium regime is by a single or multiple quenches of the system parameters. For non-integrable systems, the quench protocol often leads to thermalization Srednicki (1994); Marcos et al. (2008), although there can be subtleties depending on the choice of initial states Barthel and Schollwöck (2008); Cramer et al. (2008a).

For integrable systems, the situation is not as unequivocal. It is generally known that an isolated integrable system does not thermalize Marcos et al. (2008); Deutsch (2018), but can be described by a generalized Gibbs ensemble (GGE) Rigol et al. (2007); Vidmar and Rigol (2016); Essler and Fagotti (2016); Fioretto and Mussardo (2010); Wright et al. (2014).However, it has recently been found that the $XXX$ Heisenberg spin chain, integrable under the Bethe ansatz scheme, can exhibit weak eigenstate thermalization for a typical, but not necessarily all eigenstates Alba (2015); Alba and Calabrese (2017); Grisins and Mazets (2011). It is also known that certain finite dimensional systems exhibit comparable statistical relaxation regardless of whether they are integrable or not Santos et al. (2012). Furthermore, an integrable Jaynes-Cummings model can show rapid decoherence analogous to chaotic dynamics and the inability to recover the purity at large times increases in the question of dynamical relaxation thermodynamic limit Angelo et al. (1999). For an integrable system with infinite dimensional Hilbert space, Bogoliubov and Krylov showed that under certain assumption about the thermal reservoir, the system relaxes to thermal state in thermodynamic limit (see Chirikov (1986) and references therein). In 1d-bosonic systems, a class of special initial states thermalizes under integrable dynamics governed by Gross-Pitaevskii equation Grisins and Mazets (2011).

Various physical quantities have been used to analyze the onset of statistical relaxation following a quench. In a quantum many-body system, the reduced density matrix (RDM) carries important information regarding the relaxation dynamics Peschel and Chung (1999); Chung and Peschel (2000); Peschel (2004); Peschel and Eisler (2009); Peschel (2003); Fagotti and Essler (2013); Calabrese et al. (2012). The entanglement entropy calculated from the RDM measures the loss of information and its dynamics describes how the quantum information is spread. The one body momentum distribution, also obtained from the RDM, carries the signature of the dynamical relaxation.

Recently, the out of time order correlator (OTOC) Larkin and Ovchinnikov (1969) has gained prominence in the context of scrambling of quantum information in non-equilibrium systems Maldacena et al. (2016); Rozenbaum et al. (2017); Swingle (2018); Heyl et al. (2018); Lin and Motrunich (2018a, b); Li et al. (2017); Kukuljan et al. (2017); Luitz and Bar Lev (2017); Cotler et al. (2018). Although information scrambling is usually a property of chaotic systems, the OTOC of certain non-local operators exhibit scrambling even in an integrable Ising chain Lin and Motrunich (2018a). It has been further argued that scrambling could be independent of the integrability of the Hamiltonian Iyoda and Sagawa (2018) and that mutual information in an integrable spin chain can exhibit weak scrambling Alba and Calabrese (2019). The above list of examples, which is by no means exhaustive, indicates that there are still many open questions in the context of dynamical relaxation and information scrambling in integrable systems.

In this paper we discuss two different dynamical aspects of an integrable harmonic chain. The first point is to address the question of dynamical relaxation to a steady state under a multiple quench protocol. The second point is to analyze the dynamics of information scrambling under a similar set of quenches.

The main differences of our analysis with the existing approaches in the literature are

• •

We consider a series of multiple global quenches with no restriction to their number, in contrast to a single quench which is typical in the existing literature. As we shall show, the multiple quench protocol leads to important differences compared to a single quench.

• •

Our analysis involves systems with infinite dimensional Hilbert space at each site and local Hermitian operators. This is particularly relevant for information scrambling and OTOC where majority of the existing analyses involves finite spin systems and unitary operators.

• •

The time development is treated in an exact analytical fashion by solving TDSE using non-linear Ermakov equation. The total number of particles and the number of quenches can be arbitrarily large and the analysis is valid for indefinitely large time.

In order to achieve the above, we study the exact relaxation dynamics following multiple global quenches in a harmonic chain with $N$ oscillators, which is trivially integrable. For this system, the exact entanglement dynamics following a single quench has been studied earlier Ghosh et al. (2017). Under a multiple quench protocol, we show that the off-diagonal terms in the one-particle RDM Peschel and Chung (1999); Peschel (2003) vanish exponentially fast in time in the limit of large $N$, while the diagonal terms saturate to a steady value. The mixing of a large number of incommensurate and irrational normal mode frequencies is responsible for this feature. In contrast, for a single quench, the relaxation to the steady state is much slower. Under the same conditions, the momentum distribution and the entanglement entropy also show the signatures of a steady state.

Universality of the small time behaviour of OTOC in chaotic systems and its relation with Lyapunov exponents has been the focus of most of the studies Rozenbaum et al. (2017); Maldacena et al. (2016); Swingle (2018); Lakshminarayan and Prakash (2019). Here we want to focus on long time behaviour much after the completion of the quench protocol. The multiple quench protocol has important consequences in the quantum information scrambling and OTOC. While recurrences are characteristic feature of integrability in OTOC, in contrast for this system, OTOC saturates to a nonzero value under multiple quenches in thermodynamic limit, which is indicative of scrambling Alba and Calabrese (2019). As our exact analysis is valid for arbitrarily large number of particles, we can clearly demonstrate the finite size effect by varying the number of oscillators.

The harmonic chain under consideration can be experimentally realized in various systems such as in optical tweezers Barredo et al. (2016) and the tuning of individual coupling parameters can be done using ultracold atoms Hunger et al. (2011) or Rydberg states Buchmann et al. (2017); Macrì and Pohl (2014). In particular, it can be simulated using Bose-Hubbard model in the strong superfluid phaseCramer et al. (2008b). Various quench protocols have already been experimentally realized for the Bose-Hubbard model using cold atoms in optical lattices Vidmar et al. (2015); Kaufman et al. (2016); Islam et al. (2015). It is thus plausible that the predictions of our work can be empirically verified.

This paper is organized as follow. In Sec. 2, we setup the formalism for multiple quench protocol and obtain the solutions of TDSE using non-linear Ermakov equations. In Sec. 3, we discuss the time dependence of RDM and show that the off-diagonal elements vanish for large time and in the thermodynamic limit. In Sec. 4, we obtain the momentum distribution and entanglement entropy using the RDM and discuss their steady state properties. Quantum information scrambling is discussed in Sec. 5 and it is shown that OTOC saturates to a non-zero steady value for the case of multiple quenches in the same limit. We conclude the paper in Sec. 6 with a summary of results and remarks.

II Harmonic chain and the quench protocol

We consider an isolated harmonic oscillator chain with $N$ oscillators and with periodic boundary condition. The Hamiltonian is given by

[TABLE]

where $X$=$(x_{1},x_{2},...,x_{N})^{T}$ and $\Sigma$ is an $N\times N$ real symmetric matrix. Here the frequency $\omega$ and the nearest neighbour coupling $k$ are explicit functions of time. A solution of the corresponding time dependent Schrödinger’s equation (TDSE) Lewis and Riesenfeld (1969); Lohe (2009); Ghosh et al. (2017)

[TABLE]

can be written as

[TABLE]

where $\Omega=U^{T}\sqrt{\Sigma^{\prime D}}U$, $\Sigma_{jj}^{\prime D}=\frac{\Sigma_{jj}^{D}(0)}{b_{j}^{4}(t)}$, $\tilde{b}=U^{T}\tilde{b}^{D}U$,$\tau_{j}=\int_{0}^{t}\frac{dt^{\prime}}{b_{j}^{2}(t^{\prime})}$, $U$ is an orthogonal transformation which transforms the matrix $\Sigma$ to its diagonalized form $\Sigma^{D}$ and $\tilde{b}^{D}$ is a diagonal matrix with elements $\frac{\dot{b}_{j}(t)}{2b_{j}(t)}$. Here $b_{j}(t)$ satisfies the the nonlinear Ermakov equation Lewis and Riesenfeld (1969); Lohe (2009); Pinney (1950) given by

[TABLE]

where $\lambda_{j}$’s are the normal mode frequencies of the Hamiltonian, which have the form

[TABLE]

with $\lambda_{j}=\lambda_{N-j}$. We choose the initial condition as the ground state of the time independent pre-quenched Hamiltonian of the $N$ oscillator chain, which requires that $b(0)=1$ and $\dot{b}(0)=0$. Note that finding the solutions of the TDSE is equivalent to finding the solutions of the Ermakov equations (4). The normalization is determined by the condition that $\int_{-\infty}^{\infty}\psi^{*}(x,t)\psi(x,t)dx=1$ Lohe (2009).

Solution for series of quenches: Our quench protocol is shown in Fig. (1). At time $t=0$, the frequency $\omega(i)$ is quenched to $\omega(f)$ and $k(=2)$ is left unchanged. This defines a single quench. After a time $T$, the frequency $\omega(f)$ is changed back to $\omega(i)$, defining the second quench. This sequence is repeated till the required number of quenches is achieved. Now we shall discuss the solution of Ermakov equations under $n$ such quenches. In terms of the variable

[TABLE]

the Ermakov equations (4) can be written as

[TABLE]

The corresponding solutions are given by

[TABLE]

where $\alpha_{j,i}$, $\beta_{j,i}$ and $\gamma_{j,i}$ are the time independent constants which are obtained from the boundary conditions on $b$ or $\eta$. Furthermore, using the above equation we could express $\eta_{j}$ and its derivatives in a matrix form as

[TABLE]

where

[TABLE]

The continuity of the wavefunction across the quenches implies the continuity in $\eta_{j}(t)$ and $\dot{\eta}_{j}(t)$. $\vec{\eta}_{j}(t^{+})$ and $\vec{\eta}_{j}(t^{-})$ defined just after and before any quench respectively are related by

[TABLE]

where the form of B matrix could be found using Eq.(7) as

[TABLE]

The coefficient vectors with elements $(\alpha,\beta,\gamma)_{j}$ across the quench are related as

[TABLE]

with the boundary condition,

[TABLE]

This gives a complete solution of the Ermakov equations for the given quench protocol.

III Reduced density matrix

In this Section we shall discuss the time-dependent one-body RDM obtained under a series of quenches in the coordinate representation. This RDM would subsequently be used to obtain the one-body momentum distribution and the entanglement entropy.

In order to obtain the time dependent one-body RDM, we trace out the rest of the system except the $i^{\text{th}}$ site. The RDM has the form

[TABLE]

Using (3) in (15) and carrying out the integration, we get

[TABLE]

where $Z$, $\beta$, $\gamma$ are given by

[TABLE]

A generic off-diagonal element of the RDM has the form

[TABLE]

where $\Delta$ denotes the distance from the diagonal axis of the RDM and $\Delta=0$ corresponds to the diagonal elements of the RDM. The equation (18) shows that the off diagonal elements are shifted Gaussians with a time dependent exponent. The ratio $r$ between the diagonal and a generic off-diagonal element can be written as

[TABLE]

The exponents $\gamma$ and $(\gamma-\beta)$ are plotted in Fig. 2 and are found to be real positive for all time. This also follows from the reality of the eigenvalues of the RDM, which is Hermitian Ghosh et al. (2017). For a series of multiple quenches we find that $(\gamma-\beta)$ rapidly saturates to a small positive value irrespective of $N$. However $\gamma$ has very different behaviour depending on the value of $N$. For higher $N$, $\gamma$ has a large value which oscillates very little with time. For smaller $N$, the mean value of $\gamma$ is also high but the fluctuations are very big as well. On the other hand, for a single quench, both $\gamma$ and $(\gamma-\beta)$ show appreciable oscillations and the value of $\gamma$ is much less irrespective of the value of $N$.

In Eq.(III) we have given the ratio of the off-diagonal to the diagonal elements of the RDM. We thus find that for multiple quenches and in the thermodynamic limit, the off-diagonal matrix elements of the RDM tend to zero. For a single quench, the suppression of the off-diagonal matrix elements with time is much less.

Hence, for multiple quenches with the large system size, the RDM assumes a diagonal form at a time large compared to the duration of the quenches. This can be qualitatively understood as follows. The solution of the Ermakov equation for each normal mode contains an irrational number given by the square root of corresponding normal mode frequency $\lambda_{j}$. As the number of oscillators increases, a large number of irrational and incommensurate frequencies start contributing to the wave function and the RDM. The mixing of a large number of such modes is responsible for the statistical relaxation of the quantities $\gamma$ and $(\gamma-\beta)$ with time, which in turn reduces the RDM to a diagonal form.

IV momentum distribution and entropy

We derive the analytical expression of one body momentum distribution by taking the Fourier transform of (16), which is given by

[TABLE]

After the integration, the momentum distribution takes the form

[TABLE]

At time $t=0$ just before the quench, the distribution (21) is a Gaussian with mean at $p=0$. In Fig. 3 the central component of the momentum distribution is plotted with time for both single and multiple quenches and with $N=500$. For a single quench, in the long time limit the momentum distributions reaches a steady value but still shows appreciable oscillations. For multiple quenches the value in the long time limit is lower and the fluctuations are negligible. Thus in the long time and in the thermodynamic limit, the higher the number of quenches, the better is the relaxation of the system to a steady state.

In order to study the entanglement entropy, we bipartite the system into two parts of one oscillator versus $N-1$ oscillators. The entanglement entropy of the smaller subsystem calculated using the reduced density matrix (16) has the form Ghosh et al. (2017)

[TABLE]

where the $\xi$ has the form given by

[TABLE]

$\beta$, $\gamma$ are given in eq. (LABEL:beta-gamma).

The von Neumann entropy as a function of time is plotted in Fig. 4. In a finite dimensional system, the maximum entropy is proportional to the Hilbert space dimension. The bound on entanglement entropy is saturated for generic chaotic systems. Here the system being infinite dimensional, there is no finite limit to the maximum entanglement entropy. As in the case of the momentum distribution, higher number of quenches leads to smaller fluctuations in the entanglement entropy and a larger steady state value. This can be qualitatively understood in terms of reduced variation in eigenvalues of RDM given by $p_{n}=(1-\xi)\xi^{n},n=0,1,2,\dots$ Ghosh et al. (2017). The variation in $p_{n}$ reduces as $\xi\to 1$ or in other words as $(\gamma-\beta)\to 0$ (see Fig. 2).

V Quantum information scrambling

The OTOC between two operators $M$ and $N$ separated by a lattice distance $l$ is defined as

[TABLE]

where $l=j-i$. For well separated local operators, the OTOC starts from a zero value and then increases as the information propagates with Lieb-Robinson velocity.

Here we consider the OTOC between position and momentum operators, which are local and Hermitian. We choose $M_{i}(t)=x_{i}(t)$, $N_{i}(0)=p_{i}(0)$ which are labeled by the site index. The expression for $x_{i}(t)$ in the Heisenberg picture for a single time dependent quantum harmonic oscillator with frequency $\omega(t)$ is given by Ji et al. (1995)

[TABLE]

where $b(t)$ is the solution of corresponding Ermakov equation. Using Eq.(V) with the canonical commutation relations $[x_{i},p_{j}]=i\hbar\delta_{ij}$ we get,

[TABLE]

In Fig.(5), OTOC is plotted for a fixed system size $N=500$ and with different number of quenches. We observe that OTOC saturates to a non-zero value with very small fluctuation for the case of multiple quenches in contrast to a constant slow decrease with large fluctuation for single quench. A strong system size dependence on the OTOC is clearly seen from Fig. 6. We find that even with multiple quenches, the OTOC fluctuates enormously for smaller number of particles and does not saturate to any steady value. In Fig.(5), the quasi-recurrence in the case of multiple quenches for $N=500$ occurs at a time $t\sim\frac{N}{v_{max}}$ where $v_{max}$ is the Lieb-Robinson velocity for this system with the given set of parameters after quench. We also note that the quasi-recurrences decrease in the thermodynamic limit as well as with the higher number of quenches.

The saturation of the OTOC to a non-zero steady value is indicative of information scrambling. Even though the system under consideration is integrable, we find that under multiple quenches and in the thermodynamic limit, OTOC saturates to a non-zero steady value. In this formalism, the saturation can be attributed to the mixing of large number of modes with incommensurate frequencies. In the conventional picture of expanding the state as superposition of eigenfunctions of post-quench Hamiltonian, the multiple quenches would amount to accessing larger proportions of Hilbert space. This is also consistent with larger entanglement entropy with increasing number of quenches.

VI conclusion

In this paper we have analyzed the relaxation dynamics and quantum information scrambling in an isolated harmonic chain under multiple quenches. The various physical quantities show remarkably different non-equilibrium behaviour under the multiple quench protocol compared to a single quench. The RDM has been shown to assume a diagonal form exponentially fast compared to a single quench. The entanglement dynamics and the momentum distribution also show relaxation to a steady state. The exact analytical results obtained here are valid for arbitrary number of particles. In order to demonstrate the finite size effects, we have graphically exhibited our results for $N=20$ and $N=500$. It is clearly seen that the quasi-revivals of the physical observables, characteristic of the finite size effects, reduce remarkably as $N$ is increased.

The conclusions are similar for the quantum information scrambling and OTOC. For a single quench and for a low value of $N=20$, the OTOC shows almost complete revival and the scrambling is practically non-existent. On the other hand, for five quenches and with $N=500$, the OTOC saturates to a finite steady value with very little revival. It may be noted that we have considered Hermitian operators to evaluate the OTOC and the system has an infinite dimensional Hilbert space. This is very different compared to the usual finite spin systems where the OTOC is normally evaluated using unitary operators.

The above result indicates that for multiple quenches and in the thermodynamic limit, the integrable harmonic chain relaxes to a steady state and the OTOC exhibits quantum information scrambling. In our formalism this has been achieved using an exact solution of the TDSE valid throughout the quench protocol and all the observables have been evaluated using the exact solution. As the number $N$ of the oscillators is increased, a large number of irrational and incommensurate normal mode frequencies start contributing to the observables whose number is of the order of $N$. The mixing of these incommensurate frequencies leads to the steady state and saturating behaviour of the various physical quantities.

One question that has not been addressed here is what is the nature of the equilibrated state. For an integrable harmonic chain, the emergence of GGE has already been established Biroli et al. (2010). Whether this happens for the system under consideration with multiple quenches remains to be seen. A related question is the emergence of weak eigenstate thermalization, which has recently been observed for certain integrable systems Alba (2015); Alba and Calabrese (2017); Grisins and Mazets (2011). In this paper, we have focused only on a single state of the TDSE whereas to address the questions on thermalization, a more complete knowledge of the time dependent many-body spectrum would be needed. This is beyond the scope of the present analysis, which we hope to address these in future work.