Complex semiclassical analysis of the Loschmidt amplitude and dynamical quantum phase transitions

TL;DR

This paper introduces a complex semiclassical method to analyze the Loschmidt amplitude in spin systems, revealing dynamical quantum phase transitions through theoretical and numerical approaches.

Contribution

It develops a new semiclassical computational technique for the Loschmidt amplitude applicable to generic spin systems, enabling analysis of dynamical phase transitions.

Findings

Dynamical phase diagrams identified for the infinite-range transverse Ising model.

Dynamical transitions depend on initial conditions and quench parameters.

Semiclassical results are validated by numerical experiments.

Abstract

We propose a new computational method of the Loschmidt amplitude in a generic spin system on the basis of the complex semiclassical analysis on the spin-coherent state path integral. We demonstrate how the dynamical transitions emerge in the time evolution of the Loschmidt amplitude for the infinite-range transverse Ising model with a longitudinal field, exposed by a quantum quench of the transverse field $\Gamma$ from $\infty$ or $0$. For both initial conditions, we obtain the dynamical phase diagrams that show the presence or absence of the dynamical transition in the plane of transverse field after a quantum quench and the longitudinal field. The results of semiclassical analysis are verified by numerical experiments. Experimental observation of our findings on the dynamical transition is also discussed.

Equations

$$\hat{\mathcal{H}}=-\frac{NJ}{2}{\left(\frac{1}{N}\sum _{i=1}^N{\sigma }_i^z\right)}^2-\Gamma \sum _{i=1}^N{\sigma }_i^x-h\sum _{i=1}^N{\sigma }_i^z,$$

$$\mathcal{L}\left(t|\psi \right)=⟨\psi |e^{-it\hat{\mathcal{H}}}|\psi ⟩,$$

$$f\left(t|\psi \right)=-\frac{1}{N}\log \mathcal{L}\left(t|\psi \right).$$

$$|\theta ,\phi ⟩=e^{ib}\left(e^{-i\frac{\phi }{2}}\cos \frac{\theta }{2}|\uparrow ⟩+e^{i\frac{\phi }{2}}\sin \frac{\theta }{2}|\downarrow ⟩\right),$$

$$⟨\theta ,\phi |\mathbf{\sigma }|\theta ,\phi ⟩=\left(\sin \theta \cos \phi ,\sin \theta \sin \phi ,\cos \theta \right).$$

$${\int }_{-1}^{1}d\cos \theta {\int }_0^{2\pi }\frac{d\phi }{2\pi }|\theta ,\phi ⟩⟨\theta ,\phi |=|\uparrow ⟩⟨\uparrow |+|\downarrow ⟩⟨\downarrow |=I,$$

$$G(t|{\Omega }^{\prime },{\Omega }^{\prime \prime })={\int }_{{\Omega }^{\prime }}^{{\Omega }^{\prime \prime }}\prod _{i=1}^N\mathcal{D}\cos {\theta }_i\mathcal{D}{\phi }_i{e}^{S[\mathbf{\theta },\mathbf{\phi }]}.$$

$$S[\mathbf{\theta },\mathbf{\phi }]=i{\int }_0^{t}ds\left\{\frac{1}{2}\sum _i{\dot{\phi }}_i(s)\cos {\theta }_i(s)-\mathcal{H}\left(\mathbf{\theta },\mathbf{\phi },s\right)\right\},$$

$$\frac{1}{2}{\dot{\theta }}_i\sin {\theta }_i=\frac{\partial \mathcal{H}}{\partial {\phi }_i},\frac{1}{2}{\dot{\phi }}_i\sin {\theta }_i=-\frac{\partial \mathcal{H}}{\partial {\theta }_i}.$$

$$G(t|{\Omega }^{\prime },{\Omega }^{\prime \prime })\sim \sum _{\nu }{A}_{\nu }{e}^{S_{\mathrm{cl}}[{\bar{\mathbf{\theta }}}^{(\nu )},{\bar{\mathbf{\phi }}}^{(\nu )}]},$$

$$W[\mathbf{\theta },\mathbf{\phi }]=-\frac{1}{4}m{\int }_0^{t}ds\sum _i\left({\dot{\theta }}_i^2+{\dot{\phi }}_i^2{\sin }^2{\theta }_i\right),$$

$$\frac{1}{2}{\dot{\theta }}_j\sin {\theta }_j$$

$$=\frac{\partial \mathcal{H}}{\partial {\phi }_j}+\frac{i}{2}m\left({\ddot{\phi }}_j{\sin }^2{\theta }_j+2{\dot{\theta }}_j{\dot{\phi }}_j\sin {\theta }_j\cos {\theta }_j\right),$$

$$\frac{1}{2}{\dot{\phi }}_j\sin {\theta }_j$$

$$=-\frac{\partial \mathcal{H}}{\partial {\theta }_j}-\frac{i}{2}m\left({\ddot{\theta }}_j{\sin }^2{\theta }_j-{\dot{\phi }}_j^2\sin {\theta }_j\cos {\theta }_j\right).$$

$$\tan \left(\frac{{\bar{\theta }}_i(0)}{2}\right)e^{i{\bar{\phi }}_i(0)}=\tan \left(\frac{{\theta }_i^{\prime }}{2}\right)e^{i{\phi }_i^{\prime }},$$

$$\tan \left(\frac{{\bar{\theta }}_i(t)}{2}\right)e^{-i{\bar{\phi }}_i(t)}=\tan \left(\frac{{\theta }_i^{\prime \prime }}{2}\right)e^{-i{\phi }_i^{\prime \prime }}.$$

$${\zeta }_j(s)=\tan \left(\frac{{\theta }_j(s)}{2}\right)e^{i{\phi }_j(s)},$$

$${\eta }_j(s)=\tan \left(\frac{{\theta }_j(s)}{2}\right)e^{-i{\phi }_j(s)}.$$

$${\zeta }_j(0)={\zeta }_j^{\prime }\equiv \tan \left(\frac{{\theta }_j^{\prime }}{2}\right)e^{i{\phi }_j^{\prime }},$$

$${\eta }_j(t)={\eta }_j^{\prime \prime }\equiv \tan \left(\frac{{\theta }_j^{\prime \prime }}{2}\right)e^{-i{\phi }_j^{\prime \prime }},$$

$$\displaystyle\langle\theta_{j},\varphi_{j}|\bm{\sigma}_{j}|\theta_{j},\varphi_{j}\rangle=\frac{1}{1+\zeta_{j}\eta_{j}}\left(\begin{array}[]{c}\zeta_{j}+\eta_{j}\\ -i(\zeta_{j}-\eta_{j})\\ 1-\zeta_{j}\eta_{j}\end{array}\right),$$

