Nonadiabatic dynamics of the excited states for the Lipkin-Meshkov-Glick model

TL;DR

This paper studies how excited state quantum phase transitions influence the nonadiabatic dynamics in the Lipkin-Meshkov-Glick model, revealing state-dependent scaling behaviors and validating mean-field approaches in the thermodynamic limit.

Contribution

It demonstrates the impact of excited state quantum phase transitions on dynamics and confirms mean-field methods accurately capture scaling near these transitions.

Findings

Scaling properties depend on energy eigenstates.

Mean-field dynamics match exact results in the thermodynamic limit.

Excited state quantum phase transitions significantly affect nonadiabatic behavior.

Abstract

We theoretically investigate the impact of the excited state quantum phase transition on the adiabatic dynamics for the Lipkin-Meshkov-Glick model. Using a time dependent protocol, we continuously change a model parameter and then discuss the scaling properties of the system especially close to the excited state quantum phase transition where we find that these depend on the energy eigenstate. On top, we show that the mean-field dynamics with the time dependent protocol gives the correct scaling and expectation values in the thermodynamic limit even for the excited states.

Figures

Equations

$$\hat{H}=-h{\hat{J}}_z-\frac{{\gamma }_x}{N}{\hat{J}}_x^2,$$

$${\hat{P}}_{\pm }=\frac{1}{2}\left[1\pm \exp \left(i\pi {\hat{J}}_z\right)\exp \left(i\pi N/2\right)\right].$$

$$⟨j,m^{\prime }|\hat{H}|j,m⟩$$

$$=\left[-hm-\frac{{\gamma }_x}{2N}\left(k_{-}(j,m{)}^2+k_{+}(j,m{)}^2\right)\right]{\delta }_{m^{\prime },m}$$

$$-\frac{{\gamma }_x}{2N}{k}_{+}(j,m)k_{+}(j,m+1){\delta }_{m^{\prime },m+2}$$

$$-\frac{{\gamma }_x}{2N}{k}_{-}(j,m)k_{-}(j,m-1){\delta }_{m^{\prime },m-2}$$

$${\hat{H}}_{+}={\hat{P}}_{+}^{†}\hat{H}{\hat{P}}_{+},$$

$${\hat{H}}_{+}|i⟩=E_i^{+}|i⟩.$$

$${\hat{J}}_z|m⟩=m|m⟩.$$

$$E_{i+1}-E_i=a_i\cdot {\left|{\gamma }_x-{\gamma }_x^{cr,i}\right|}^{-c_i}$$

$${\gamma }_x(t)=\frac{t}{Q},$$

$$i\dot{|\psi (t)⟩}={\hat{H}}_{+}(t)|\psi (t)⟩,|\psi (t=t_{\mathrm{in}})⟩=|i(t_{\mathrm{in}})⟩.$$

$$E_i^{res}(t)\equiv ⟨\psi (t)|\hat{H}(t)|\psi (t)⟩-⟨i(t)|\hat{H}(t)|i(t)⟩,$$

$$|\psi (t)⟩=\sum _i{\alpha }_i(t)|i(t)⟩,$$

$$F(t)={\left|⟨\psi (t)|i(t)⟩\right|}^2.$$

$$\sigma (t)\equiv \frac{1}{N}\sqrt{\sum _j{E}_j^{res}(t{)}^2{\left|{\alpha }_j(t)\right|}^2},$$

$$\sigma (Q)\sim Q^{\alpha },$$

$${|\theta ,\phi ⟩}_t$$

$$\times \exp \left[\exp (i\varphi )\tan (\theta /2){\hat{J}}_{-}\right]|N/2⟩$$

$${|\theta ,\phi ⟩}_t$$

$$\times \sum _{m=-N/2}^{N/2}\sqrt{\left(\genfrac{}{}{0ex}{}{N}{m+N/2}\right)}\cdot \sin {\left(\frac{\theta }{2}\right)}^{\frac{N}{2}+m}$$

$$\times \cos {\left(\frac{\theta }{2}\right)}^{\frac{N}{2}-m}\cdot \exp \left[-i\left(\frac{N}{2}+m\right)\cdot \phi \right]|m⟩.$$

$$|\Psi (t)⟩=\int d\theta d\phi \beta (\theta ,\phi ){|\theta ,\phi ⟩}_t,$$

