Return amplitude after a quantum quench in the XY chain

TL;DR

This paper derives an exact formula for transition amplitudes in the quantum XY chain, analyzes return amplitudes for specific states, and explores finite-size effects and dynamical phase transitions.

Contribution

It provides a new exact analytical expression for transition amplitudes and investigates finite-size effects and dynamical phase transitions in the XY chain.

Findings

Exact formula for transition amplitudes between eigenstates.

Analytical expression for return amplitude of polarized and Nél states.

Finite-size effects and traversal phenomena depend on initial state and system size.

Abstract

We determine an exact formula for the transition amplitude between any two arbitrary eigenstates of the local $z$-magnetization operators in the quantum XY chain. We further use this formula to obtain an analytical expression for the return amplitude of fully polarized states and the N\'eel state on a ring of length $L$. Then, we investigate finite-size effects in the return amplitude: in particular quasi-particle interference halfway along the ring, a phenomenon that has been dubbed traversal~\cite{FE2016}. We show that the traversal time and the features of the return amplitude at the traversal time depend on the initial state and on the parity of $L$. Finally, we briefly discuss non-analyticities in time of the decay rates in the thermodynamic limit $L\rightarrow\infty$, which are known as dynamical phase transitions.

Figures

Equations

$$F_{\Psi }(t)=|⟨\Psi |e^{-i\mathcal{H}t}|\Psi ⟩|,$$

$${\bar{\Gamma }}_{\Psi }\equiv \underset{t\to \infty }{\lim }\underset{L\to \infty }{\lim }-\frac{1}{L}\log F_{\Psi }$$

$${\mathcal{H}}_{\text{free}}={\text{c}}^{†}\text{A}\text{c}+\frac{1}{2}{\text{c}}^{†}\text{B}{\text{c}}^{†}+\frac{1}{2}\text{c}{\text{B}}^T\text{c}-\frac{1}{2}\mathrm{Tr}\text{A},$$

$$e^{i{\mathcal{H}}_{\text{free}}t}\mathbf{\gamma }{e}^{-i{\mathcal{H}}_{\text{free}}t}=\text{T}\mathbf{\gamma },$$

$$\text{T}=e^{-it\left(\begin{array}{cc}\text{A}& \text{B}\\ \text{-B}& \text{-A}\end{array}\right)}\equiv \left(\begin{array}{cc}{\text{T}}_{11}& {\text{T}}_{12}\\ {\text{T}}_{21}& {\text{T}}_{22}\end{array}\right).$$

$$F_{\Psi }(t)=\left|⟨\Psi |e^{-i{\mathcal{H}}_{\text{free}}t}|\Psi ⟩\right|.$$

$$e^{-i{\mathcal{H}}_{\text{free}}t}=e^{-\frac{1}{2}\mathrm{Tr}\text{Y}}{e}^{\frac{1}{2}{\text{c}}^{†}\text{X}{\text{c}}^{†}}{e}^{{\text{c}}^{†}\text{Y}\text{c}}{e}^{\frac{1}{2}\text{c}\text{Z}\text{c}},$$

$$\text{X}={\text{T}}_{12}[{\text{T}}_{22}{]}^{-1},\text{Z}=[{\text{T}}_{22}{]}^{-1}{\text{T}}_{21},e^{\text{-Y}}=[{\text{T}}_{22}{]}^T.$$

$$|\mathbf{\xi }⟩=e^{-\mathbf{\xi }\cdot {\mathbf{c}}^{†}}|\Omega ⟩,$$

$$\int d\bar{\mathbf{\xi }}d\mathbf{\xi }{e}^{-\bar{\mathbf{\xi }}\cdot \mathbf{\xi }}|\mathbf{\xi }⟩⟨\bar{\mathbf{\xi }}|=\mathbf{I},⟨\bar{\mathbf{\xi }}|\mathbf{\xi }⟩=e^{\bar{\mathbf{\xi }}\cdot \mathbf{\xi }},$$

$$f_{\Psi }(t)=⟨\Psi |e^{-it{\mathcal{H}}_{\text{free}}}|\Psi ⟩,$$

$$f_{\Psi }(t)=e^{-\frac{1}{2}\text{Tr}(\text{Y})}\int d\bar{\mathbf{\xi }}d\mathbf{\xi }\int d\bar{\mathbf{\eta }}d\mathbf{\eta }{e}^{-\bar{\mathbf{\xi }}\cdot \mathbf{\xi }-\bar{\mathbf{\eta }}\cdot \mathbf{\eta }+\frac{1}{2}\bar{\mathbf{\eta }}\text{X}\bar{\mathbf{\eta }}+\frac{1}{2}\mathbf{\xi }\mathbf{Z}\mathbf{\xi }}⟨\Psi |\mathbf{\eta }⟩⟨\bar{\mathbf{\eta }}|e^{{\mathbf{c}}^{†}\mathbf{Y}\mathbf{c}}|\mathbf{\xi }⟩⟨\bar{\mathbf{\xi }}|\Psi ⟩.$$

$$f_{\Psi }(t)=e^{-\frac{1}{2}\text{Tr}(\mathbf{Y})}\int d\bar{\mathbf{\xi }}d\mathbf{\xi }\int d\bar{\mathbf{\eta }}d\mathbf{\eta }{e}^{-\bar{\mathbf{\xi }}\cdot \mathbf{\xi }-\bar{\mathbf{\eta }}\cdot \mathbf{\eta }+\frac{1}{2}\bar{\mathbf{\eta }}\mathbf{X}\bar{\mathbf{\eta }}+\frac{1}{2}\mathbf{\xi }\mathbf{Z}\mathbf{\xi }+\bar{\mathbf{\eta }}{\mathbf{e}}^{\mathbf{Y}}\mathbf{\xi }}\prod _{j\in S_{\psi }}{\eta }_j\prod _{j\in S_{\psi }}{\bar{\xi }}_j.$$

$$f_{\Psi }(t)$$

$$=e^{-\frac{1}{2}\text{Tr}(\text{Y})}\text{Pf}(\tilde{\text{M}}),$$

$$\text{M}=\left(\begin{array}{cc}\text{X}& e^{\text{Y}}\\ -e^{{\text{Y}}^T}& \text{Z}\end{array}\right),$$

$$F_{\Psi }(t)=\left|e^{-\frac{1}{2}\text{Tr}(\text{Y})}\text{Pf}(\tilde{\text{M}})\right|,$$

$$|\Psi ⟩=\int d{\xi }_{i_m}\dots d{\xi }_{i_1}|{\mathbf{\xi }}_{\Psi }⟩,$$

$$\displaystyle\ \mathcal{H}_{XY}=-\frac{J}{2}\sum_{j=1}^{L}\Big{[}\left(\frac{1+\gamma}{2}\right)\sigma_{j}^{x}\sigma_{j+1}^{x}+\left(\frac{1-\gamma}{2}\right)\sigma_{j}^{y}\sigma_{j+1}^{y}\Big{]}-\frac{h}{2}\sum_{j=1}^{L}\sigma_{j}^{z},$$

$${\mathcal{H}}_{XY}=\frac{J}{2}\sum _{j=1}^L(c_j^{†}{c}_{j+1}+\gamma c_j^{†}{c}_{j+1}^{†}+\text{H.c.})-\frac{h}{2}\sum _{j=1}^L(2c_j^{†}{c}_j-1),$$

$$P\equiv \prod _{j=1}^L{\sigma }_j^z.$$