$${\dot{\zeta }}_j=-i(1+{\zeta }_j{\eta }_j{)}^2\frac{\partial \mathcal{H}}{\partial {\eta }_j},$$

$${\dot{\eta }}_j=-i(1+{\zeta }_j{\eta }_j{)}^2\frac{\partial \mathcal{H}}{\partial {\zeta }_j}.$$

$$\displaystyle e^{S_{\rm cl}[\bar{\bm{\zeta}},\bar{\bm{\eta}}]}=\prod_{j=1}^{N}\Biggl{\{}\sqrt{\frac{(1+\bar{\zeta}_{j}(0)\bar{\eta}_{j}(0))(1+\bar{\zeta}_{j}(t)\bar{\eta}_{j}(t))}{(1+\zeta_{j}^{\prime}\eta_{j}^{\prime})(1+\zeta_{j}^{\prime\prime}\eta_{j}^{\prime\prime})}}$$

$$\displaystyle\times\left(\frac{\zeta_{j}^{\prime}\eta_{j}^{\prime}\zeta_{j}^{\prime\prime}\eta_{j}^{\prime\prime}}{\bar{\zeta}_{j}(0)\bar{\eta}_{j}(0)\bar{\zeta}_{j}(t)\bar{\eta}_{j}(t)}\right)^{\frac{1}{4}}\Biggr{\}}$$

$$\displaystyle\times\exp\int_{0}^{t}ds\Biggl{\{}\frac{1}{4}\sum_{j=1}^{N}\frac{(1-\bar{\zeta}_{j}\bar{\eta}_{j})(\dot{\bar{\zeta}}_{j}\bar{\eta}_{j}-\bar{\zeta}_{j}\dot{\bar{\eta}}_{j})}{\bar{\zeta}_{j}\bar{\eta}_{j}(1+\bar{\zeta}_{j}\bar{\eta}_{j})}$$

$$\displaystyle-i\mathcal{H}(\bar{\bm{\zeta}},\bar{\bm{\eta}},s)\Biggr{\}}.$$

$$\dot{\zeta }=i\Gamma \left(1-{\zeta }^2\right)-2i\zeta \left(h+J\frac{1-\zeta \eta }{1+\zeta \eta }\right),$$

$$\dot{\eta }=-i\Gamma \left(1-{\eta }^2\right)+2i\eta \left(h+J\frac{1-\zeta \eta }{1+\zeta \eta }\right).$$

$$\mathcal{L}(t|{\Omega }_{\mathrm{b}})=⟨{\Omega }_{\mathrm{b}}|e^{-i\hat{\mathcal{H}}t}|{\Omega }_{\mathrm{b}}⟩\sim \sum _{\nu }{A}_{\nu }{e}^{-Nf[{\bar{\zeta }}^{(\nu )},{\bar{\eta }}^{(\nu )}]},$$

$$f[\bar{\zeta },\bar{\eta }]=-\frac{1}{2}\log \frac{(1+{\zeta }_{\mathrm{b}}\bar{\eta }(0))(1+\bar{\zeta }(t){\eta }_{\mathrm{b}})}{(1+{\zeta }_{\mathrm{b}}{\eta }_{\mathrm{b}}{)}^2}$$

$$-i{\int }_0^{t}ds\left(\frac{\Gamma }{2}(\bar{\zeta }+\bar{\eta })+h+\frac{J}{2}\frac{1+2\bar{\zeta }\bar{\eta }-3{\bar{\zeta }}^2{\bar{\eta }}^2}{(1+\bar{\zeta }\bar{\eta }{)}^2}\right).$$

$$\bar{\zeta }(s)=\exp \left(-2is\frac{(1+C)h+(1-C)J}{1+C}\right),$$

$$\bar{\eta }(s)=C\exp \left(2is\frac{(1+C)h+(1-C)J}{1+C}\right).$$

$$C\exp \left(2it\frac{(1+C)h+(1-C)J}{1+C}\right)=1,$$

$$\tau =\frac{\pi }{h}.$$

Full text

I Introduction

Triggered by experiments using ultracold atomic systems, dynamics of a closed quantum many-body system has been one of the fascinating topics in condensed matter physics Polkovnikov:11 . In particular, the time evolution after a sudden change of the Hamiltonian has attracted a lot of attention as a basic setting of a problem on the out-of-equilibrium quantum state. One of the interesting phenomenon associated with this so-called quantum quench is the dynamical quantum phase transition (DQPT). While the equilibrium quantum phase transition is usually associated with a singularity of the ground-state energy in the axis of a parameter contained in the Hamiltonian, the DQPT involves a singularity in time. The present paper focuses on such a dynamical singularity appearing in the return probability to the initial state, which is directly related to the Loschmidt amplitude defined below Heyl:13 .

The phenomena of the DQPT are observed not only in the Loschmidt amplitude but also in the time average of local physical quantities such as order parameters. Although a certain correspondence is pointed out Zunkovic:16 , these two kinds of quantities are generally different. The local physical quantities represent the properties of the steady state in the long time limit after a quantum quench. They bring a clear physical consequence and are easy to access by experiments. The DQPT of them corresponds to a phase transition with the parameter in the Hamiltonian after a quantum quench. The Loschmidt amplitude, on the other hand, can be seen as an extension of the partition function on the imaginary axis corresponding to time. The DQPT here is defined as a singular behavior with time in the rate function of it, as an analogy with the thermodynamic phase transition accompanied by the singularity of the free energy as a function of the temperature. However, the Loschmidt amplitude involves delicate points in several aspects: physical meaning of the singularities, experimental implementations, and even technicalities for theoretical computations. Several recent works Heyl:13 ; Heyl:14 ; Heyl:15 ; Zunkovic:15 ; Zunkovic:16 ; Karrasch:13 ; Vajna:14 ; Sharma:15 ; Divakaran:16 ; Gambassi:11-1 ; Gambassi:11-2 ; Chiocchetta:15 ; Chiocchetta:16 ; Maraga:15 ; Smacchia:13 ; Gambassi:12 have been devoted to resolve the first delicate point based on statistical mechanical concepts such as renormalization group, symmetry breaking, universality, and scaling NishimoriBook . They have provided solid advances. For instance, the singularity has been tied with a behavior of the order parameter and entanglement production in systems with symmetry-broken phases Heyl:14 ; Jurcevic:16 . However, a general comprehension including the relation of the singularities to other local quantities with a generic initial state is still lacking. One of the origins of the difficulty in obtaining a general description lies, in our opinion, in the limitation on theoretical techniques to compute the Loschmidt amplitude. Most of theoretical works so far depend on the result of specific models being analytically tractable, and generic properties of the Loschmidt amplitude’s singularity have been speculated from the result. Hence, a more versatile computational method will be a great help to understand the Loschmidt amplitude.