$$\beta (\theta ,\phi )=⟨\Psi (t)|\theta (t),\phi (t)⟩.$$

$${\dot{J}}_x(t)$$

$${\dot{J}}_y(t)$$

$${\dot{J}}_z(t)$$

$${\sigma }_M(Q)\sim Q^{\alpha }.$$

Full text

I Introduction

Universality and feature analogy connect different branches of physics. One of the most famous examples are phase transitions which span from cosmology to the quantum world Del Campo and Zurek (2014); Zurek (1985); Hohenberg and Halperin (1977). Here one of the intriguing questions are scaling properties around the phase transition in both equilibrium and non-equilibrium setups. One way to create the latter is a continuous quench from one state to another Polkovnikov et al. (2011); Yukalov (2011). The connection between both can be characterized by the Kibble-Zurek mechanism, which connects the quench-induced defect formation with the universal scaling exponents in equilibrium such as relaxation time and the healing length Kibble (1980); Zurek (1996).

A quantum phase transition (QPT) is a paradigm at the small-scale domain, its famous property is the non-analyticity of the ground-state energy and ground-state wave-function at some point in the system configuration space. The validity of the Kibble-Zurek mechanism has been applied to QPTs in the Ising-like models with a finite dimension Zurek et al. (2005); Dziarmaga (2005); Acevedo et al. (2014). Later it was shown that the scaling connection still holds even in zero-dimensional models like the Dicke Dicke (1954), where the length-scale is absent Hwang et al. (2015); Acevedo et al. (2015). Furthermore, its connection to Landau-Zener transitions has been shown Landau (1932); Zener (1932); Del Campo and Zurek (2014); Damski (2005). The scaling connection is relevant e.g. for quantum computation algorithms in systems with a QPT Mostame et al. (2010, 2007). Note, that the Kibble-Zurek mechanism has been verified in various experimental setups e.g. using quantum atomic gases Weller et al. (2008); Corman et al. (2014); Braun et al. (2015); Yukalov et al. (2015), superfluids Bäuerle et al. (1996) or Coulomb Crystals Mielenz et al. (2013).

Here we extend the study of scaling properties toward excited state quantum phase transition (ESQPT) Caprio et al. (2008); Pérez-Fernández et al. (2011a); Brandes (2013); Stránský et al. (2014, 2015). Similar to the ground state QPT, the ESQPT is defined by a discontinuity in the density of states in the excited part of the spectrum, it is imprinted in different observables and was experimentally found e.g. in the bending spectra of different molecular systems Winnewisser et al. (2005); Larese and Iachello (2011); Larese et al. (2013).

For our investigation we consider the so-called Lipkin-Meshkov-Glick (LMG) model, which describes the dynamics of $N$ two-level systems with an internal all-to-all coupling. This model was first introduced in the context of nuclear physics to check the validity of different predictions Lipkin et al. (1965); Meshkov and Glick (1965); Glick et al. (1965). Recently, it has attracted much attention again in the context of cold atoms and many-body physics, especially as a rather simple model with an experimental realization for studying QPTs Zibold et al. (2010); Gross et al. (2010); Pan and Draayer (1999); Ribeiro et al. (2008); Richerme et al. (2014).

The QPT-induced scaling properties, the excited state impact to the wave-package dynamics or the transitionless driving in the time-dependent LMG model in presence of a linear or a sudden quench has been theoretically investigated in Ref. Caneva et al. (2008); Solinas et al. (2008); Pérez-Fernández et al. (2011b); Campbell et al. (2015); Puebla and Relano (2015); Santos and Pérez-Bernal (2015); Santos et al. (2016), but with the ground state as the initial condition. The ground state is a common choice, especially in the context of the quantum computation and annealing problems Puri et al. (2016); Schaller et al. (2006). In contrast, here we choose several excited states as the initial condition to probe the scaling effects induced by the ESQPT and find that the scaling parameter decreases and converges to a constant value if the energy of the chosen initial excited state increases. Note, that our study has experimental relevance, as the excited states can be prepared with coherent states which are commonly used in the existing experimental setup of the LMG model Zibold et al. (2010); Gross et al. (2010).