$$[\text{A}{]}_{ij}=-{\delta }_{i,j}h+\frac{J}{2}({\delta }_{i+1,j}+{\delta }_{i-1,j})-\frac{J\mathcal{N}}{2}({\delta }_{i,1}{\delta }_{j,L}+{\delta }_{i,L}{\delta }_{j,1}),$$

$$[\text{B}{]}_{ij}=\frac{J\gamma }{2}({\delta }_{i+1,j}-{\delta }_{i-1,j})+\frac{J\gamma \mathcal{N}}{2}({\delta }_{i,1}{\delta }_{j,L}-{\delta }_{i,L}{\delta }_{j,1}).$$

$$\mathcal{N}{e}^{iL{\varphi }_k}=-1.$$

$${\lambda }_k^{\text{A}}=-h+\cos ({\varphi }_k);{\lambda }_k^{\text{B}}=i\gamma \sin ({\varphi }_k);$$

$$\text{exp}\left[-it\left(\begin{array}{cc}{\lambda }_k^{\text{A}}& {\lambda }_k^{\text{B}}\\ -{\lambda }_k^{\text{B}}& -{\lambda }_k^{\text{A}}\end{array}\right)\right]=\left(\begin{array}{cc}{\lambda }_k^{{\text{T}}_{11}}& {\lambda }_k^{{\text{T}}_{12}}\\ {\lambda }_k^{{\text{T}}_{21}}& {\lambda }_k^{{\text{T}}_{22}}\end{array}\right).$$

$${\lambda }_k^{\text{X}}=-\frac{i{\lambda }_k^{\text{B}}}{i{\lambda }_k^{\text{A}}+{\epsilon }_k\cot (t{\epsilon }_k)},{\lambda }_k^{\text{Z}}=-{\lambda }_k^{\text{X}},{\lambda }_k^{\text{Q}}=\frac{1}{\cos (t{\epsilon }_k)+i\frac{{\lambda }_k^{\text{A}}}{{\epsilon }_k}\sin (t{\epsilon }_k)},$$

$${ϵ}_k=\sqrt{(\cos {\varphi }_k-h{)}^2+{\gamma }^2{\sin }^2{\varphi }_k}.$$

$$[\text{X}{]}_{lj}=\frac{1}{L}\sum _{k=1}^L{\lambda }_k^{\text{X}}{e}^{i{\varphi }_k(l-j)};[\text{Q}{]}_{lj}=\frac{1}{L}\sum _{k=1}^L{\lambda }_k^{\text{Q}}{e}^{i{\varphi }_k(l-j)}.$$

$${\Gamma }_{\mathcal{B}}(t)=-\frac{1}{L}\log F_{\mathcal{B}}(t),$$

$$F_{\uparrow }(t)=\left|e^{-\frac{1}{2}\text{Tr}(\text{Y})}\right|=\prod _{k=1}^L|{\lambda }_k^{\text{Q}}{|}^{-1/2},$$

Full text

I Introduction

A prominent theoretical tool to investigate non-equilibrium phenomena in many-body systems is the quantum quench. In the original formulation of the problem Silva , a quantum system is prepared in the ground state of a translation invariant Hamiltonian depending on a control parameter. Such a parameter is suddenly modified (quenched) so that the unitary time evolution of the quantum state is governed by the post-quench Hamiltonian.

In low dimensions, tailored field theoretical CC2006 ; CC2007 ; FM ; SE ; BES ; Delfino2014 ; DV2017 ; KZ ; TM , free fermionic Sac ; Muk07 ; CEF2011 ; CEF2011a ; CEF2011b ; Igloi2012 ; Bucciantini2014b and Bethe Ansatz techniques QA ; DeNa ; DeNa2 ; FCEC14 ; Poz ; Proz ; MPC ; BPC have been exploited for analytical investigation of the unitary dynamics after the quench. For reviews we refer to the special issue CEM and Eis_rev ; Pol_rev .

There is a firm consensus CEM that macroscopic subsystems of a closed infinite system can equilibrate at late times to a stationary quantum density matrix. For one-dimensional integrable models, such a density matrix is a reduced density matrix obtained from a Generalized Gibbs Ensemble (GGE) GGE_rigol ; GGE_rigol2 . This recent theoretical study—especially after KW —has been motivated and boosted by the experimental implementation with ultra-cold atoms of one-dimensional systems which are in almost perfect isolation from the external environment coldatoms .

In the framework pedagogically reviewed in FE2016 , subsystem relaxation towards equilibrium can only occur if the volume of the system is infinite. In a finite system, the energy spectrum might be discrete and the wave function recur Bocchieri1957 , although with periods growing exponentially with the system size. More importantly, in an integrable model, finite size effects can prevent relaxation to a stationary state much earlier Igloi2000 than a recurrence takes place. Suppose that an integrable spin chain is quantized on a ring of length $L$ and $v$ is a maximal LR effective propagation velocity for the quasi-particles. Then, qualitatively, quasi-particles emitted in pairs CC2006 can interfere halfway along the ring at instants of time multiples of $T_{\text{tr}}=\frac{L}{2v}$, which is $O(L)$. Following Sec. 5.2 in FE2016 , we will call the time $T_{\text{tr}}$ around which quasi-particle interference in finite volume prevent relaxation to a time-independent stationary state, the traversal time.

An observable that diagnoses finite-size effects and can provide a quantitative estimate of $T_{\text{tr}}$ is the return amplitude or fidelity Quan2006 ; Yuan2007 ; Rossini2007 ; S2008 ; Zhong2011 ; Zanardi2011 ; MC12 ; Santos14 ; DeLuca14 ; PS14 ; Santos15 ; Mazza2016 ; Poz_echo ; Poz_echo2

[TABLE]

where $|\Psi\rangle$ is the many-body initial state and $\mathcal{H}$, the post-quench Hamiltonian. In particular, destructive or constructive quasi-particle interference is detectable in the logarithm of the return amplitude $\Gamma_{\Psi}(t)=-\frac{1}{L}\log F_{\Psi}(t)$; hereafter, we will refer to $\Gamma_{\Psi}$ as the decay rate. For the quenches considered in this paper, the decay rate relaxes in infinite volume at later times—or rather in a window $1\ll t\ll L$—to a constant Z2010 , that is the limit

[TABLE]

exists and can be calculated. For finite large $L$ instead, when approaching the traversal times, the amplitude of the decay rate oscillations around $\bar{\Gamma}_{\Psi}$ increases; see Fig. 1 for an example that will be discussed in detail. At each traversal time, moreover, the decay rate exhibits local maxima or minima, which signal that the system, being finite, cannot relax to a time-independent stationary state. Actually, after many traversals, the time average of the finite volume return probability $R_{\Psi}\equiv F_{\Psi}^{2}$ converges to its value in the diagonal ensemble Z2010 ; F13 ; LV ; Mazza2016 , $R_{\Psi}^{\text{DE}}$. A correspondent decay rate in infinite volume can be extracted which generally differs Z2010 ; F13 from $\bar{\Gamma}_{\Psi}$ defined earlier; cf. Eq. (48) in Sec. IV of this paper.