Under this circumstance, here we propose a new theoretical framework for computing the Loschmidt amplitude for a generic spin system on the basis of a semiclassical computation. This can be regarded as a mean-field method and is expected to be exact in the infinite dimension, though it is still applicable as an approximation to a generic spin system in any dimension with an arbitrary state. The static approximation is often used with the mean-field method and is known to give a correct result for quantities in the equilibrium in the system with an infinite-range interaction SuzukiBook . However, the static approximation does not work for the computation of the Loschmidt amplitude. In this sense, our method goes beyond the static approximation and can be useful for computation of out-of-equilibrium quantities.

Our semiclassical method is essentially the same as the one used in Refs. Zunkovic:15 ; Sciolla:10 ; Sciolla:11 , but their analysis has been only on local physical quantities. This is presumably due to the lack of general prescriptions to compute the Loschmidt amplitude so far. The present work complements this point. The key difference of our method from the preceding studies Zunkovic:15 ; Sciolla:10 ; Sciolla:11 lies in the determination of the semiclassical path that follows the initial and final conditions properly. In our method, the range of dynamical variables is extended from real to complex numbers, and matching the semiclassical path with the boundary conditions is achieved in the complex space. This idea has been proposed in Refs. Klauder:79 ; Alscher:99 for single-spin systems and we extensionally apply it to many-spin systems. Accordingly, there emerge multiple solutions in the boundary value problem, and the solution that gives the largest return probability is selected. We find that it is this selection that gives the singularity in the Loschmidt amplitude. In this sense, the singularity of Loschmidt amplitude is very similar to an equilibrium phase transition.

Using the complex semiclassical method, in the present paper, we study the infinite-range transverse Ising model Hamiltonian with uniform coupling $J$ and longitudinal field $h$:

[TABLE]

where $\sigma_{i}^{\alpha}$ ($i=1,2,\dots,N$; $\alpha=x,z$) is the Pauli matrix and $N$ is the number of spins. For this model, we consider a quantum quench of the transverse field $\Gamma$ from $\Gamma_{\rm i}$ to $\Gamma_{\rm f}$ at $t=0$. As shown in Fig. 1, this system shows two different phases in equilibrium SuzukiBook and both inter- and intra-phase protocols of quench are examined.

The Loschmidt amplitude is defined by

[TABLE]

where the state $\left|{\psi}\right\rangle$ is chosen as the ground state of the Hamiltonian with $\Gamma=\Gamma_{\rm i}$. The Loschmidt amplitude is expected to exhibit a large deviation nature, and hence its rate function at $N\to\infty$ is the primary object of our analysis. The rate function is defined as

[TABLE]

Note that its real part, $f_{\rm r}=\Re{f}$, accounts for the return probability $P(t|\psi)=|\mathcal{L}\left(t|\psi\right)|^{2}$ as $2f_{\rm r}=-\frac{1}{N}\log P(t|\psi)$, while the imaginary part has no direct physical consequence.

The rest of the paper is organized as follows. In Sec. II, we describe the formulation and procedures needed to make the problem computationally tractable. In Sec. III, the analytical solutions computed from the invented method are shown and are compared to numerical experiments on finite size systems. Exact derivation of the rate function, available only on some specific parameters, is also given to justify the result. Section IV is devoted to discussion and summary. The relevance of the present work to experiments, quantum engineering, and computation is discussed there.

II Formulation

II.1 Spin coherent states and path integrals

We start from reviewing the path integral formulation for spin systems. An arbitrary state of a single spin is represented by a spin-coherent state as

[TABLE]

where $\left|{\uparrow}\right\rangle$ and $\left|{\downarrow}\right\rangle$ are the eigenstates of $\sigma^{z}$ with eigenvalues $+1$ and $-1$, respectively. Hereafter the gauge $b$ is fixed to be [math] and is disregarded, since it does not affect any physical consequences. As is well known, the average of spin variables $\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ over a spin-coherent state corresponds to three-dimensional polar representation as

[TABLE]

The spin-coherent state constitutes an overcomplete basis:

[TABLE]

where $I$ denotes the $2\times 2$ unit matrix. Note that two states with different $(\theta,\varphi)$ are not orthogonal in general.

We apply this spin-coherent state formulation to $N$-spin systems and write the variables as $(\bm{\theta},\bm{\varphi})=\{(\theta_{i},\varphi_{i})\}_{i=1}^{N}$. Using the spin-coherent states, we write any propagator with arbitrary time-dependent Hamiltonian $\hat{\mathcal{H}}(t)$ as $G(t|\Omega^{\prime},\Omega^{\prime\prime})\equiv\left\langle{\Omega^{\prime\prime}}\right|\mathcal{T}e^{-i\int^{t}_{0}ds\hat{\mathcal{H}}(s)}\left|{\Omega^{\prime}}\right\rangle$, where $\mathcal{T}$ is the time-ordering operator, and $\Omega^{\prime}=(\bm{\theta}^{\prime},\bm{\varphi}^{\prime})$ and $\Omega^{\prime\prime}=(\bm{\theta}^{\prime\prime},\bm{\varphi}^{\prime\prime})$ are initial and final states respectively. This propagator is rewritten in a path integral form as Klauder:79

[TABLE]

This is an integral over all possible paths of the variables $(\bm{\theta}(s),\bm{\varphi}(s))$. The action functional $S[\bm{\theta},\bm{\varphi}]$ is given by

[TABLE]

where the dot symbol denotes the time derivative and $\mathcal{H}\left(\bm{\theta},\bm{\varphi},s\right)=\langle\bm{\theta}(s),\bm{\varphi}(s)|\hat{\mathcal{H}}|\bm{\theta}(s),\bm{\varphi}(s)\rangle$.

II.2 Complex semiclassical analysis

The path integral formalism gives the exact result if we can perform the integration over all paths literally. However, this is difficult in general, and the semiclassical approximation is here employed.

II.2.1 Boundary value problem

The basic idea of the semiclassical method is to take into account only the dominant stationary paths among all the paths. The stationary condition in the action $S$ leads to the following equations of motion (EOMs):