Our work is organized as follows. In Sec. II we introduce the model, in Sec. III we review its spectral characteristics and investigate the gap closing scaling for the excited states. In Sec. IV we numerically investigate the scaling properties due to the ESQPT and in Sec. V the state dynamics for different speeds of the quench protocol by solving the time-dependent Schrödinger equation. In Sec. VI we show that the results from the previous chapters can be obtained form the corresponding mean-field model. We summarize our findings in Sec.VII.

II Model

The anisotropic LMG Hamiltonian reads Orús et al. (2008)

[TABLE]

where $\hat{J}_{\eta}=\frac{1}{2}\sum_{i=1}^{N}\hat{\sigma}_{\eta}^{(i)}$, $\eta=\{x,y,z\}$, are collective angular momentum operators and $\hat{\sigma}_{x,y,z}$ are the common Pauli matrices. The magnetic field $h$ defines the level splitting in the non-interacting case, whereas $\gamma_{x}$ describes the coupling strength between the $N$ two-level systems. Note, that here the angular momentum $j$ is directly connected to the system size via $N=2j$, this suggests to use the angular momentum basis of the $\hat{J_{z}}$ operator $\left|j,m\right>$ for the matrix representation of the Hamiltonian Eq. (1).

The Hamiltonian Eq. (1) preserves the length of the total spin $\hat{J}^{2}=\hat{J}_{x}^{2}+\hat{J}_{y}^{2}+\hat{J}_{z}^{2}$ Leyvraz and Heiss (2005) and has a parity symmetry which divides $\hat{H}$ into two different independent energetic subspaces, one even ’+’ and one odd ’–’ Lerma-H and Dukelsky (2014a); Sorokin et al. (2014). The corresponding parity operator reads

[TABLE]

Note, the parity conservation is directly imprinted in the matrix representation of the Hamiltonian Eq. (1)

[TABLE]

as a decomposition into two decoupled regions with even or odd state $\left|j,m\right>$. Here $\hat{J}_{\pm}\left|j,m\right>=k_{\pm}(j,m)\left|j,m\pm 1\right>$ with $\hat{J}_{\pm}=\hat{J}_{x}\pm i\hat{J}_{y}$ is used.

As parity is conserved, we will restrict ourselves only to the even ’+’ parity and will perform all calculations with the even Hamiltonian

[TABLE]

which contains the ground state of the system. In the following, we will denote the eigenstates of the Hamiltonian $\hat{H}_{+}$ by $\left|i\right>$ for fixed $\gamma_{x}$ values, i.e.

[TABLE]

Note that the spectra of both even and odd Hamiltonians become identical in the $N\to\infty$ limit. The lower part of the even-parity spectrum is shown in Fig. 1 as a function of the coupling strength $\gamma_{x}$.

III Spectral Properties

A lot of intriguing physical properties are incorporated in Eq. (1). Most prominently, a second order quantum phase transition (QPT) appears in the thermodynamic limit $N\to\infty$ at $\gamma_{x}/h=\gamma_{x}^{cr,1}/h\equiv 1$ Gilmore and Feng (1978); Gilmore et al. (1978); Ribeiro et al. (2008), which is visible in the energy spectrum. There, the ground state energy becomes non-analytic and the gap to the first excited state closes as a function of $\gamma_{x}/h$, see the two blue (solid and dashed) lines in Fig. 1 (left panel). The energy distance between the ground and the first excited state is shown in Fig. 2 (top panel) by blue colored lines for different particle numbers $N$. For $N\to\infty$ the ground state gap closes and the position of the minimum of the $i=1$ curves converges to $\gamma_{x}^{cr,1}/h=1$. Note that the energy axis in Fig. 2 is not scaled with $1/N$ , unlike the energy axis in Fig. 1. The energy spectrum was determined by a numerical diagonalization of $\hat{H}_{+}$ in Eq. (4) represented in the $\hat{J}_{z}$ basis $\{\left|m\right>\}$, i.e.

[TABLE]