Finite-size effects in the return amplitude have been investigated in the last years Happola ; Montes2012 ; Sharma2012 ; Rajak2014 ; KR2017 ; Jafari2017 ; Damski ; Jafari2019 seeking for a revival—i.e. an instant at which $F_{\Psi}=O(1)$—of the initial wave function. However, as pointed out in FE2016 and stressed once more in the Sec. IV of this paper, at the traversal time there is usually not a bona fide revival of the quantum state. The return amplitude remains always exponentially small with respect to the system size, i.e. $O(e^{-L})$ even for reasonably large but finite systems. Indeed a necessary condition for a revival in a time which is $O(L)$ of a state that has an exponential number of non-zero overlaps with the eigenenergy basis, is the fact that the many-body spectrum, in a certain energy unit $\omega=O(1/L)$, is integer spaced. This can happen, for instance, in the continuum limit if non-interacting quasi-particles quantized on a cylinder have a linear dispersion relation. In a lattice model however, the quasi-particles dispersion is never linear on the whole Brillouin zone. The absence of revivials in the critical XY chain was investigated in detail in KR2017 . As a matter of fact, up to now, one can observe physically relevant B2017 examples of revivals only within the context of conformal invariant field theories DS2011 ; Cardy_echo .

The decay rate of the return amplitude in the thermodynamic limit can also show non-analyticities in time that have been dubbed dynamical phase transitions Heyl2013 and are the focus of an intense research activity Pollmann2010 ; Karrasch2013 ; AS ; SB ; Gullo2015 ; Lupo2016 ; Pollmann2016 ; Heyl2017 ; H1 ; H2 ; H3 ; S16 ; BD18 ; TH18 ; Sirker18 ; H4 ; for a recent comprehensive survey we refer to the review Heyl_rev . The nature of such singularities and a possible link with the statistical mechanics are still unsettled AS ; Poz_echo , although they have been observed experimentally DPT_exp2 ; DPT_exp ; see also the viewpoint in G17 .

In this paper, we will study both the finite size effects and the non-analyticities in the thermodynamic limit for the return amplitudes in the XY chain LMS . We will focus on initial states $|\Psi\rangle$ that are eigenstates of the local $z$-magnetization operator; examples will include fully polarized states and the Néel state. These initial states are Gaussian—i.e. the Wick theorem holds KRV2018a ; BTC2018 —and allow local relaxation in the thermodynamic limit to a stationary state SC ; BTC2018 ; moreover, they are experimentally relevant DPT_exp ; DPT_Monroe .

Our analytical results are based on a new formula for the matrix elements of the evolution operator associated with a quadratic fermionic Hamiltonian, which we will derive in Sec. II; see Eq. (14). In Sec. III, we specialized this result to calculate return amplitudes in the XY chain. In Sec. IV, we provide detailed physical applications of the formalism and analyze the finite-size effects and the singularities in the thermodynamic limit in the decay rate of a fully polarized state and of the Néel state. We will then calculate the traversal times and show that they are initial state dependent; we will also determine the instants of time at which singularities in the thermodynamic limit of the decay rate show up. Our conclusions are gathered in Sec. V, and an appendix completes the paper.

II Return amplitude in free fermionic chains

The aim of this section is the derivation of a formula to calculate the Return Amplitude (RA) for an arbitrary many-body state time-evolving with a free fermionic Hamiltonian. We start by considering the following quadratic Hamiltonian in the real space ($T$ denotes transposition below)

[TABLE]

where A and B are $L\times L$ real symmetric and antisymmetric matrices, respectively. The $L$-dimensional column vector c contains the annihilation operators $c_{j}$ with $j=1,\dots,L$. Their Hermitian conjugates are $c_{i}^{\dagger}$ and canonical anticommutation relations hold. Transposition of column vectors is understood if required to define a fermionic quadratic form, such as in Eq. (1) or Eqs. (10)-(12). Following Balian and Brezin Balian1969 , we introduce the $2L$-dimensional column vector $\boldsymbol{\gamma}\equiv(\boldsymbol{c},\boldsymbol{c}^{\dagger})$. The evolution operator acts linearly on the vector $\boldsymbol{\gamma}$ as

[TABLE]

where the $2L\times 2L$ matrix T is obtained by applying the Baker-Campbell-Hausdorff formula (cf. appendix A)

[TABLE]

The $L\times L$ blocks $\textbf{T}_{\alpha\beta}$ ($\alpha,\beta=1,2$) can be written down explicitly if the free fermionic system is translation invariant; in the next section we will provide a physical example. As anticipated in Sec. I, given a quantum state $|\Psi\rangle$, its RA is defined by

[TABLE]

If $|\Psi\rangle$ is an eigenstate of the local fermion occupation numbers, a determinant representation for $F_{\Psi}$ can be obtained through the following route. First, we factorize the evolution operator $e^{-i\mathcal{H}_{\text{free}}t}$ à la Balian-Brezin Balian1969 as

[TABLE]

where X, Y, Z can be calculated from the four blocks of the matrix T in Eq. (3)

[TABLE]

The Balian-Brezin factorization is briefly outlined in the  appendix A; here, we only remark that X and Z are complex antisymmetric matrices. The decomposition (5) is suitable for a coherent state representation of the RA, for analogous manipulations of other quantum matrix elements see also Mizusaki13 . Coherent fermionic states $|\boldsymbol{\xi}\rangle$ are introduced in the usual way through a set of $L$ anticommuting Grassmann numbers $\xi_{i}$, $i=1,\dots L$; namely

[TABLE]

where $|\Omega\rangle$ is the vacuum defined by $\mathbf{c}|\Omega\rangle=0$. Analogously, we can define the left coherent state as $\langle\boldsymbol{\bar{\xi}}|=\langle\Omega|e^{\boldsymbol{\bar{\xi}}\cdot\mathbf{c}}$, where the set of Grassmann numbers $\boldsymbol{\bar{\xi}}$ is independent from the set $\boldsymbol{\xi}$. Recalling that Grassmann numbers anticommute with creation/annihilation operators, it is easy to verify that $\mathbf{c}|\boldsymbol{\xi}\rangle=\boldsymbol{\xi}|\boldsymbol{\xi}\rangle$ and analogously $\langle\boldsymbol{\bar{\xi}}|\mathbf{c^{\dagger}}=\langle\boldsymbol{\bar{\xi}}|\boldsymbol{\bar{\xi}}$. We finally also mention the resolution of the identity and the scalar product between two coherent states, i.e.

[TABLE]

where I is the $2^{L}\times 2^{L}$ identity matrix, acting on the fermionic Fock space and $d\boldsymbol{\bar{\xi}}d\boldsymbol{\xi}\equiv d\bar{\xi}_{n}\dots d\bar{\xi}_{1}d\xi_{1}\dots d\xi_{L}$. Let us now consider the amplitude

[TABLE]

where the state $|\Psi\rangle$ is an eigenstate of the local fermion occupation numbers, that is the overlaps $\langle\boldsymbol{\bar{\xi}}|\Psi\rangle=\prod_{j\in S_{\Psi}}\bar{\xi}_{j}$, where $S_{\Psi}\equiv\{i_{1},\dots,i_{m}\}$ contains the lattice sites with fermion occupation one, for $1\leq i_{1}<i_{2}\dots<i_{m}\leq L$. Inserting two resolutions of the identity in the Eq. (9), and using Eq. (5), we obtain

[TABLE]

The non-trivial matrix element $\langle\boldsymbol{\bar{\eta}}|e^{\mathbf{c^{\dagger}}\bf{Y}\mathbf{c}}|\boldsymbol{\xi}\rangle$ in Eq. (10), can be calculated by observing that, see also Mizusaki13 , $e^{\mathbf{c^{\dagger}}\bf{Y}\mathbf{c}}|\boldsymbol{\xi}\rangle=|e^{\bf{Y}}\boldsymbol{\xi}\rangle$ and then using Eq. (8). Eventually, we end up with the Grassmann integral