[TABLE]

We naively expect that the solution of these EOMs, satisfying the boundary conditions $\left(\bm{\theta}(0),\bm{\varphi}(0)\right)=\Omega^{\prime}=\left(\bm{\theta}^{\prime},\bm{\varphi}^{\prime}\right)$ and $\left(\bm{\theta}(t),\bm{\varphi}(t)\right)=\Omega^{\prime\prime}=\left(\bm{\theta}^{\prime\prime},\bm{\varphi}^{\prime\prime}\right)$, is the desired semiclassical path. If there are multiple semiclassical paths, we give indices to them as $\{(\bar{\bm{\theta}}^{(\nu)},\bar{\bm{\varphi}}^{(\nu)})\}_{\nu}$ where the symbol $\bar{\cdot}$ represents a generic semiclassical path hereinafter. Semiclassical actions corresponding to those paths are defined as $S_{\rm cl}[\bar{\bm{\theta}}^{(\nu)},\bar{\bm{\varphi}}^{(\nu)}]=S[\bar{\bm{\theta}}^{(\nu)},\bar{\bm{\varphi}}^{(\nu)}]$. They give an approximation of the propagator as

[TABLE]

where $A_{\nu}$ denotes a possible amplitude factor.

Unfortunately, this procedure does not work in the present problem. The solution of Eq. (9) cannot satisfy, in general, both the boundary conditions $\left(\bm{\theta}(0),\bm{\varphi}(0)\right)=\Omega^{\prime}$ and $\left(\bm{\theta}(t),\bm{\varphi}(t)\right)=\Omega^{\prime\prime}$. For a given initial condition $\left(\bm{\theta}(0),\bm{\varphi}(0)\right)=\Omega^{\prime}$, the time evolution of the system is uniquely determined by the EOMs, and the final values $\left(\bm{\theta}(t),\bm{\varphi}(t)\right)$ do not necessarily coincide with the boundary one $\Omega^{\prime\prime}$. This is the reason why the Loschmidt amplitude has been difficult to be evaluated by the semiclassical or similar methods, though some exceptions are found when the semiclassical path is constant in time Zunkovic:15 ; Obuchi:12 ; Takahashi:13 . To overcome this problem, following the prescription in Refs. Klauder:79 ; Alscher:99 , we below introduce the so-called Wiener regularization term making the path integral well defined in the action, and deal with the unregularized action as the vanishing limit of the regularization term. This yields a different boundary condition.

II.2.2 Wiener regularization and modified boundary condition

By using the prescription by Klauder Klauder:79 , Alscher and Grabert demonstrated that the exact propagator can be computed in single-spin systems with arbitrary time-dependent magnetic fields Alscher:99 . Here we apply this to many-spin systems.

The Wiener regularization is defined as

[TABLE]

where $m$ represents a constant. Adding this term to the action, $S[\bm{\theta},\bm{\varphi}]\to S[\bm{\theta},\bm{\varphi}]+W[\bm{\theta},\bm{\varphi}]$, and taking the stationary condition, we obtain the modified semiclassical EOMs as

[TABLE]

Due to the regularization term, the higher-order derivatives appear in the EOMs and its general solution has more arbitrary constants, which naturally enables us to have a solution connecting to both the boundary values $\Omega^{\prime}$ and $\Omega^{\prime\prime}$. Meanwhile, the terms coming from the regularization introduce the imaginary number into the EOMs. Hence the corresponding semiclassical path becomes complex in general and loses a clear physical interpretation. Bloch sphere representation is not applicable to visualize the semiclassical path. From a formal correspondence, the Wiener regularization can be regarded as a kinetic energy of spins with a pure imaginary mass.

To recover the original action, we take the zero mass limit $m\to 0$. For small $m$, the time span $s\in[0,t]$ is divided into three characteristic regions Alscher:99 : $T_{1}=[0,m]$, $T_{\rm cl}=[m,t-m]$, and $T_{2}=[t-m,t]$. In $T_{\rm cl}$, the mass terms proportional to $m$ become irrelevant and the time evolution is essentially driven by the original unregularized EOMs. In $T_{1}$ and $T_{2}$, the trajectory is strongly hinged by the mass terms to match the boundary conditions. As a result, in the $m\to 0$ limit, we observe jumps at $s=0$ and $s=t$ from the boundary values to the edges of the semiclassical path in $T_{\rm cl}\to[0,t]$. These jumps give a condition for the values at the boundary $(\bar{\bm{\theta}}(0),\bar{\bm{\varphi}}(0))$ and $(\bar{\bm{\theta}}(t),\bar{\bm{\varphi}}(t))$, which has a simple explicit form:

[TABLE]

This condition implies that there can be multiple semiclassical paths to satisfy Eq. (13) and that they can be complex even in the $m\to 0$ limit. We note again that, for single-spin systems, it was shown in Ref. Alscher:99 that the solution of the unregularized EOMs (9) under the condition (13) gives the exact propagator.

II.2.3 Solving the boundary value problem

The boundary value problem becomes well-defined now and we can find solutions matching both the boundary values $\Omega^{\prime}$ and $\Omega^{\prime\prime}$ in a generic situation. A practical way for solving the problem is to employ the following variable transformation Alscher:99 :

[TABLE]

These variables are, if $(\theta_{j}(s),\varphi_{j}(s))$ are real, a stereographic representation of a point on the unit sphere projected from the south pole onto the equatorial plane. Hence we call them stereographic variables. The boundary condition is now written as

[TABLE]

and the remaining boundary values, $\zeta_{i}(t)$ and $\eta_{i}(0)$, are not specified. The spin variables in the Hamiltonian are converted to the stereographic variables through the relation

[TABLE]

and the semiclassical EOMs (9) are

[TABLE]

Using the solution of the EOMs, $(\bar{\zeta}_{j},\bar{\eta}_{j})$, we can write the semiclassical action as Alscher:99

[TABLE]

II.2.4 Spatially uniform solutions

A problem arises when we compute the semiclassical paths satisfying Eq. (15). We need to fix both the initial conditions on $\zeta_{i}$ and the final ones on $\eta_{i}$. The initial conditions on $\eta_{i}$ must be selected so as to meet the final conditions. This requires us to solve the EOMs many times, and results in a bottleneck of the present method to compute the propagator. This is because the computational cost for searching such an initial condition grows exponentially with the number of spins. Therefore, in practice, we need an assumption that reduces the degree of freedom, namely, the computational cost of searching the initial value of $\eta_{i}$.