In the thermodynamic limit the excited spectrum possesses an additional phase transition - the excited state quantum phase transition (ESQPT) Leyvraz and Heiss (2005); Caprio et al. (2008); Pérez-Fernández et al. (2009); Bastidas et al. (2014); Engelhardt et al. (2015). In the case of the LMG model, the ESQPT is visible as a logarithmic divergence of the spectral level density at the energy $E=-0.5/(h\cdot N)$ of the so called separatrix, which existence is typical for a ESQPT Cejnar and Stránský (2008); Cejnar and Jolie (2009), for a fixed $\gamma_{x}>\gamma_{x}^{cr,1}$ value Ribeiro et al. (2007). In the right panel of Fig. 1 we plot the numerically determined level density for two different $\gamma_{x}$ values. Whereas for $\gamma_{x}<\gamma_{x}^{cr,1}$ the density of states is approximately flat (dashed curve), for $\gamma_{x}>\gamma_{x}^{cr,1}$ the density divergence is clearly visible (solid curve) at $E=-0.5/(h\cdot N)$ A look into the level clustering for finite particle number $N$ visualizes the formation of the ESQPT, see Fig. 1. It is build by an infinite series of avoided level crossings of the neighboring excited states, as their gap closes at the energy of the ESQPT in the thermodynamic limit Solinas et al. (2008); Caprio et al. (2008); Stránský et al. (2015); Šindelka et al. (2017); Lerma-H and Dukelsky (2014b). Some of such neighboring states are marked by solid and dashed colored lines in the Fig. 1.

We visualize this gap closing by plotting the energy distance between the neighboring excited states $\left|i\right>$ and $\left|i+1\right>$ which are colored in Fig. 2 (top panel) for different values of $N$. Qualitatively, the form of the curves for the excited states and their finite-size scaling behavior is similar to the ground state curves for $\left|i=1\right>$. The cusp-like minimum of the curve $i$ is at a position $\gamma_{x}^{cr,i}/h$ which increases with $i$. It marks the region where level density around the state $\left|i\right>$ diverges in the thermodynamic limit and the ESQPT takes place. In the $E(\gamma_{x})$ representation, the flow of the states $\left|i\right>$ and $\left|i+1\right>$ will bring them to the energy of $E(\gamma_{x}^{cr,i})=-0.5/N$ at $\gamma_{x}^{cr,i}/h$, see Fig. 1 Stránský et al. (2014). For a further analysis of the gap closing in the excited states we will later-on define an adiabatic protocol approaching the point $\gamma_{x}^{cr,i}$ from the left. Therefore, we fit the left part of each curve $i$ in Fig. 2 (top panel) close to the cusp minimum by Sachdev (2007)

[TABLE]

for different $N$ values, here the $a_{i},c_{i}$ are the fitting parameters and $\gamma_{x}^{cr,i}$ is determined numerically by finding the energetic minimum of the difference $E_{i+1}-E_{i}$ by varying $\gamma_{x}$. The obtained numbers for the critical scaling $c_{i}$ of the states $\left|i+1\right>$ and $\left|i\right>$ are shown in Fig. 2(lower panel). The scaling behavior is quantitatively similar for different particle numbers $N$. The ground state scaling parameter $c_{i=1}$ approaches $1/2$ for $N\to\infty$, which is typical for a second order QPT Sachdev (2007). For the excited states $\left|i>1\right>$ $c_{i}$ decays till a constant value. Hence, the energy distance of higher neighboring states $i$ and $i+1$ develops close to the ESQPT region in a similar manner.

IV Adiabatic Dynamics

We now analyse the system behavior by continuously quenching it through the ESQPT, starting in an excited state of $\hat{H}_{+}$. To probe the influence of the ESQPT on the adiabatic system dynamics, we assume a time dependent interaction strength $\gamma_{x}(t)$, linearly increased from some starting value $\gamma^{(i)}$ at time $t_{\rm in}$ to the final value $\gamma^{(f)}$ at time $t_{f}$

[TABLE]

where the quench rate $Q$ determines the inverse speed of the protocol, and $t\in[t_{\rm in},t_{f}]$. Note that the ground state QPT for our model, Eq. (1), has been investigated with similar schemes Caneva et al. (2008); Campbell et al. (2015); Solinas et al. (2008).

Thus, starting in a time-local eigenstate $\left|i(t_{\rm in})\right>$ of the Hamiltonian $\hat{H}_{+}(t)$ Eq. (5) at time $t=t_{\rm in}$ with $\gamma=\gamma_{x}(t_{\rm in})$, we solve the time dependent Schrödinger equation

[TABLE]