[TABLE]

Integrating Eq. (11) first over the Grassmann variables $\bar{\xi}$ and $\eta$ and then over $\xi$ and $\bar{\eta}$ results in

[TABLE]

where the antisymmetric matrix $\tilde{\bf{M}}$ is obtained from

[TABLE]

keeping only the $2m$ lines and columns $\{i_{1},\dots,i_{m}\}$ and $\{i_{1}+L,\dots,i_{m}+L\}$. In practice, these are the lines and columns in correspondence with the lattice sites occupied by a fermion. We can finally conclude (cf. Eq. (4))

[TABLE]

which is an explicit expression for the RA and the main result of this section. The determinant representation in Eq. (12) can be also checked by observing that the initial state can be represented as

[TABLE]

where the coherent state $|\boldsymbol{\xi}_{\Psi}\rangle$ is obtained from the vector $\boldsymbol{\xi}_{\Psi}$ of components $\xi_{j}$ if $j\in S_{\Psi}$ and zero otherwise. The substitution of the expression of $|\Psi\rangle$ discussed above into Eq. (9), leads to Eq. (12) .

Finally, we emphasize that the technique outlined in this section allows calculating the transition amplitude between two arbitrary fermionic states $|\Psi_{i}\rangle$ and $|\Psi_{f}\rangle$, once they are decomposed into eigenstates of the local occupation number operators $c^{\dagger}_{j}c_{j}$. Indeed, let $|\Psi\rangle$ and $|\Psi^{\prime}\rangle$ to be two of such eigenstates, then, $S_{\Psi}=\{i_{1},\dots,i_{m}\}$ and $S_{\Psi^{\prime}}=\{i_{1}^{\prime},\dots,i_{k}^{\prime}\}$ are the set of the lattice sites occupied by one fermion in $|\Psi\rangle$ and $|\Psi^{\prime}\rangle$, respectively. Consequently, to calculate $\langle\Psi^{\prime}|e^{-i\mathcal{H}_{\text{free}}t}|\Psi\rangle$, we can use Eq. (14) where the matrix $\tilde{\boldsymbol{M}}$ is obtained by keeping only the $m+k$ lines and columns $\{i_{1}^{\prime},i_{2}^{\prime},\dots,i_{k}^{\prime}\}$ and $\{i_{1}+L,i_{2}+L,\dots,i_{m}+L\}$. An example will be discussed in Sec. IV (see Eq. (54)) for $m=L$ and $k^{\prime}=0$.

III The XY chain: Balian-Brezin factorization

A systematic application of the Eq. (14) is given by the calculation of the return amplitude in the XY spin chain LMS . The Hamiltonian of the model is

[TABLE]

where $\sigma_{j}^{\alpha}$ $(\alpha=x,y,z)$ are the Pauli matrices and $\gamma$ and $h$ are real parameters conventionally called anisotropy and magnetic field. In the following, we will assume $J>0$ (ferromagnetic model) and $h>0$, the latter restriction is not essential but it simplifies the discussion DR14 when considering odd values of $L$. Although Eq. (14) is valid also for open boundary conditions and can be used for numerical calculations in that case, from now on, we only focus on the periodic case: $\sigma_{j}^{\alpha}=\sigma_{j+L}^{\alpha}$.

The XY chain reduces to the Ising spin chain and free fermions hopping on a lattice for $\gamma=1$ and $\gamma=0$, respectively. Using the Jordan-Wigner transformation, the XY spin chain can be mapped into the quadratic fermionic Hamiltonian in Eq. (1) LMS . To this end, one introduces the fermionic creation operators as $c_{j}^{\dagger}=\prod_{l=1}^{j-1}\sigma_{l}^{z}\sigma_{j}^{+}$, and after nowadays standard manipulations Eq. (16) can be recast in the form

[TABLE]

where we defined $c_{L+1}^{\dagger}\equiv-\mathcal{N}c_{1}^{\dagger}$ and $\mathcal{N}=\pm 1$ is the eigenvalue of the conserved parity operator

[TABLE]

Consequently, the Hamiltonian in Eq. (17) can be rewritten as Eq. (1), if we consider

[TABLE]

It is worth noticing that up to a phase, the quantum state $|\Psi\rangle$, introduced in Sec. II can be identified with a state where spins at positions $i_{1},\dots,i_{m}$ have positive $z$-component. In the following, although this nomenclature is not widespread CEF2011a , we will call the even-parity ($\mathcal{N}=1$) subspace of the Hilbert space the Neveu-Schwarz (NS) sector and the odd-parity ($\mathcal{N}=-1$) subspace, the Ramond sector (R). In order to apply the formalism introduced in Sec. II we must spell out the Balian-Brezin factorization of the evolution operator in the XY chain. This is not difficult because the matrices $\boldsymbol{A}$ and $\boldsymbol{B}$ commute. Indeed, they can be diagonalized by the unitary matrix $[\textbf{U}]_{mk}=\frac{1}{\sqrt{L}}e^{im\phi_{k}}$ with the quantization condition for the momenta

[TABLE]

The above formula can be easily derived after specializing the eigenvalue equation for $\mathbf{A}$ and $\mathbf{B}$ to the first component $i=1$ in Eqs. (19)-(20). Assuming from now on $J=1$, a simple calculation shows that the eigenvalues of these matrices denoted by $\lambda^{\textbf{A}}_{k}$ and $\lambda^{\textbf{B}}_{k}$ are

[TABLE]

with $k=1,\dots,L$ and $\phi_{k}=\frac{2\pi}{L}\left(k-\frac{\mathcal{N}+1}{4}\right)$. Then $2L\times 2L$ matrix equation, (3), becomes a $2\times 2$ matrix equation for the eigenvalues, denoted by $\lambda_{k}^{\textbf{T}_{\alpha\beta}}$, of the four mutually commuting blocks $\boldsymbol{T}_{\alpha\beta}$

[TABLE]

The matrix exponential in Eq. (23) can be explicitly calculated. From the Eqs. (6), we obtain the eigenvalues of the matrices X, Z and $\textbf{Q}\equiv e^{\textbf{Y}}$,

[TABLE]

where $k=1,\dots,L$ and $\varepsilon_{k}$ is a function that coincides with the Bogoliubov quasi-particle dispersion relation Sach_book ; Fran_book ,

[TABLE]

Notice that the eigenvalues of the matrices X and Z and Q in the Eq. (24) are not affected by the choice of the branch of the square root in Eq. (25), which we assume is positive. To apply Eq. (14), which involves the matrix $\tilde{\textbf{M}}$, obtained by removing rows and columns from M, we also need an expression for both X and Q. This is possible since they are both diagonalized by the discrete Fourier transform U and therefore can be written as

[TABLE]

It is useful to observe, cf. Eq. (13), that the matrix Q is symmetric, see also KRV2018a .

It should be emphasized that the formalism discussed here and in Sec. II does not require the determination of the fermionic representation of the spectrum of $\mathcal{H}_{XY}$ at finite $L$, which is actually not a straightforward task XYfacchi ; DR14 . Finally, we observe that unlike the case of the ground state DR14 , the parity sector, to which the initial states considered in this paper belong, will be easily inferred from Eq. (18).

IV Examples: finite-size effects and singularities in the thermodynamic limit