In the present paper, we assume the spatial uniformity. Our Hamiltonian (1) has infinite-range interactions and the mean-field ansatz gives the exact result for static systems. Although it is not evident whether the spatial uniformity holds for dynamical systems, we examine this ansatz in the following. The boundary values of $(\zeta_{i}(s),\eta_{i}(s))$ are identical for all $i$’s, so that $(\zeta^{\prime}_{i},\eta^{\prime}_{i})=(\zeta^{\prime},\eta^{\prime}),$ and $(\zeta^{\prime\prime}_{i},\eta^{\prime\prime}_{i})=(\zeta^{\prime\prime},\eta^{\prime\prime})$. Then, only two functions, $\zeta(s)$ and $\eta(s)$, are sufficient to describe the dynamics, and the exhaustive search of $\eta(0)$ is now a reasonable task. Moreover, as far as the Loschmidt amplitude is concerned, the initial and final boundary values are common: $\zeta^{\prime}=\zeta^{\prime\prime}=\zeta_{\rm b}$ and $\eta^{\prime}=\eta^{\prime\prime}=\eta_{\rm b}$. Summarizing these particular conditions, we obtain the explicit formulas of the EOMs as

[TABLE]

For a given $t$, these EOMs are solved under the conditions $\zeta(0)=\zeta_{\rm b}$ and $\eta(t)=\eta_{\rm b}$. The other boundary values $\zeta(t)$ and $\eta(0)$ are not specified and are determined uniquely from the above conditions. We also note that the relation $\zeta(s)=\eta^{*}(s)$ does not necessarily hold in general.

The solutions $(\bar{\zeta}^{(\nu)}(s),\bar{\eta}^{(\nu)}(s))$ are not unique and we can represent the Loschmidt amplitude as

[TABLE]

where

[TABLE]

The time derivative terms are eliminated by performing the integration by parts or using the EOMs. We also note that the amplitude $A_{\nu}$ is not important to calculate the rate function in Eq. (3) at $N\to\infty$.

II.2.5 Dominant semiclassical paths and a heuristic search procedure

Equation (13) has a countably infinite number of solutions, and the EOMs do as well. Among those many semiclassical solutions, the one that makes the real part of $f[\bar{\zeta}^{(\nu)},\bar{\eta}^{(\nu)}]$ the smallest gives the rate function in Eq. (3). How can we find such a dominant solution? The exhaustive search of $\bar{\eta}(0)$ in the whole complex space is not plausible even under the spatial uniformity. To overcome the situation, we here give a heuristic procedure to obtain such a dominant path. Since the correct initial condition $\bar{\eta}(0)$ depends on the end time $t$, we hereafter use a notation $C(t)=\bar{\eta}(0;\bar{\eta}(t)=\eta_{\rm b})$. The basic idea of the heuristic is starting from a trivial solution at a specific time $t^{*}$ and extending it with changing the time $t$ from $t^{*}$ gradually.

The first trivial solution is obtained at $t^{*}=0$, where $C(0)=\eta_{\rm b}$. Then, for a small time step $\Delta t$, $C(\Delta t)$ is obtained as follows. We examine several values as the initial condition for $\bar{\eta}(s)$ around $\eta_{\rm b}$ and solve the EOMs. We select the best one for $C(\Delta t)$ that makes the final value $\bar{\eta}(s=\Delta t)$ closest to $\eta_{\rm b}$. For the next time step $t=2\Delta t$, we examine the values around $C(\Delta t)$ and repeat the same procedures, giving $C(2\Delta t)$. We repeat this procedure until we reach a desired end time $t$, yielding the sequence of the initial condition. We write this sequence as $C_{1}(t)$.

To obtain the second trivial solution, an important observation is that the dynamics is periodic at most of parameters Sciolla:11 . There exists a specific period $\tau$ and the order parameters at $t_{n}=t_{0}+n\tau$ are identical for $\forall{n}\in\mathbb{N}$. This implies that at $t^{*}=\tau$ the final condition $\bar{\eta}(s=t^{*})=\eta_{\rm b}$ is realized by having the initial condition $\bar{\eta}(0)=\eta_{\rm b}$, yielding $C(t^{*}=\tau)=\eta_{\rm b}$. Extending $C(t)$ back from $t=\tau$ to $t=0$ based on the same procedure for $C_{1}(t)$, we get another sequence of the initial condition, and write it as $C_{2}(t)$. In Fig. 2, a schematic picture of this heuristic is given.

The question is whether these two sequences of initial conditions, $C_{1}(t)$ and $C_{2}(t)$, are identical or not. If they are different, they give two different semiclassical paths. In such a situation, there should be a switch between two paths at a certain critical time $t_{\rm c}$ in the period $[0,\tau]$, that yields a singularity of the Loschmidt amplitude. Meanwhile, if they are identical, only one dominant semiclassical path exists and is analytic with respect to $t$.

For longer time $t>\tau$, we repeat the above procedure. For the next period $[\tau,2\tau]$, $C_{3}(t)$ is obtained by extending $C(t)$ from $t=\tau$ to $2\tau$ with the trivial value $C(\tau)=\eta_{\rm b}$, and $C_{4}(s)$ is given by an extension from $t=2\tau$ to $\tau$ with $C(2\tau)=\eta_{\rm b}$. We note that by construction $C_{2}(s)$ and $C_{3}(s)$ are continuously connected. The solutions for the whole time axis are obtained along this way.

We adopt the above scenario to search the solution. This may give a wrong result in general, but, as far as we have investigated, the result shows a good agreement with numerical experiments as we see in the following. Our heuristic procedure is constructed under the assumption that the system shows a periodic behavior and only one transition at most in one cycle. As long as this assumption is true, our heuristic can find the correct dominant path. For more general cases, e.g. spin glasses without periodicity Obuchi:12 ; Takahashi:13 , other heuristics should be tailored. Investigation of such cases is beyond the scope of this paper and will be an interesting future work.

III Result