Due to the time dependent Hamiltonian $\hat{H}_{+}(t)$ the wave package $\Psi(t^{\prime})$ at time $t^{\prime}>t_{\rm in}$ will not necessarily be the $i$th eigenstate $\left|i(t^{\prime})\right>$ of the Hamiltonian any more, it will strongly differ especially in the case of the non-adiabatic evolution. We solve this equation numerically by representing the operators in the ${\left|j,m\right>}$-Basis.

To analyze the wave package $\left|\psi(t)\right>$, we first introduce the so-called residual energy Caneva et al. (2008); Hwang et al. (2015)

[TABLE]

which is the difference between the energy expectation value at time $t$ and the adiabatic eigenenergy of the Hamiltonian $\hat{H}_{+}(t)$ at time $t$. Obviously, this is zero at $t_{\rm in}$ but changes if the adiabatic condition is not fulfilled. If the initial state is the ground state, the residual energy remains positive, but in case of excited states it can have both signs which can switch in time. The top panels of Fig. 3 show the residual energy for different initial eigenstates of the Hamiltonian $\hat{H}_{+}(t)$ at $t_{\rm in}$ and a fixed $Q$ value (left) or for different $Q$ values but a fixed initial eigenstate (right). Note, that we usually choose $t_{\rm in}=-10/h$, therefore $E_{i}^{res}$ can deviate from [math] at $t=0$, especially for smaller $Q$ values. For higher $Q$ values the deviation is negligible at the begin of the protocol, therefore we show the dynamics only for $\gamma_{x}(t)/h>0$ in Fig. 3. The $E_{i}^{res}$ shows a strongly oscillating behavior which crosses the zero value for all initial states except the ground state which is labeled by ’1’. The beginning of the oscillation phase indicates the breakdown of adiabaticity. This is due to the onset of the ESQPT which causes the gap closing between the neighboring excited states, compare with Fig. 2.

The wave package $\left|\psi(t)\right>$ can be represented in the time-local eigenbasis at each time $t$

[TABLE]

where $\left|\alpha_{i}(t)\right|^{2}$ gives the probability that the state $\left|i(t)\right>$ is occupied. Note, that the $\alpha_{i}(t)$ as the wave function $\left|\psi(t)\right>$ depend on the choice for the initial condition $\left|i(t_{\rm in})\right>$, see Eq. (9). For a better analysis of the adiabaticity loss we need a positive valued signature instead of the residual energy. The simplest choice is the fidelity Campbell et al. (2015); Santos et al. (2016)

[TABLE]

It measures the overlap of the state $\left|i(t)\right>$ with the wave package and deviates from unity in the non-adiabatic regime. The onset of deviation is shifted to higher $\gamma_{x}(t)$ values for higher initial eigenstates and coincides with the beginning of the oscillating phase in the residual energy, see Fig. 3 (middle left panel). The protocol quench scale $Q$, Eq. (8), has a strong effect on the adiabaticity loss, too (Fig. 3, middle right panel). For a slower protocol (larger $Q$) the $\gamma_{x}$ position for the adiabaticity loss shifts closer to the ESQPT value $\gamma_{x}^{cr,i}$, whereas for a faster protocol the adiabaticity break down shifts to a smaller $\gamma_{x}$ values, for an initial state $\left|i(t_{\rm in})\right>$.

However, in the ongoing discussion we use a quantity $\sigma$, which is not only positive but contains information about all eigenstates,

[TABLE]

which describes the energy width of the wave package $\left|\Psi(t)\right>$ in the time-local basis $\{\left|i(t)\right>\}$. 111Note: The standard deviation of a distribution ${\varepsilon_{j}}$ with probability $p_{j}$ is given by $\sqrt{\sum_{j}(\varepsilon_{j}-\bar{\varepsilon})^{2}p_{j}}$ and $\bar{\varepsilon}=\sum_{j}\varepsilon_{j}p_{j}$. Here, $\bar{\varepsilon}=\left<\psi\right|\hat{H}_{+}\left|\psi\right>$. Similar as in case of the residual energy, Eq. (10) and the fidelity, (12), the presence of the ESQPT in the excited states affects the dynamical evolution of the $\sigma(t)$, see Fig. 3 (lower panels). Its value increases if the adiabaticity is lost, which depends on the position of the ESQPT $\gamma_{x}^{cr,i}$ and the protocol quench scale $Q$, as we have discussed above. Summarizing, all quantities plotted in Fig. 3 consistently display the adiabaticity change.