We present some analytic calculations of RA in the XY chain, based on Eq. (14). We will focus on the initial states that are eigenstates of the local $\sigma_{i}^{z}$ operators. We consider initial states denoted by $|\boldsymbol{\mathcal{B}}\rangle$ that are built by repeating an elementary block $\mathcal{B}$ of contiguous spins along the chain. For instance, the state $|\boldsymbol{\uparrow}\rangle\equiv|\uparrow\uparrow\dots\uparrow\uparrow\rangle$ is a fully polarized state in the positive $z$-direction while the state $|\boldsymbol{\downarrow}\boldsymbol{\uparrow}\rangle\equiv|\downarrow\uparrow\dots\downarrow\uparrow\rangle$ is the Neél state.

It is also useful to introduce the decay rate

[TABLE]

which, for our protocols, remains finite in the thermodynamic limit; see Viti_in ; JM_return for important exceptions in the context of inhomogeneous initial states. We note that a decay rate can also be defined for any quantum transition probability between two states $|\Psi_{i}\rangle$ and $|\Psi_{f}\rangle$ on a chain of length $L$ as: $-\frac{1}{L}\log|\langle\Psi_{f}|\Psi_{i}(t)\rangle|$.

IV.1 Return Amplitude for the Fully Polarized Initial State

The fully polarized initial state belongs to the NS sector of the Hilbert space of the XY chain, independently of the parity of $L$. It is therefore understood that $\phi_{k}=\frac{2\pi}{L}(k-1/2)$ and $k=1,\dots,L$ whenever needed. Applying Eq. (14) and observing that $|\text{Pf}(\textbf{M})|=1$, it follows

[TABLE]

which implies

[TABLE]

In the limit $L\rightarrow\infty$, we can formally replace $\sum_{k=1}^{L}\rightarrow\frac{L}{2\pi}\int_{0}^{2\pi}d\phi$ and obtain

[TABLE]

IV.1.1 Stationary Value of the Decay Rate

As anticipated in Sec. I, the thermodynamic limit of the decay rate in the Eq. (30) approaches a constant for large times

[TABLE]

The value of the constant in the Eq. (31) can be obtained explicitly: for $h>1$ the logarithm in the Eq.(30) can be expanded in a convergent power series for any value of $t$. The infinite time limit is then calculated by retaining only the time-independent terms in the power series of the logarithm or more formally applying the Riemann-Lebesgue lemma to drop all the oscillating contributions. This is of course equivalent to replace the large time limit with the time average in Eq. (31); it turns out after resuming all the time-independent terms

[TABLE]

A similar result has been obtained in S2008 ; Z2010 for a different decay rate in the Ising spin chain. Although derived for $h\geq 1$, the Eq.(32) holds for $h<1$ as well. In the time window $1\ll t\ll L$, the decay rate $\Gamma_{\boldsymbol{\uparrow}}$ in Eq. (29) gets closer to Eq. (32); an example is given for $L=1000$, $\gamma=1/2$ and $h=1$ in Fig. 1: compare the blue oscillating curve with the dashed horizontal black line for times sufficiently smaller than the system size.

IV.1.2 Traversals

For large but finite $L$, the decay rate relaxes to its stationary value in the thermodynamic limit until a time $T_{\text{tr}}=O(L)$. Within a qualitative quasi-particle picture CC2006 , quasi-particles of opposite momentum are emitted in pairs at $t=0$ and travel around the ring of length $L$ with an effective maximal velocity $v$. They can then interfere constructively or destructively at the time $T_{\text{tr}}=\frac{L}{2v}$ when traversing half of the ring. According to the nomenclature introduced in Sec. 5.2 of Ref. FE2016 , we will refer to $T_{\text{tr}}$ as a traversal time.

A formal explanation of this effect and a quantitative estimation for $T_{\text{tr}}$ for reasonably large $L$ is provided by the following argument. If $L$ is large enough, we can approximate the quasi-particle dispersion $\varepsilon_{k}$, given in Eq. (25) by a continuous function $\varepsilon(\phi)$. Now suppose that $\varepsilon(\phi)$ has an inflection point at $\phi=\phi_{\text{in}}\in[0,2\pi]$, i.e.

[TABLE]

where we defined $v_{\text{in}}=\varepsilon^{\prime}(\phi_{\text{in}})$.

For large $L$ , the interval $\Delta_{\phi}\equiv(\phi-\phi_{\text{in}})$ might still contain a macroscopic fraction of modes. For all the discrete modes $\phi_{k}$ in the interval $\Delta_{\phi}$, we have $\phi_{k}-\phi_{\text{in}}=\frac{2\pi n_{k}}{L}+\delta_{0}$, where $0\leq\delta_{0}<\frac{2\pi}{L}$ and $n_{k}$ is an integer. Extracting the discrete modes in the interval $\Delta_{\phi}$ from the sum in the Eq. (29), we can rewrite their contribution up to the first order in $\Delta_{\phi}$ as

[TABLE]

At times which are positive integer multiples of $T_{\text{tr}}=\frac{L}{2v_{\text{in}}}$, the elements of the sum in the Eq. (34) are all equal, allowing larger fluctuation of the decay rate around the constant in Eq. (32). Rightly at the traversal time, the decay rate displays local maxima or minima which can be more pronounced depending on the validity of the linear approximation in the Eq. (33) and the values of $\varepsilon(\phi_{\text{in}})$ and $\delta_{0}$ in the Eq. (34).

For instance, if $\varepsilon(\phi_{\text{in}})=0$ and $\delta_{0}=\frac{\pi}{L}$, as long as $\lim_{\phi\rightarrow\phi_{\text{in}}}|\lambda^{\textbf{B}}(\phi)|/\varepsilon(\phi)=1$, at time $t=T_{\text{tr}}$, each term in the Eq. (34) diverges logarithmically. In the actual spin chain, this logarithmic divergence is converted into a maximum of the decay rate, which is clearly visible along the critical line $h=1$ (here $\phi_{\text{in}}=0$, $\varepsilon(0)=0$, $v_{\text{in}}=\varepsilon^{\prime}(0)=\gamma$); examples are shown in the Fig. 1 for $\gamma=1/2$ and for $\gamma=1$ in Fig. 2. Such a maximum of the decay rate indicates (see below the Eq. (54) for an explanation) that the time-evolved state is closer to a state fully polarized in the negative $z$-direction.

Eq. (29) for even $L$ and discrete momenta $\phi_{k}$ in the NS sector can also be applied to a fully negative polarized initial state, denoted hereafter by $|\downarrow\dots\downarrow\downarrow\rangle\equiv|\boldsymbol{\downarrow}\rangle$. However, if $L$ is odd, this is no longer true. In such a case, the fully negative polarized initial state belongs to the R sector of the $XY$ chain (recall Eq. (18)) and the Eq. (29) gives its decay rate provided the discrete momenta $\phi_{k}$ are chosen accordingly. In particular, for odd $L$, along the critical line $h=1$, the decay rate $\Gamma_{\boldsymbol{\downarrow}}$ is minimum, contrary to $\Gamma_{\boldsymbol{\uparrow}}$, when $t$ is close to $T_{\text{tr}}=\frac{L}{2\gamma}$. Indeed, for the R sector $\delta_{0}=0$ and all the terms in the Eq. (34) are now vanishing. We refer to Fig. 2 where this prediction has been checked against exact diagonalization for the critical Ising spin chain with $L=15$ spins; see the blue curve.