We present the results of our semiclassical computation. We study two cases: quenches from $\Gamma_{\rm i}=\infty$ (Sec. III.1) and quenches from $\Gamma_{\rm i}=0$ (Sec. III.2). The first case is the quench from $\Gamma_{\rm i}=\infty$ to a finite value $\Gamma_{\rm f}<\infty$, where the boundary condition is the ground state at $\Gamma_{\rm i}=\infty$, namely, $\left|{\Omega_{\rm b}}\right\rangle=\otimes_{i}\left|{\rightarrow}\right\rangle_{i}$ with $\left|{\rightarrow}\right\rangle_{i}$ being the eigenstate of $\sigma_{i}^{x}$ for eigenvalue $+1$. The other case is the opposite quench, from $\Gamma_{\rm i}=0$ to $\Gamma_{\rm f}>0$. We set $h\geq 0+$ and thus $\left|{\Omega_{\rm b}}\right\rangle=\otimes_{i}\left|{\uparrow}\right\rangle_{i}$. Since exact calculation is possible for a quench from $\Gamma_{\rm f}=\infty$ to $\Gamma_{\rm i}=0$, we show its result in Sec. III.1 as well. We also show the results of numerical studies in Sec. III.3 to confirm that the complex semiclassical analysis gives a reasonable result.

III.1 Quench from $\Gamma_{\rm i}=\infty$

In this case, the boundary condition is given by $(\theta^{\prime},\varphi^{\prime})=(\theta^{\prime\prime},\varphi^{\prime\prime})=(\pi/2,0)$, that is $(\zeta_{\rm b},\eta_{\rm b})=(1,1)$. With this boundary condition, if $h=0$, the state does not evolve and the semiclassical path is written as $\bar{\zeta}(s)=\bar{\eta}(s)=1$ for $\forall{s}$. Hence, we consider the case $h>0$ where, as we show below, a finite periodicity $0<\tau<\infty$ is present. In fact, we see several patterns of the rate function and DQPT as well. We obtain the corresponding phase diagram.

III.1.1 A solvable case: $\Gamma_{\rm f}=0$

We first investigate the quench to $\Gamma_{\rm f}=0$. In this case, the state is evolved under the classical Ising Hamiltonian and an analytical solution of Eq. (22) is available. We solve the equation under the conditions $\zeta(0)=1$ and $\eta(t)=1$. Putting the initial condition as $(\zeta(0),\eta(0))=(1,C)$, we get the explicit solution of the dynamics as

[TABLE]

Then the condition $\eta(t)=1$ gives us

[TABLE]

which yields $C(t)$.

This example clearly shows the presence of multiple paths satisfying the boundary condition. As declared in Sec. II.2.5, we investigate two paths associated with the initial conditions $C_{1}(t)$ and $C_{2}(t)$, each of which is continuously extended from $C(0)=\eta_{\rm b}=1$ and from $C(\tau)=\eta_{\rm b}=1$ respectively, where $\tau$ is the period of the dynamics. The period $\tau$ can be obtained by putting $C=1$ and $t=\tau$ in the solution (25) as

[TABLE]

The solutions of Eq. (26) connecting to $C(0)=1$ and $C(\tau)=1$, $C_{1}(t)$ and $C_{2}(t)$ respectively, are shown in Fig. 3 for $h/J=0.1$. As a reference, the solutions in the next period $[\tau,2\tau]$, $C_{3}$ and $C_{4}$, are also displayed.

Given a time $t$, the rate function $f_{k}(t)$ corresponding to the initial condition $C_{k}(t)$ is evaluated by inserting Eq. (25) with $C=C_{k}(t)$ into Eq. (24) and performing the integration with respect to $s$. The result is shown in the right panel of Fig. 3. This exhibits the DQPT at $t=\tau/2$ where the switch from $f_{1}$ to $f_{2}$ occurs. Similarly, the switch from $f_{3}$ to $f_{4}$ occurs at $t=\frac{3}{2}\tau$, showing another DQPT.

Apart from the analytical solution of the EOMs, the rate function itself can be computed exactly for the case $\Gamma_{\rm f}=0$. Dashed lines in the right panel of Fig. 3 represent the result. In terms of the total spin operator $\hat{S}_{z}=\frac{1}{2}\sum_{i}\sigma_{i}^{z}$, our Hamiltonian after the quench is written as

[TABLE]

The eigenvalue of this Hamiltonian is characterized by that of $\hat{S}_{z}$ denoted by $M$ taking the value $\frac{N}{2}-k$ with $k=0,1,\dots,N$. For a given $M$, there are $\binom{N}{k}$ degenerate states. Let us define a normalized vector $\left|{\frac{N}{2}-k}\right\rangle$ in this subspace, which is the equal-weight sum of the $\binom{N}{k}$ basis vectors. Using this basis, we can write the initial state as

[TABLE]

Applying the time-evolution operator $e^{-i\hat{\mathcal{H}}t}$ only gives a phase factor for each term. The Loschmidt amplitude is written as

[TABLE]

Using an approximation valid for $N\gg 1$

[TABLE]

we write the amplitude as

[TABLE]

For $N\gg 1$ the range of the sum can be safely extended from $k=-\infty$ to $\infty$ to yield

[TABLE]

where the Poisson summation formula is used in the first line NishimoriBook . We then define $n$ that minimizes $f_{\rm r}(t;n)$ as $n^{*}(t)=\mathop{\rm arg~{}min}\limits_{n\in\mathbb{Z}}f_{\rm r}(t;n)$, where

[TABLE]

The contribution from $n=n^{*}$ dominates the sum in Eq. (33) and we obtain the real part of the rate function as $f_{\rm r}(t)=f_{\rm r}(t,n^{*}(t))$. $f_{\rm r}(t)$ exhibits singularities because $n^{*}(t)$ changes discretely as $t$ grows. Thus the transition time $t_{\rm c}$ is obtained by equating two neighboring values $f_{\rm r}(t;n)=f_{\rm r}(t;n+1)$ as

[TABLE]

Hence the period is given by $\tau=\pi/h$. The branches of $f_{\rm r}(t;n)$ with $n=0,1,2$ are shown in the right panel of Fig. 3, which exhibits the perfect agreement with $f_{k}(t)$ evaluated by the integration in Eq. (24) with the solution of EOMs (25) and $C_{k}(t)$.

Two noteworthy consequences are provided by this analytical solution. One is that the DQPT always exists for any $h>0$, while it does not for $h=0$. Some earlier works have pointed out that a DQPT appears when quench crosses an equilibrium quantum phase transition Heyl:13 ; Karrasch:13 . The present results reveal the existence of the opposite situation. The other is that the real part of the rate function $f_{\rm r}(t)=f_{\rm r}(t;n^{*}(t))$ shrinks by the speed of $O(t^{-2})$ as $t$ grows and finally vanishes in the limit $t\to\infty$. The vanishing rate function may be thought to imply $|\mathcal{L}(t)|\to 1$, but this is not the case because of the presence of the factor $1/\sqrt{1-iJt}$ in Eq. (33). The modulus of this factor decreases as $t$ grows, so that $|\mathcal{L}(t)|$ goes to zero. This implies that there exists a crossover time $t_{\times}$ determined by comparing the $O(1)$ factor and the exponentially scaling one $e^{-Nf}$. For $t>t_{\times}$ the $O(1)$ factor dominates the Loschmidt amplitude. However, the crossover time $t_{\times}$ is expected to a unbounded increasing function of $N$. Hence in the large size limit our computation of the rate function is meaningful in the whole time region.