For a deeper analysis we now investigate the scaling of the energy spread $\sigma(t)$ by plotting it for different $Q$ values, keeping $\gamma_{x}(t)$ fixed, see Fig. 4 (top panel). For each initial state $\left|i(t_{\rm in})\right>$ we choose $\gamma_{x}(t)$ to be close to their specific transition at $\gamma_{x}^{cr,i}$, see Fig. 2. Qualitatively, the scaling of $\sigma(Q)$ is similar for both excited and ground states, as there are two regions for $Qh^{2}<1$ and $Qh^{2}\gg 1$ where the scaling is linear in the double logarithmic representation, but with different slopes. A similar separation has been found recently for the Jaynes-Cummings model Hwang et al. (2015). Quantitatively, the scaling exponent $\alpha$ for $Qh^{2}\gg 1$ , which is the slope in the double-logarithmic representation, i.e.

[TABLE]

changes for different initial states, see Fig. 4 (lower panel, orange dots). For the ground state $i=1$ the exponent is around -0.35 and agrees with Ref. Hwang et al. (2015) and Silvi et al. (2016) in the thermodynamic limit. Note, that the Ref. Silvi et al. (2016) predicts the ’-1/3’ scaling in the classical regime, which applies for the LMG model with the well known mean-field, see Sec. VI. For higher excited states the slope decreases and tends to a constant value around -1. Interestingly, the slope variation can be achieved even for the ground state in quantum critical systems. Thus by varying the effective system size it can be decreased till the value of -2, which is typical of an adiabatic evolution in presence of a second order QPT Hwang et al. (2015); Caneva et al. (2008). In our case the system size $N$ is fixed and the change is caused by the excited states dynamics which suggests a correspondence between the dynamics around the ESQPT for higher states and the finite size dependent dynamics around the ground state QPT. Note, the form of the $\sigma(Q)$ curve echoes the gap scaling of the excited states, see Fig. 2. The constant exponent for $i/N>0.05$ again shows the dynamical indistinguishably of the excited states close to the ESQPT.

V State dynamics

For a further analysis of our protocol, Eq. (8), we plot the system dynamics exemplary for two different $Q$ values in Fig. 5 by plotting the occupation probability $P_{j}(t)=\left|\alpha_{j}(t)\right|^{2}$ for each time-local eigenstate $\left|j(t)\right>$. For a small $Qh^{2}=1$ value (left), the adiabaticity is lost shortly after $\gamma_{x}(t)/h=0$ and the distribution becomes spread between all eigenstates of $H_{+}(t)$, whereas for a higher $Qh^{2}=10$ value (right) the adiabaticity is lost at a much later stage $\gamma_{x}(t)/h\approx 1.25$, and the wave function is located mostly between the 15.th and 25.th eigenstate. We also recognize, that the width increases continuously for some time after the position of the ESQPT (green dotted line) is crossed Solinas et al. (2008). This representation suggests, that the ESQPT signatures are visible for higher $Q$ values and become smoothed for the smaller ones, as the wave package becomes then strong delocalized over many eigenstates.

Next, we investigate the dynamics in terms of the spin-coherent states which can be visualized on the Bloch sphere. An instantaneous spin coherent state $\left|\theta,\varphi\right>_{t}\equiv\left|\theta(t),\varphi(t)\right>$ rotates the Dicke state $\left|m=N/2\right>$ by the angles $\theta$ and $\phi$ and is defined as Dowling et al. (1994); Zhang et al. (1990)

[TABLE]

and can be rewritten using all eigenstates of the $\hat{J}_{z}$ as  Engelhardt et al. (2015)

[TABLE]

Thus, the wave function can be expanded in this new basis according to

[TABLE]

where

[TABLE]

The quantity $\left|\beta(\theta,\varphi)\right|^{2}$ gives the contribution of the corresponding coherent state to the wave function $\left|\psi(t)\right>$. Figure 6 shows snapshots of this probability for discrete $t\sim\gamma_{x}(t)$ values on the upper half of the Bloch sphere which is projected to its great circle. The gray lines show the position of coherent states with the same energy. Each of those lines, if the eigenstates were presented directly on the Bloch sphere, would mark the position of the center of a time-local eigenstate with the corresponding energy of the Hamiltonian $\hat{H}$ at time $t$. Compare therefore yellow dots and colored halo in the top left panel of Fig. 6.