Notice that, for odd $L$ and $h=1$, at the first traversal, $F_{\boldsymbol{\downarrow}}$ has a maximum; see the inset in Fig. 2. Such a maximum, however, should not be confused with a revival of the quantum state $|\boldsymbol{\downarrow}\rangle$, which would happen FE2016 , if the RA were of order one, independently of the system size. Instead, the Eq. (34) suggests that for large enough $L$ at the traversal time, $F_{\boldsymbol{\downarrow}}(T_{\text{tr}})=O(e^{-L})$, unless only a finite set of momenta do not fall in the interval $\Delta_{\phi}$; again see the inset in Fig. 2. The latter requirement is met, in physical models, only if the quasi-particle energy levels are all equally spaced, which is not the case for the $XY$ spin chain.

Finally, we remark that quantum revivals in the $XY$ chain of a fully polarized initial state are only possible in the limiting case $h=0$, $\gamma=1$. Indeed, for $h=0$ and $\gamma=1$, the Hamiltonian in Eq. (16) trivializes and it is an exercise in the quantum mechanics to obtain

[TABLE]

which agrees with the Eq. (27) for the same parameter values and shows that the fully polarized state is reproduced at the revival times $T_{\text{rev}}=2\pi n$, $n=1,2,\dots$.

IV.1.3 Dynamical Phase Transitions

At finite $L$, the RA of a quantum state vanishes when the time-evolved initial state is orthogonal to the final one (and vice-versa). For instance from the Eq. (28), it follows that a necessary and sufficient condition for $F_{\boldsymbol{\uparrow}}=0$ is that it exists a $t_{k}$ such that $|\lambda_{k}^{\textbf{Q}}(t_{k})|^{-1}=0$. However, for any value of $\gamma$,

[TABLE]

can vanish only for $h<1$ and $\phi_{k_{c}}=\arccos(h)$ when

[TABLE]

The finite-size logarithmic divergences at $t^{(n)}_{k_{c}}$ are downgraded, for $L\rightarrow\infty$, to non-differentiable points of the decay rate in the Eq. (30), which have been termed dynamical phase transitions Heyl2013 . Indeed the same values for $t^{(n)}_{k_{c}}$ in the Eq. (37) were already obtained for the fully polarized initial state in the Ising spin chain ($\gamma=1$) in Heyl2013 , relating them to so-called Fisher zeros Fisher of an analytically continued RA. In summary, Eq. (37) shows that the generalization of the results in Heyl2013 for a fully polarized initial state to $\gamma\not=1$ is straightforward and leads to a qualitatively similar behaviour of the decay rate.

IV.2 Return Amplitude for the Néel Initial State

From the formalism of Sec. II, it is also possible to determine analytically the RA for the Néel state. We mention that the RA for a quench from the Néel state in the XY chain has been calculated in AS ; Mazza2016 at $\gamma=0$. It has been also worked out in the thermodynamic limit for the gapless XXZ spin chain in Poz_echo ; Poz_echo2 by a Quantum Transfer Matrix approach K_rev ; P14 ; an earlier related numerical study was carried out, again for the XXZ spin chain, in AS .

According to the notations introduced at the beginning of this section, the Néel state is $|\boldsymbol{\downarrow\uparrow}\rangle$ and, assuming $L=2\ell$, it will belong to the NS sector if $\ell$ is even (i.e. $L=4n$) and to the R sector if $\ell$ is odd (i.e. $L=4n+2$). In order to apply Eq. (14), we observe that the matrices $\tilde{\boldsymbol{X}}$, $\tilde{\boldsymbol{Z}}=-\tilde{\boldsymbol{X}}$ and $\tilde{\boldsymbol{Q}}$ obtained from $\boldsymbol{X}$ and $\boldsymbol{Q}$ removing the odd lines and columns are again simultaneously diagonalized by the unitary matrix $[\tilde{\boldsymbol{U}}]_{mk}=\frac{1}{\sqrt{\ell}}e^{2i\phi_{k}m}$, with $k,m=1,\dots,\ell$. This can be proven as follows:

From Eq. (26), we obtain the matrix $[\tilde{\textbf{X}}]_{nm}=\frac{1}{2\ell}\sum_{k=1}^{2\ell}\lambda_{k}^{\textbf{X}}e^{2i\phi_{k}(n-m)}$. Let us apply such a matrix to the vector $\tilde{u}^{(k^{\prime})}_{m}\equiv[\tilde{\boldsymbol{U}}]_{mk^{\prime}}$, then, we obtain

[TABLE]

where we used $u^{(k^{\prime})}_{n}=u^{(k^{\prime}+\ell)}_{n}$. An analogous calculation can be performed with the matrix $\tilde{\textbf{Q}}$ in Eq. (26), thus proving that $\tilde{\textbf{X}}$ and $\tilde{\textbf{Q}}$ commute.

Coming back to the determination of the decay rate for the Néel state in the XY chain, we recall that ${\rm Pf}(\tilde{\textbf{M}})=\pm\sqrt{\det{\tilde{\textbf{M}}}}$ and

[TABLE]

where the matrix $\tilde{\textbf{Z}}+\tilde{\textbf{Q}}\tilde{\textbf{X}}^{-1}\tilde{\textbf{Q}}$ in the Eq. (39) is the Schur complement Sc of the block matrix $\tilde{\textbf{M}}$. Since all the matrices in the Eq. (39) commute (cf. Eq. (38) and the discussion above), we obtain, from the Eq. (14), the decay rate

[TABLE]

where $\lambda_{k}^{\tilde{\textbf{X}}/\tilde{\textbf{Q}}}=\frac{1}{2}(\lambda_{k}^{\textbf{X}/\textbf{Q}}+\lambda^{\textbf{X}/\textbf{Q}}_{k+\ell})$ and $\Gamma_{\boldsymbol{\uparrow}}$ is given in the Eq. (29). After some lengthy but straightforward trigonometric manipulations, the Eq. (40) can be eventually rewritten as

[TABLE]

where we introduced the notations

[TABLE]

When $h\leq 1$, for $\phi_{k}<\arccos(h)$ or $\phi_{k}>2\pi-\arccos(h)$ the angle $\phi^{b}_{k}$ in the Eq. (42) is defined as $\phi^{b}_{k}=\arctan\left(\frac{\gamma\sin(\phi_{k})}{-h+\cos(\phi_{k})}\right)\mp\pi$ respectively. Eq. (41) with $L=4n$ and $h>1$ simplifies for $\gamma=0$ (i.e. free fermions hopping on the lattice), leading to AS

[TABLE]

which is independent of $h$; see also Mazza2016 for an analogous expression with open boundary conditions. We have checked the Eq. (41) with exact diagonalization in the Ising spin chain up to $L=14$ spins; see the inset in Fig. 3 where the red dots are obtained by exact diagonalization and the continuous black line is Eq. (41) with momenta in the R sector. In the thermodynamic limit Eq. (41) becomes

[TABLE]

where $\varepsilon^{\pm}(\phi)=\frac{1}{2}(\varepsilon(\phi)\pm\varepsilon(\phi+\pi))$ and $\bar{\phi}^{b}(\phi)=\frac{1}{2}(\phi^{b}(\phi)+\phi^{b}(\phi+\pi))$.

IV.2.1 Stationary Value of the Decay Rate

At late times, the infinite volume Néel decay rate approaches a constant Z2010

[TABLE]

The value of such a constant can be calculated as briefly outlined above the Eq. (32) for the fully polarized initial state. It turns out

[TABLE]

where ${}_{2}F_{1}$ is the Gauss hypergeometric function. Although convergence is rather slow, the sum and the integral in the Eq. (46) can be evaluated numerically. At $\gamma=0$, Eq. (46) simplifies drastically and $\bar{\Gamma}_{\downarrow\uparrow}=\log(2)/2$. It is instructive to compare this result with the one that could be obtained taking the time-average of the finite volume return probability in Eq. (43). Such a time average corresponds to the value of the return probability in the diagonal ensemble Mazza2016