III.1.2 General $\Gamma_{\rm f}>0$

Let us proceed to general final values $\Gamma_{\rm f}>0$. The analytical solution of the EOMs is not available in this case. Hence we numerically search the initial conditions $C_{1}(t)$ and $C_{2}(t)$, and evaluate the corresponding rate functions $f_{1}(t)$ and $f_{2}(t)$.

Our heuristic procedure starts from evaluating the period $\tau$ of the dynamics. For this purpose, we run the numerical simulation of the EOMs (22) using the naive initial condition, $\zeta(0)=\zeta_{\rm b}=1$ and $\eta(0)=\eta_{\rm b}=1$. We employ the Runge-Kutta method of the fourth order. As an example, we show the result for the case with $h/J=0.1$ and $\Gamma_{\rm f}/J=0.6$ in the left panel of Fig. 4.

The period $J\tau\approx 5.8$ is easily read from this panel. We again stress that this dynamics with the naive condition $\eta(0)=\eta_{\rm b}=1$ does not satisfy the boundary condition (13) for a generic end time $t$. Given an end time $t$, we need to estimate the appropriate initial condition $\eta(0)=C(t)$, and then compute the path $(\bar{\zeta}(s),\bar{\eta}(s))$. As a result, the semiclassical paths satisfying Eq. (13) are very different from the naive ones. Putting the end time as $Jt=J\tau\approx 5.8$, we plot the real parts of such paths in the center panel of the same figure. As explained in Sec. II.2.5, we have two different sequences of the initial conditions, yielding two different paths $(\bar{\zeta}_{1}(s),\bar{\eta}_{1}(s))$ and $(\bar{\zeta}_{2}(s),\bar{\eta}_{2}(s))$. Both paths satisfy $\bar{\zeta}(0)=\zeta_{\rm b}=1$ and $\bar{\eta(t)}=\eta_{\rm b}=1$ as they should. The real parts of the corresponding two rate functions are plotted in the right panel. The smaller branch at each time corresponds to the true rate function, leading to the DQPT at $Jt_{\rm c}\approx 4.42$ as a switch from $f_{1}$ to $f_{2}$. Note that this panel is plotted against the end time $t$ while the center one is plotted against the dummy time $s$, given the end time $t=\tau$.

The DQPT observed here has the nature of the first order transition, in a sense that the first order time derivative of the rate function jumps at the transition time. By examining the several parameters, we have realized that this first order nature tends to be stronger as $\Gamma_{\rm f}$ increases, but suddenly vanishes at a certain critical value $\Gamma_{\rm fc}(h)$. For $\Gamma>\Gamma_{\rm fc}(h)$, the curve of the rate function has a smooth peak without singularity. In Fig. 5, we plot $C_{1},C_{2},f_{1},$ and $f_{2}$ for $h/J=0.1$ with slightly different two values of $\Gamma_{\rm f}$, $\Gamma_{\rm f}/J=1.5$ and $1.6$.

They clearly show that the critical value $\Gamma_{\rm fc}(h)$ is present between these two values of $\Gamma_{\rm f}$. In the same way, computing the rate function in a range of $h$ and $\Gamma_{\rm f}$, we draw a phase diagram in the case of quench from $\Gamma_{\rm i}=\infty$ in Fig. 6.

The phase boundary approaches to the equilibrium transition point $\Gamma_{\rm c}=J$ in the limit $h\to 0$. This is reasonable because the period of the dynamics $\tau$ diverges as $h\to 0$ at $\Gamma_{\rm f}<\Gamma_{\rm c}$, and DQPTs do not exist according to the present scenario.

III.2 Quench from $\Gamma_{\rm i}=0$

We next study the opposite quench from $\Gamma_{\rm i}=0$. The boundary condition is now given by $\zeta_{\rm b}=\eta_{\rm b}=0$.

As in the previous case, the numerical search of $C_{1}(t)$ and $C_{2}(t)$ brings the behavior of the rate function and the DQPT in this case. However the results are rather different. In the previous case, the DQPT was the first order like and there was a prominent cusp in a period $[0,\tau]$. When going across the DQPT boundary, the cusp turned into a smooth peak and the bifurcation or merge of the two initial conditions $C_{1}(t)$ and $C_{2}(t)$ occurs in the middle of the period $[0,\tau]$. For the quench from $\Gamma_{\rm i}=0$, however, this is not the case and the DQPT emerges in a more delicate form.

Figure 7 is the plots of $C_{1},C_{2},f_{1},$ and $f_{2}$ for $h/J=0.1$ with slightly different two values of $\Gamma_{\rm f}$, $\Gamma_{\rm f}/J=0.6$ and $0.7$.

This figure demonstrates that the bifurcation of the two initial conditions $C_{1}(t)$ and $C_{2}(t)$ occurs around $t\approx\tau$ in a rather continuous manner. As a result, discriminating the two branches of the solution is harder than the quench from $\Gamma_{\rm i}=\infty$. This tendency holds for the range of $h$ and $\Gamma_{\rm f}$ we have searched, which requires us to conduct a more precise numerics to obtain the phase diagram. Moreover, as we see from the bottom panels ($\Gamma_{\rm f}/J=0.7$), $C_{2}(t)$ tends to show a rather singular behavior: a smooth curve suddenly changes into a plateau as $t$ decreases and finally it vanishes for small $t$. Although we cannot completely reject a possibility that these behaviors are caused by certain numerical errors, we have carefully checked and confirmed that the shown $C_{2}(t)$ satisfies the boundary condition in a good precision for the intermediate and large $t$ region, and no branches exists continuously connected to the plateau for small $t$. Hence, these singular behaviors are expected to be true. Fortunately, they are irrelevant for locating the DQPT point since the DQPT occurs at larger $t$ where no pathological behavior appears. The resultant phase diagram for $\Gamma_{\rm i}=0$ is given in Fig. 8.