The QPT of the ground state is visible in the structure change of the iso-energetic lines, thus for $\gamma_{x}(t)/h>1$ two different energetic structures are present which are separated by a separatrix shaped as the symbol ’$\infty$’. At the separatrix the ESQPT takes place. The separatrix grows with increasing $\gamma_{x}(t)$. The upper and lower figure panels show the evolution of the wave function in terms of spin coherent states for two different $Q$ parameters (the same as in Fig. 5) , one smaller (top panel, fast quench) and larger (lower panel, slow quench). The initial state is in both cases the same as in Fig. 5, thus the 17th excited eigenstate of the Hamiltonian $\hat{H}_{+}(\gamma_{x}(t=0))$ and coincides with one gray line in Fig. 6. Note that the width of the $\beta(\theta,\varphi)$ distribution becomes smaller with larger $N$. The adiabaticity loss due to the appearance of the QPT and the ESQPT at some later point in the dynamical evolution is visible by a deviation of the coherent state distribution from the corresponding gray line. For a faster protocol, the difference is clearly visible at $\gamma_{x}(t)/h=1$ (second figure, top panel), whereas for a slower protocol the wave package approximately follows the corresponding gray energy line even at $\gamma_{x}(t)/h=1.5$ (second figure, lower part). However, the protocol speed is essential for a further evolution. For a fast evolution (top panel) the wave package does not have enough time to adapt and keeps its shape till $\gamma_{x}(t)/h\approx 2$ (first four snapshots). Meanwhile the separatrix goes through the wave package and embeds it for higher $\gamma_{x}(t)$ values. For a slower change of $\gamma_{x}(t)$ (lower part) the wave package approximately still follows the corresponding energy line which collects it around the north pole on the Bloch sphere along the separatrix line (see third and fourth snapshot). Due to this fact the ESQPT signatures should be present for example in observables like $\left<J_{z}\right>(t)$ or $\left<J_{x}^{2}\right>(t)$. For a higher $\gamma_{x}(t)$ value the wave package is completely embedded by the separatrix.

VI Mean-field connection

The coherent state representation is closely connected to the mean-field treatment, which in case of the closed LMG model coincides with the quantum results in the thermodynamic limit Ribeiro et al. (2008). In this section we apply the protocol Eq. (8) at the mean-field level and compare the resulting semi-classical with the numerical quantum results from the previous sections. However, the mean-field provides the system information based on the operator expectation values and it is not possible to fix some energy level just by counting from the ground state, which we have explicitly used in the previous section.

To derive the mean-field equations we use the Heisenberg equation $\partial_{t}\left<\hat{J}_{\eta}\right>=i\left<\left[\hat{H},\hat{J}_{\eta}\right]\right>$ and assume that the factorization $\left<\hat{J}\eta\hat{J}_{\eta}^{\prime}\right>=\left<\hat{J}_{\eta}\right>\left<\hat{J}_{\eta}^{\prime}\right>$ for $\eta\in{x,y,z}$ holds in the thermodynamic limit, this yields Morrison and Parkins (2008); Dusuel and Vidal (2005)

[TABLE]

where $J_{i}(t)\equiv\frac{1}{N}\left<\hat{J}_{i}\right>$ denotes the rescaled average operator value.

To compare the results between the fully quantum and the mean-field calculations, we have to start with a similar initial condition in both cases. However, the initial condition in case of the mean-field calculation is a point on the Bloch sphere and corresponds to a spin-coherent state. In contrast, in the quantum case we have chosen an eigenstate of the Hamiltonian $\hat{H}_{+}$ at $t=t_{\rm in}$ as the initial condition. But an eigenstate is according to Eq. 17 a superposition of spin coherent states, see Fig. 6 for the visualization. Therefore, in the mean-field case we start with a set of initial states, whose energies correspond to the eigenenergy of the initial eigenstate in the quantum case we want to compare with. The dynamical evolution of coherent states is shown in Fig. 6 by yellow dots. We see that the quantum dynamics of the wave functions is captured by the mean-field calculation very well.