[TABLE]

Then we can extract the infinite volume decay rate

[TABLE]

which is half of Eq. (46) at $\gamma=0$. As anticipated in the Introduction, this example shows clearly that the thermodynamic limit and the infinite time limit do not commute when calculating the return probability and its decay rate.

IV.2.2 Traversals

Qualitatively FE2016 , at the first traversal time $T_{\text{tr}}=O(L)$ the decay rate in Eq. (41) stops relaxing towards the constant in Eq. (46). We can estimate $T_{\text{tr}}$ as follows: According to the discussion around Eq. (33), we look for inflection points of the functions $\varepsilon^{\pm}$ in Eq. (44). It turns out that $\varepsilon^{-}$ has, irrespectively of $\gamma$ and $h$, an inflection point at $\phi_{\text{in}}=\pi/2$ where $\bar{\phi}^{b}=0$. Therefore, quite interestingly, Eq. (34) can be applied as if on the critical line $h=1$ for the fully polarized initial state, upon replacing $\varepsilon\rightarrow\varepsilon^{-}$ and

[TABLE]

Then for $L=4n$ and momenta $\phi_{k}$ in the NS sector $\delta_{0}=\frac{\pi}{L}$, the Eqs. (34) and (49) predict a maximum of the decay rate for $\Gamma_{\boldsymbol{\downarrow\uparrow}}$ at the traversal time $T_{\text{tr}}=\frac{L}{2v_{\text{in}}}$. The same conclusion applies for $L=4n+2$ and momenta $\phi_{k}$ in the R sector. These predictions have been checked in the Fig. 3 for $\gamma=1$ and $h=2$. In particular, we note that the location of the first traversal time for the Néel state is retarded with respect to the one of the fully polarized state. We can then conclude that $T_{\text{tr}}$ is state-dependent, and bounded by $\frac{L}{2v_{\text{max}}}$ with $v_{\text{max}}$ the maximum allowed quasi-particle dispersion given in Eq. (25).

IV.2.3 Dynamical Phase Transitions

We can determine the instants of time at which the finite-size decay rate $\Gamma_{\boldsymbol{\downarrow\uparrow}}(t)$ is logarithmically divergent. We observe that the argument of the logarithm in the Eq. (41) can be rewritten as

[TABLE]

which can vanish if and only if each of the two squares in the sum vanishes. In particular,

• •

If $h>1$, $-\pi/2<\bar{\phi}^{b}_{k}<\pi/2$, then we can have $\bar{\phi}_{k_{c}}^{b}=0$ for $\phi_{k_{c}}=\{0,~{}\pi/2,~{}\pi\}$ in the interval $[0,\pi]$. Since $\varepsilon^{-}(\phi)$ is zero for $\phi=\pi/2$ and $\varepsilon^{-}(\pi)=-\varepsilon^{-}(0)=1$ for $h>1$, we obtain a first series of singularities at

[TABLE]

with $\phi_{k_{c}}=\pi$. However, this is not enough, there might exist values of the angle $\phi_{k^{*}}\in[0,\pi]$ such that

[TABLE]

and thus leading to a dynamical phase transition at $t=t_{nm}>0$. This is illustrated for the Ising spin chain with $h=2$ in the Fig. 4: black dashed lines denote the first two terms in the series (51) $(n=0,1)$ while red dashed lines indicate the first four singularities predicted by Eq. (52).

• •

If $h\leq 1$, we can have $-\pi/2\leq\bar{\phi}^{b}_{k}\leq\pi/2$. However $\bar{\phi}^{b}_{k_{c}}=0$ only for $\phi_{k_{c}}=\pi/2$ while $\bar{\phi}^{b}_{k_{c}}=\mp\pi/2$ for $\phi_{k_{c}}=0,\pi$ respectively. Moreover, since for $h\leq 1$, $\varepsilon^{-}(\pi/2)=0$ and $\varepsilon^{+}(0)=\varepsilon^{+}(\pi)=1$, we obtain a series of dynamical phase transitions at times

[TABLE]

which are the same as the ones given in Eq. (51) (again $\phi_{k_{c}}=\pi$). If $h\leq 1$, then, the argument leading to Eq. (52) also remains intact.

This analysis supports the original claim AS that dynamical phase transitions in the Néel decay rate are not specific of quenches crossing the critical line $h=1$; they rather occur both quenching into the paramagnetic ($h>1$) or ferromagnetic phase ($h\leq 1$). Moreover, the location of the singularities depends non-trivially on the quasi-particle spectrum, see the Eq. (52).

Finally, we observe that for finite $L$, the critical momenta $\phi_{k_{c}}$ and $\phi_{k^{*}}$ in the Eqs. (51) and  (52) might not be allowed by the NS or R quantization. In such a case the finite size Neel decay rate $\Gamma_{\downarrow\uparrow}$ will not be singular at the corresponding critical times. This is again shown in the inset of the Fig. 3 when $\gamma=1$ and $h=2$. For $L=14$, in the Ramond sector one can verify that $\phi_{k_{c}}=\pi$ is an allowed angle and therefore the Néel decay rate will be logarithmically divergent at times given in Eq. (51). In the inset of the Fig. (3), there is also a visible peak of the the decay rate at $t_{13}=5.349\dots$ given in Eq. (52), see also the Fig. 4. It is indeed possible to show that $\phi_{k=5}$ is close (less than half of a percent) to the angle $\phi_{k^{*}}$, i.e. solution of the Eq. (52), for $n=1$ and $m=3$.

IV.3 The Flip Amplitude

The method introduced in Sec. II, allows also determining the decay rate of the transition probability between the two fully polarized initial states: $|\boldsymbol{\uparrow}\rangle$ and $|\boldsymbol{\downarrow}\rangle$. The transition probability is zero for odd $L$ since the two states belong to different sectors of the Hilbert space, whereas for $L$ even and $\phi_{k}$ in the NS sector (observe that $\lambda_{k}^{\mathbf{Z}}=\frac{\varepsilon_{k}\lambda_{k}^{\textbf{X}}}{i\lambda_{k}^{\textbf{B}}\sin(t\varepsilon_{k})}$)

[TABLE]

The analysis of the traversals parallels the one for the fully polarized initial state given in Sec. IV.1.2, therefore we will not repeat it here. We limit to observe that along the critical line $h=1$, at the first traversal time $T_{\text{tr}}=\frac{L}{2\gamma}$ the decay rate $\Gamma_{\text{flip}}(T_{\text{tr}})$ has a minimum. This signals that the time evolved state $e^{-i\mathcal{H}_{XY}t}|\boldsymbol{\uparrow}\rangle$ is closer, although still exponentially far for large $L$, to a state fully polarized in the negative $z$-direction.

IV.3.1 Dynamical Phase Transitions

We finally discuss dynamical phase transitions for the flip amplitude; to our best knowledge, this case has been never studied analytically before. At finite $L$, the decay rate $\Gamma_{\text{flip}}$ in the Eq. (54) is logarithmically divergent at

[TABLE]

and $k=1,2,\dots L$; see the blue curve in the Fig. 5, obtained for $\gamma=2$, $h=1/2$ and $L=100$. A similar behaviour has been also found in AS for a different physical setting: a quench from the Néel state in a free fermionic chain with staggered magnetic field; cf. Fig. 4 there.