The phase boundary $\Gamma_{\rm fc}(h)$ approaches to $\Gamma_{\rm d}/J=1/2$ in the limit $h\to 0+$. This is reasonable because $\Gamma_{\rm d}$ is the dynamical transition point of an order parameter, $m_{z}=\left\langle{\Psi(t)}\right|\sigma^{z}\left|{\Psi(t)}\right\rangle$ where $\left|{\Psi(t)}\right\rangle=e^{-it\hat{\mathcal{H}}}\left|{\Omega_{\rm b}}\right\rangle$, with this particular choice of $\Gamma_{\rm i}$ Sciolla:11 . Upon approaching to $\Gamma_{\rm d}$, the period $\tau$ of the dynamics diverges and DQPTs should vanish. Note that this dynamical value $\Gamma_{\rm d}$ does not have any meaning for the equilibrium transition. This is in contrast to the quench from $\Gamma_{\rm i}=\infty$ where the equilibrium transition point $\Gamma_{\rm c}$ works as the DQPT transition point at $h=0$.

Unlike in the $\Gamma_{\rm i}=\infty$ case, the dynamics does not stop even at $h=0$. This enables us to see an interesting behavior of the Loschmidt amplitude at $h=0$ and $\Gamma_{\rm f}=\Gamma_{\rm d}$. This point is on the separatrix in the phase space and the order parameter monotonically decreases as $t$ grows. No periodicity exists (or $\tau=\infty$). Hence, we only examine the first sequence of the initial condition for $\eta(s)$ and $C_{1}(t)$, and compute the corresponding rate function. The result is shown in Fig. 9.

This figure shows that the rate function asymptotically vanishes as $t\to\infty$, but this does not necessarily imply $|\mathcal{L}(t)|\to 1$ as pointed out at the end of Sec. III.1.1.

III.3 Comparison with numerical experiments

To validate our semiclassical computations, we here show the results of numerical experiments and compare them with the semiclassical results for several parameters. Our Hamiltonian (1) commutes with the squared total spin operator $\hat{\bm{S}}^{2}=\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}$. For both the quenches from $\Gamma_{\rm i}=0$ and $\infty$, the initial state is in the subspace of the total spin $S=N/2$. Hence the state of our system preserves this total spin and we may consider the time-dependent state inside this subspace. In the basis of the eigenvalues of $\hat{S}_{z}$, our Hamiltonian is represented in a tridiagonal matrix form and we can easily evaluate the time evolution of the state $\left|{\Psi(t)}\right\rangle=e^{-it\hat{\mathcal{H}}}\left|{\Omega_{\rm b}}\right\rangle$ by the LU decomposition. The dimension of the subspace is $N+1$ and we can treat fairly large size systems. However, the computation requires us to take a lot of sums of complex numbers and the numerical precision tends to be degraded as $N$ becomes large. This computational difficulty sensitively depends on the parameters and below the simulated system sizes are adaptively changed for this reason.

Figure 10 is the plots of the rate functions for the quench from $\Gamma_{\rm i}=\infty$.

The results of numerical simulation show a good agreement with the theoretical curve denoted by the solid black line, both below and above the transition point $\Gamma_{\rm fc}(h/J=0.1)/J\approx 1.53$. This justifies our semiclassical computation. The upper left panel in Fig. 10 for $\Gamma_{\rm f}=0$ and $h/J=0.1$ is compared with the result in Fig. 3 where the period is given by $J\tau=\frac{\pi J}{h}\approx 31.4$. We see the consistent agreement between the numerical and semiclassical computations. The deviation for the whole time and the oscillating behavior at large $t$ are considered to be due to the finite size effect.

Figure 11 represents the result of a quench from $\Gamma_{\rm i}=0$ at $h=0+$.

Again, the numerical results show a fairly good agreement with the semiclassical curve. At the separatrix, $\Gamma_{\rm f}/J=1/2$, the monotonic decay of the rate function after a single peak is well reproduced by the numerics, validating our semiclassical computation even at a special point of the dynamics.

IV Discussion and summary

In this paper, we have invented a computational method for the Loschmidt amplitude based on the complex semiclassical approach, and applied it to the transverse field Ising model with a symmetry breaking field in the infinite dimension. Two quantum quenches, from zero and infinite transverse fields, have been examined. From the behavior of the rate function, the presence or absence of the DQPTs have been captured. The phase diagrams have been mapped out in the plane of final transverse field and symmetry breaking field. These results have been examined by numerical simulations independently that solve the Schrödinger equation literally, which fully supports our semiclassical computations.

Although our computational method has succeeded in unveiling several properties of the Loschmidt amplitude, its physical implications are still unclear. Žunkovič et al. have pointed a connection between the Loschmidt amplitude and the order parameter in the steady state long after the quench Zunkovic:16 . However we have not found such a connection as far as the quench from $\Gamma_{\rm i}=\infty$ is concerned. Therefore the presented result might add a further mystery on the DQPT. Disentangling DQPTs of the Loschmidt amplitude and an order parameter may open a new comprehension on quantum dynamics.

An experimental observation of a DQPT is a fascinating topic. A very recent work Flaschner:16 has actually observed DQPTs using a certain topological nature of the singularity Vajna:15 ; Budich:16 ; Sharma:16 . Unfortunately, this is possible only in non-interacting systems and its generalization to interacting systems is unclear. Although there are some other experiments Zhang:09 ; Zangara:16 observing the Loschmidt amplitude, their methods rely on the smallness of the system or certain locality of the phenomena. The application of their methods to global phenomena in many-spin systems is again nontrivial. Our model, the Ising model with long range interactions, itself can be realized in a trapped ion system Islam:13 . Another recent experiment on this system has observed nontrivial cusps in the probability to return to the ground-state manifold, giving a clear evidence of the DQPT Heyl:14 ; Jurcevic:16 . Their setup corresponds to $\Gamma_{i}=0$ and $h=0+$ in the present paper, and we expect that further nontrivial results can be obtained in other setups according to our findings. Such additional experiments are encouraged.

A more direct application of our method might be found in quantum engineering or computing. In those disciplines, it is an important problem to estimate the probability achieving a desired state in certain quantum processes. For example in quantum annealing QA1 ; QA2 ; QA3 , the probability to find the ground state is an important object to be calculated. Using techniques from the spin glass theory Obuchi:12 ; Takahashi:13 combined with the present method, its typical value might be evaluated. This will provide a theoretical challenge for both quantum mechanics and random spin systems.

Acknowledgements

The authors thank Makoto Negoro for discussions about experimental relevance. This work was supported by KAKENHI Nos. 26870185 (TO), 26400402 (SS), and 26400385 (KT).