In case of the adiabatic evolution the energy of each coherent state in the set would evolve in the same way, thus the width of the energy distribution for a set of coherent states would be zero. However, in case of the adiabaticity breakdown the set of coherent states will evolve in a different manner and as a consequence the width of the energy distribution becomes finite. We will characterize the width by the standard deviation $\sigma_{M}$ and calculate it as a function of $Q$ for sets of coherent states with different initial energy $E_{\rm in}\equiv E(t=t_{\rm in})$ at a special time $\gamma_{x}^{cr,E_{\rm in}}(t)$. The latter is equal to $\gamma_{x}^{cr,i}$ if $\left<i(t_{\rm in})\right|\hat{H}_{+}\left|i(t_{\rm in})\right>=E_{\rm in}$ and gives the value when in case of the fully adiabatic evolution the initial set will cross the ESQPT region. Same as in the quantum case, the curve $\sigma_{M}(Q)$ has a linear part for $Qh^{2}\gg 1$ in the double-logarithmic representation (not shown here), which we fit by a power law

[TABLE]

The corresponding exponents $\alpha$ for different initial energies are shown in Fig. 4 by blue squares and matches pretty well to the purely quantum calculations. However the perfect agreement can not be expected due to the finite-size effects on the one hand and the factorization assumption at the lowest level on the other.

The ongoing protocol deforms the energy landscape (see Fig. 6), induces a symmetry breaking and creates a separatrix for $\gamma_{x}(t)/h>1$. Especially for a slower protocol (large $Q$) the set of initial states has a small ’energy width’ and goes around the same time through the growing separatrix. However, this corresponds to a system which goes through the ESQPT. This is visible in the observables behavior as a typical peak, in case of $J_{z}$ the peak arises around the value of ’0.5’. In Fig. 7 we show the $\bar{J}_{z}$ expectation value averaged over the set of different initial conditions with the same energy $E_{\rm in}$, thus over the $J_{z}$ values of yellow dots in Fig. 6, for different protocol speeds $Q$ as a function of $\gamma_{x}(t)$. Such averaging helps to suppress protocol-induced oscillations. One can see clearly a sharp increase (yellow dotted line) of $\bar{J}_{z}$ which occurs when the wave package comes close to the separatrix. For lower $Q$ the adiabaticity breaks down and the peak is moved to higher $\gamma_{x}$ values becoming more flat. For higher $Q$ values the position of the peak tends to $\gamma_{x}(t)/h=1.4$, which corresponds to the ESQPT position for the given initial energy $E_{\rm in}$. We could fit the $Q$ dependent position of the peak by $Q_{peak}=1/(\gamma_{x}(t)-1.4)^{2}$, where the exponent ’-2’ was the fitting parameter. Note that we have also checked that the quantum calculation leads to the quantitatively same results.

VII Summary

In this work we have investigated the role of the ESQPT for the adiabaticity loss of the excited states dynamics in the LMG model. Similar to the adiabaticity loss of the ground state evolution close to the QPT, the energy spectrum around the ESQPT in case of the LMG model is composed of a series of avoided crossings, which will make the dynamics of excited states non-adiabatic. To probe the adiabaticity loss, we have linearly increased the internal coupling in the model and crossed the ESQPT boundary by starting in different excited states. We argue, that the usually used residual energy to quantify the adiabaticity loss is not useful for the excited states, therefore we investigate the energy width of the wave package. Our results show that even for the excited states the energy width of the wave package close to the ESQPT scales with the power law as a function of the inverse protocol speed $Q$. Though, the corresponding exponent depends on the excited state. Whereas for the lower part of the spectrum there is a significant change in the exponents, for the higher part of the energy spectrum this exponent tends to a constant value.

Additionally, we showed that the mean-field and quantum based results agree. Therefore, we linked the eigenstates of the LMG Hamiltonian with a set of coherent states, which are the closest ones to the mean-field states. Thus, we investigated the evolution of the mean-field states set with initially same energy. Due to their different position on the Bloch sphere, each of those states evolves in a different manner under a protocol action. The deviation of the fully adiabatic evolution of this mean-field dynamics gives on average the same scaling behavior as a function of the protocol speed. In the Bloch sphere representation the protocol action can be understood as passing of the protocol-caused separatrix growth cutting through the wave package. Especially for a slower protocol evolution the ESQPT signatures can be identified in the dynamical evolution of the wave package.

Acknowledgments

We thank Georg Engelhardt for useful discussions. The authors gratefully acknowledge financial support from the DFG Grants BR 1528/9-1, SFB 910, GRK 1558.