In the thermodynamic limit $L\rightarrow\infty$, we can formally replace $\sum_{k=1}^{L}\rightarrow\frac{L}{2\pi}\int_{0}^{2\pi}$ in the Eq. (54); the result is shown by the red curve in the Fig. 5. We then conjecture that non-analyticities in the integral $\lim_{L\rightarrow\infty}\Gamma_{\text{flip}}(t)$ occur at discrete times which are multiples of $t_{\min}\equiv\frac{\pi}{\max_{\phi\in[0,2\pi]}\varepsilon(\phi)}$, which is the first singularity in the series of Eq. (55), and of $t_{\text{edge}}\equiv\frac{\pi}{\varepsilon(\pi)}$. The latter are in correspondence with a single-particle energy at the edge of the spectrum. We have checked the last statement numerically in the XY chain, however we do not have a satisfactory analytical understanding of such a claim. An example is again given in the Fig. 5 for $\gamma=2$ and $h=1/2$: dashed red lines denote multiples of $t=nt_{\min}$ ($n=1,2,3$), while dashed black lines correspond to $t=nt_{\text{edge}}$ ($n=1,2$). Finally, notice that if the initial state is a cat state $|\text{cat}\rangle=2^{-1/2}\left(|\boldsymbol{\uparrow}\rangle+|\boldsymbol{\downarrow}\rangle\right)$, then

[TABLE]

where $L=2\ell$. One can find the Eq. (56) simply by observing that the interference term is bounded for large $L$ by $2e^{-L(\Gamma_{\boldsymbol{\uparrow}}+\Gamma_{\text{flip}})}$, which is always exponentially small in the system length compared to $e^{-2L\min\{\Gamma_{\boldsymbol{\uparrow}},\Gamma_{\text{flip}}\}}$. This non-equilibrium protocol was considered in the experiment DPT_exp and proposed theoretically in S16 for an Ising spin chain with the long-range couplings. However, no fully analytical prediction was presented so far.

V Conclusions

In this paper, we studied the return amplitude in the XY spin chain for the eigenstates of the local $z$-component magnetization operator. We obtained analytic expressions for the fully polarized states and the Néel state, which generalize the previous results in Quan2006 ; AS ; Mazza2016 to arbitrary values of the anisotropy parameter $\gamma$. These formulas have been derived from a determinant representation for the matrix elements of the evolution operator of a quadratic fermionic Hamiltonian, see Eq. (14). We then focussed on the analysis of the finite-size effects in the return amplitude showing that they are signalled by the so-called traversals FE2016 , whose features depend also on the parity of the length of the spin chain. In particular, at the first traversal, the decay rate might show a maximum or a minimum depending on the quantization sector (NS or R) to which the initial state belongs. Analogously, we provided evidence that traversal times are also initial state dependent. Our results have been tested with exact diagonalization methods for the Ising spin chain up to $L=15$. We hopefully made it clear that at a traversal time the return amplitude is expected to be $O(e^{-L})$ and therefore such a finite-size effect does not lead to a revival of the quantum state; in agreement with the discussion in FE2016 . Finite-size effects in the evolution of the entanglement entropy Fagotti2011 ; KRV2018a ; BTC2018 could be in principle studied similarly. Our approach is also suitable to extract overlaps between the initial states examined in this paper and the eigenstates of the XY Hamiltonian. Analytical results for the overlaps obtained in Mazza2016 for the Néel state could be extended to $\gamma\not=0$ as well.

We have moreover analyzed the thermodynamic limit of the logarithm of the return amplitude, namely the decay rate, and identified analytically the instants of time at which it might develop non-analyticities. For the Néel decay rate, singularities in time do not follow a periodic pattern and they appear independently of the final Hamiltonian. Our results complement in this respect the numerical study of AS and the Bethe Ansatz analysis in Poz_echo ; Poz_echo2 .

The analytical results for $f_{\Psi}(t)=\langle\Psi|e^{-i\mathcal{H}_{XY}t}|\Psi\rangle$ contained in this paper, could be also exploited to extract, upon analytic continuation to imaginary time $t=i\tau$, universal boundary entropies along the critical line $h=1$ AL . We hope to come back on this problem in the near future.

Acknowledgements. KN acknowledges the supports by National Science Foundation under Grant No. PHY-1620555 and DOE grant DE-SC0018326. MAR acknowledges the support from CNPq. JV thanks Rodrigo Pereira for a discussion.

Appendix A Balian Brezin factorization

In this Appendix we review the Balian Brezin factorization Balian1969 . Consider a complex quadratic fermionic form (again transposition of column vectors is understood)

[TABLE]

where the $2L$-dimensional column vector $\boldsymbol{\gamma}=(\mathbf{c},\mathbf{c}^{\dagger})$ and the operators $c^{\dagger}_{i}$ and $c_{j}$ obey canonical anticommutation relations. The components of the vector $\boldsymbol{\gamma}$ satisfy

[TABLE]

It is also useful to observe that H which in (57) is a complex matrix can be taken complex antisymmetric. If H is not antisymmetric we can write $\textbf{H}=\textbf{H}_{a}+\textbf{H}_{s}$, being $\textbf{H}_{a/s}$ the antisymmetric/symmetric parts of it. Now substituting into Eq. (57) we get

[TABLE]

and all the formulas that follow have to be modified accordingly. The transformation $\mathcal{F}=e^{\mathcal{H}}$ acts linearly on the fermion $\boldsymbol{\gamma}$. Indeed, we can apply Baker-Campbell-Hausdorff formula and the commutation relations to prove that

[TABLE]

being $\textbf{T}=e^{\sigma\textbf{H}}$. The matrix T satisfies $\textbf{T}\boldsymbol{\sigma}\textbf{T}^{T}=\boldsymbol{\sigma}$.

Let us now consider two transformations $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ of the same type as in the Eq. (57). Their composition applied to $\gamma$ yields

[TABLE]

However Balian1969 , since complex antisymmetric fermionic quadratic forms form a Lie algebra by Eq. (60), there should exist complex antisymmetric matrices $\textbf{H}_{1}$ and $\textbf{H}_{2}$ such that $\textbf{T}_{1,2}=e^{\sigma\textbf{H}_{1,2}}$ and $e^{\sigma\textbf{H}}=e^{\sigma\textbf{H}_{2}}e^{\sigma\textbf{H}_{1}}$, provided $\mathcal{F}=\mathcal{F}_{1}\mathcal{F}_{2}$. Let us then consider the operator $\mathcal{F}$ and its associated matrix

[TABLE]

According to the discussion above, we can find a factorization of $\mathcal{F}$ in the form

[TABLE]

such that $\mathcal{F}_{3}$ (resp. $\mathcal{F}_{1}$) only contains $\mathbf{c}$ (resp. $\mathbf{c}^{\dagger}$) operators. This implies the equation among matrices

[TABLE]

which can be solved with the following result

[TABLE]

The antisymmetry of the matrix $\bf{Z}$ follows from the antisymmetry of $\textbf{T}_{21}\textbf{T}_{22}^{T}$, while the antisymmetry of $\bf{X}$ from the antisymmetry of $\textbf{T}_{12}\textbf{T}_{11}^{T}$. Furthermore we can also check that $\textbf{T}_{12}\textbf{T}_{21}^{T}+\textbf{T}_{11}\textbf{T}_{22}^{T}=\bf{1}$. The Eq. (65) is the so-called Balian-Brezin factorization Balian1969 . In the specific case discussed in Sec. II one can then apply Eq. (65) to the complex antisymmetric matrix

[TABLE]