Excited-state Quantum Phase Transitions in the Anharmonic Lipkin-Meshkov-Glick Model: Dynamical Aspects

PHYSICAL REVIEW E 107, 064134 (2023) Excited-state quantum phase transitions in the anharmonic Lipkin-Meshkov-Glick model: Dynamical aspects J. Khalouf-Rivera ,1,2,*J. Gamito ,3F. Pérez-Bernal ,2,4J. M. Arias ,3,4and P. Pérez-Fernández 1,4 1Departamento de Física Aplicada III, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, 41092 Sevilla, Spain 2Departamento de Ciencias Integradas y Centro de Estudios Avanzados en Física, Matemáticas y Computación, Universidad de Huelva, Huelva 21071, Spain 3Departamento de Física Atómica, Molecular y Nuclear, Facultad de Física, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain 4Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, Fuentenueva s/n, 18071 Granada, Spain (Received 22 February 2023; accepted 7 June 2023; published 26 June 2023) The standard Lipkin-Meshkov-Glick (LMG) model undergoes a second-order ground-state quantum phase transition (QPT) and an excited-state quantum phase transition (ESQPT). The inclusion of an anharmonic term in the LMG Hamiltonian gives rise to a second ESQPT that alters the static properties of the model [Gamito et al., Phys.Rev.E106, 044125 (2022)]. In the present work, the dynamical implications associated to this new ESQPT are analyzed. For that purpose, a quantum quench protocol is defined on the system Hamiltonian that takes an initial state, usually the ground state, into a complex excited state that evolves on time. The impact of the new ESQPT on the time evolution of the survival probability and the local density of states after the quantum quench, as well as on the Loschmidt echoes and the microcanonical out-of-time-order correlator (OTOC) are discussed. The anharmonity-induced ESQPT, despite having a different physical origin, has dynamical consequences similar to those observed in the ESQPT already present in the standard LMG model. DOI: 10.1103/PhysRevE.107.064134 I. INTRODUCTION The use of toy models has been fundamental for important advances in all branches of physics. These are nontrivial models but still simple enough to be solved analytically and they can be used either to look into limiting situations in complex systems or to check and better understand different approximation techniques. Some relevant examples of solvable models are Elliott’s rotational su(3) model [1] and the interacting boson model [2–5] in nuclear physics, the Rabi [6,7], Jaynes-Cummings [8], and Dicke models [9] in quantum optics, or the Lipkin-Meshkov-Glick (LMG) model in many-body physics [10–12], just to mention a few of them. In many cases, such models were originally introduced in a particular branch of physics and they were later used in completely different fields. In particular, the LMG model was originally proposed to test many-body approximations such as the time-dependent Hartree-Fock or perturbation methods in nuclear systems [10–12], but it has demonstrated to be very useful for the study of quantum phase transitions (QPTs) [13–16] and has been realized experimentally with optical cavities [17], Bose-Einstein condensates [18], nuclear magnetic resonance systems [19], trapped atoms [20–23], and cold atoms [24]. For instance, the LMG model has been used to test the possible existence of excited-state quantum phase transitions (ESQPTs) [25] and relations between ESQPTs and quantum entanglement [26,27], or quantum decoherence [28]. The ESQPT concept was introduced in [29] and an excellent review on this topic has been recently published [30]. *Corresponding author: [email protected] It is worth noting that phase transitions are well defined for macroscopic systems, however, the same ideas can be applied in mesoscopic systems where one can observe phase transition precursors even for moderate system sizes [31]. When dealing with mesoscopic systems, the study of their mean-field or large-size limit is a valuable reference to connect the precursors with the nonanaliticities expected in a QPT. Toy models, such as the LMG model, are simple enough to be solved for a large number of particles, allowing for a clear connection with the aforementioned large-size limit. This work is part of a more complete study on the anharmonic LMG (ALMG) model. The additional anharmonic term induces, in addition to the already known ESQPT [28,32], an anharmonicity-induced ESQPT that needs to be well understood. In a previous publication [33], the static aspects of both the ground-state QPT and the two ESQPT’s in the ALMG model were characterized. A mean-field analysis in the large-Nlimit was performed and different observables were used to characterize the different quantum phase transitions involved: The energy gap between adjacent levels, the ground-state QPT order parameter, the participation ratio, the quantum fidelity susceptibility, and the level density. In this work, we concentrate on the influence of the two ESQPTs on the dynamics of the ALMG model. With this aim, a quantum quench protocol that consists of an abrupt change in one of the control parameters in the ALMG Hamiltonian is defined. Then, the local density of states (LDOS, also known as strength function) together with the evolution of the system after the quench are studied using the time evolution of the survival probability, Loschmidt echoes, and an out-of-timeorder correlator (OTOC). 2470-0045/2023/107(6)/064134(13) 064134-1 ©2023 American Physical Society J. KHALOUF-RIVERA et al. PHYSICAL REVIEW E 107, 064134 (2023) The present paper is organized as follows. In Sec. II,the ALMG model is introduced, its algebraic structure reviewed, and the relevant matrix elements for the calculations in the u(1) basis are explicitly given. Section III is devoted to the analysis of a quantum quench protocol. Particularly, the time evolution of the survival probability when the system undergoes a quantum quench is discussed to understand how this quantity is influenced by the presence of the ESQPTs in the system. In Sec. IV, the ESQPTs impact on the evolution of an OTOC is explored. Finally, some conclusions are presented in Sec. V. II. THE MODEL The LMG model can be used to describe one-dimensional spin-1/2 lattices with infinite-range interactions [10–12]. For an array of Nsites, the Hamiltonian is written in terms of collective spin operators ˆ Sβ=N i=1ˆsi,β with β=x,y,zand where ˆsi,β is the βcomponent of the spin operator for a particle in site i. Therefore, the usual LMG Hamiltonian is written as ˆ H=(1 −ξ)(S+ˆ Sz)+2ξ SS2−ˆ S2 x,(1) with S=N/2. The operator ˆ Sxcan be written in terms of the usual ladder operators ˆ S+and ˆ S−, defined as ˆ S±=ˆ Sx±ıˆ Sy, and ξ∈[0,1] is a control parameter that drives the system from one phase to the other one. Indeed, from an algebraic point of view, the Eq. (1)LMG Hamiltonian presents a u(2) algebraic structure with two limiting dynamical symmetries: u(2) ⊃u(1) and u(2) ⊃ so(2) [34]. Each dynamical symmetry is associated with a different phase of the physical system. For ξ=0the system reduces to the u(1) dynamical symmetry and this phase is usually referred to as the normal (or symmetric) phase, whereas for ξ=1theso(2) dynamical symmetry is realized and the corresponding phase is called the deformed (or broken-symmetry) phase [34]. Inspired by the works in Refs. [35–37], we have included in the Eq. (1) Hamiltonian a second-order Casimir operator of u(2), S2 z, ˆ H=(1 −ξ)(S+ˆ Sz)+2ξ SS2−ˆ S2 x +α 2S(S+ˆ Sz)(S+ˆ Sz+1).(2) Again, the Hamiltonian depends on the ξcontrol parameter which drives the system between phases. In addition, a new control parameter, α, is introduced. The purpose of this work is to explore the influence of this new term and the corresponding control parameter on the dynamics of the system. It is worth noticing that for α=0, the original Hamiltonian, Eq. (1), is recovered, and for αdifferent from zero, the ξ=0 limit is transformed from a truncated one-dimensional harmonic oscillator to an anharmonic oscillator. That is the reason why Hamiltonian (2) is referred to as the anharmonic LMG model. Moreover, we observe that the so(2) limit is not longer recovered for ξ=1 unless αis zero. The Hilbert space for this system has dimension 2N,but due to the conservation of the total spin, [ ˆ S2,ˆ H]=0, we can focus on the sector of maximum irrep of the system, so the total spin quantum number S=N/2 through the work. This leads to a drastic reduction of Hilbert space dimension that now becomes N+1. However, the basis for the Hilbert space given by the subalgebra u(1), |S,Mzwith Mz= −N/2,...,0,...,N/2 (the projection of the total spin Son the zdirection), is used along this work. The matrix elements of Hamiltonian (2)intheu(1) basis are given by S,M z|ˆ Sz|S,Mz=MzδM z,Mz, S,M z|ˆ S2 z|S,Mz=M2 zδM z,Mz, ˆ S2 x=1 4(ˆ S2 ++ˆ S2 −+ˆ S+ˆ S−+ˆ S−ˆ S+), S,M z|ˆ S+ˆ S−+ˆ S−ˆ S+|S,Mz=NN 2+1−2M2 zδM z,Mz, S,M z|ˆ S2 +|S,Mz=N 2N 2+1−Mz(Mz+1)N 2N 2+1−(Mz+1)(Mz+2)δM z,Mz+2, S,M z|ˆ S2 −|S,Mz=N 2N 2+1−Mz(Mz−1)N 2N 2+1−(Mz−1)(Mz−2)δM z,Mz−2.(3) In addition, Hamiltonian (2) conserves parity (−1)S+Mzand the operator matrix can be split into two blocks, the first one including even parity states and the second one with odd parity states, with dimensions N/2+1 and dimension N/2 for an even Nvalue. A complete mean-field analysis of the semiclassical limit for Hamiltonian (2) has been carried out using spin coherent states in Ref. [33], revealing for α<0 a second-order groundstate QPT as well as two critical lines corresponding to two ESQPTs and both marked by a high density of states. A recently published work by Nader and collaborators focuses on a general LMG Hamiltonian that can be easily connected with our ALMG realization [38]. One of these high density of states critical lines was already known for the LMG model 064134-2 EXCITED-STATE QUANTUM PHASE TRANSITIONS IN … PHYSICAL REVIEW E 107, 064134 (2023) 0 0.2 0.4 0.6 0.8 1 ξ -0.2 0 0.2 0.4 0.6 0.8 ε ξ1 ξ2 εgs(ξ1) εc1(ξ2) εt=dεgs(ξ)/dξ|ξ1 (ξ−ξ1)+εgs(ξ1) εc1=ξ εc2=0.4 FIG. 1. Illustration of the tangent method discussed in the text for α=−0.6. Energy spectrum of the system in the plane ε×ξwhere the two ESQPT critical lines are highlighted with red and yellow dashed lines. The dashed blue line is the tangent for the ground-state curve εgs(ξ) at the point ξ1. This shows schematically the graphical determination of the critical quench ξ1→ξ2for a given initial state. In general, the intersection of the tangent line with the critical lines provides the critical ξ2value for which the system reaches the ESQPT critical energy after the quench. The dashed blue line stands for the tangent for the highest excited-state curve. [28,32]. Here, we pay heed to the other one, that we call anharmonicity-induced ESQPT critical line [33]. Particularly, it is worth exploring whether this critical line is of a similar nature as the other one and to what extent it has an impact on the system dynamics. For this purpose, the dynamics of the system is studied by means of the survival probability once the system undergoes a quantum quench and an out-of-time-order correlator (OTOC). III. QUENCH DYNAMICS The evolution of the system described by Hamiltonian (2) after a quantum quench should be sensitive to the presence of ESQPTs [28,38–41]. We explore the ESQPT influence on the system dynamics with a quantum quench protocol, starting from an eigenstate of the Hamiltonian, typically the ground state, and following the system evolution once a control parameter in ˆ His abruptly modified. The quenching brings the system to an excited state that evolves with time. The analysis of the ensuing system dynamics is a valuable tool to detect and explore ESQPTs in physical systems [28,39]. Let us just note that, from a mathematical point of view, this quenching analysis can be put in relation to the survival probability or a particular realization of the Loschmidt echo. The Hamiltonian in Eq. (2) depends on two control parameters, ξand α. In general, for negative αvalues, there exist two different ESQPTs and each one of them has a critical energy line marked by a high level density [33,38]. Since we are interested in characterizing both ESQPTs, a fixed value of α<0is selected and the time evolution of the system is explored after an abrupt change in the control parameter ξ. The aim is to study how the system dynamics is modified by the existence of two critical lines. In the followed quantum quench protocol, the system is initially prepared in a certain normalized eigenstate |0of ˆ H1=ˆ H(ξ1). At time t=0 a quantum quench takes place, changing ξfrom ξ1to ξ2. Thus, the Hamiltonian for the system is now given by ˆ H2=ˆ H(ξ2) and the initial state, |0, is no longer an eigenstate of ˆ H2and, consequently, evolves with time in a non trivial way. The probability amplitude of finding the evolved state, |0(t), in the initial state, |0, can be evaluated easily. The expression for this probability amplitude, denoted as a(t), is a(t)=0|0(t).The survival probability, F(t), also called nondecay probability or fidelity, is given by the absolute square of a(t), F(t)=|a(t)|2=|0|0(t)|2=|0|e−ıˆ H2t|0|2.(4) Since our goal is to evince the effect on the system dynamics of the external quench when reaching one of the ESQPTs’ critical lines, the determination of suitable ξ2values is very important, since the quenched system has to reach the corresponding critical energies. This can be achieved using the method of the tangent, developed in Ref. [40]. In Fig. 1,a typical evolution of the energy levels, ε, of the Hamiltonian in Eq. (2) is plotted as a function of the control parameter ξ, for a value of α=−0.6. In this figure, there is a change in the ground state at around ξ=0.2 that corresponds to the ground-state QPT. In addition, two lines of high level density in the excitation spectra are immediately apparent (separatrices, see Ref. [33]). These lines mark the critical energy of 064134-3 J. KHALOUF-RIVERA et al. PHYSICAL REVIEW E 107, 064134 (2023) the ESQPTs and separate the phases in such transitions. In the case shown in this figure, the separatrices occur at the critical energies εc1=ξ(yellow dashed line) and εc2=0.4 (red dashed line). A detailed discussion on this structure, including their dependence of the control parameters in the mean-field limit, can be found in Ref. [33], where the static properties of the ALMG model are presented. From Fig. 1, it is clear that for analyzing the three phases one has to start from the deformed phase ξ>ξ c=0.2. Due to the structure of our Hamiltonian, changing ξfrom an initial value ξ1implies that the system is excited along a straight line tangent to the energy line at ξ1. Thus, if the initial state is the ground state |0=|gsfor a particular ξ1value (ξ1>0.2), one needs to find the value of the ξparameter, ξ2, for which the tangent of the initial energy level ε1(ξ)atξ1crosses the critical line εc(ξ)atξ2in the plane ε×ξ. This is illustrated in Fig. 1for the case in which the eigenstate |0=|gsis the ground state of H1=H(ξ1). It is worth noticing that, within the range of values defined for ξ, using this method it is not possible to cross both ESQPTs lines from a given initial state. Indeed, for those values of ξabove the value of the critical ξcfor the QPT, it is only possible to reach the first ESQPT, εc1=ξ,(it can be seen plotting the tangent to the ground-state line). It is worth mentioning that if one uses the same tangent method starting from the symmetric phase (ξ<ξ c=0.2), one can reach the second ESQPT critical line, εc2=0.4, (yellow line), but it would be impossible to explore properly its impact on the dynamics of the system since one is forced to move over the first ESQPT critical line (red dashed line). Consequently, the tangent method from the system ground state is suitable for the study of the first ESQPT (red dashed line), but not the second one (yellow dashed line). Let us first examine the εc1=ξcritical line (red dashed line), that can be reached using the tangent method from the ξ1 ground state. On the one hand, the energy of the corresponding initial ground state is εgs(ξ1) and the equation for the tangent line at ξ1for the curve described by the ground state of the system in the ε×ξplane reads εt=m(ξ−ξ1)+εgs(ξ1), where mis the slope of the tangent to the ground-state curve at ξ1. On the other hand, the line of the first ESQPT (red dashed line) is εc1=ξ. Therefore, both lines cross at ξ2=ξc1=mξ1−εgs(ξ1) m−1,(5) where εgs(ξ1)=gs|ˆ H1|gs/N(ground-state energy per particle at ξ1) and the slope mof the tangent line is obtained making use of the Hellman-Feynman theorem in Eq. (2). Indeed, m=gs|ˆ H|gs/N=dεgs(ξ)/dξ|ξ=ξ1, where ˆ H=2 S(S2−ˆ S2 x)2−(S+ˆ Sz). A similar analysis can be performed for the anharmonicityinduced critical line. However, as we noticed above, the tangent to any point along the ground-state line with ξ>ξ c never crosses the second critical line (dashed yellow line) for the range of values of ξconsidered in this model. Hence, to explore this separatrix one should start from a more appropriate ˆ H1eigenstate. In particular, we have selected the highest excited state (denoted as |∗). As with the ground state, the most excited state of our system is well-defined in the thermodynamic limit by a coherent state [42]. Then, our initial state is now |0=|∗of ˆ H1. Let us denote the slope of the tangent to the energy line of the highest state at ξ1as m2. Then, this tangent line will reach the anharmonicity-induced ESQPT line given by εc2=ε0which is a constant. In the α=−0.6, the value of ε0=0.4 was computed with a mean-field formalism [33]. Therefore, the value for the critical ξ,ξc2, reads ξc2=m2ξ1+ε0−ε∗(ξ1) m2 ,(6) where ε∗(ξ1)=∗|ˆ H1|∗/Nand m2=∗|ˆ H|∗/N= dε∗(ξ)/dξ|ξ=ξ1is the slope of the corresponding tangent line. Once a way of crossing both ESQPT lines is available, the dynamic evolution of the system and the effect of crossing an ESQPT line can be examined. This can be accomplished computing the survival probability F(t)Eq.(4). Results for F(t)as a function of time are shown in Fig. 2for N=300, α=0(left column) and −0.6 (center and right columns) and different initial states (either the ground state |0=|gsin the left and central columns or the most excited state |0=|∗in the right column) for selected ξvalues. The α=0 case in the leftmost panels is included for the sake of completeness and reference. The panels in this column depict the time evolution of the survival probability for decreasing values of ξ2, starting always from the ground state |gsfor ξ1=0.6. The calculated ξ2at the crossing with the ESQPT is ξc=0.3. In general, the survival probability has a regular oscillatory behavior except in the region close to the ESQPT critical energy, ξ2=0.3, where the system undergoes an ESQPT and the survival probability suddenly drops down to zero and starts to oscillate randomly with small amplitudes. Once the critical energy for the ESQPT is crossed, the survival probability starts to oscillate in a regular way again. This phenomenon was reported for the first time in Ref. [28]. In the central and rightmost columns the same observable is plotted including a nonzero anharmonic term (α=−0.6). In the panels of the second column of Fig. 2, the survival probability is depicted for decreasing values of ξ2and starting always from the ground state of a Hamiltonian with ξ1=0.6 and α=−0.6. Due to the negative αvalue, the system undergoes two ESQPTs, displayed in the spectrum by means of critical lines with a noteworthy accumulation of energy levels (see Fig. 1). One of the two critical lines (red dashed line) can be traced back to the ESQPT already present in the α=0 case [25]. However, the second one (yellow dashed line) is linked to the presence of the anharmonic term in the Hamiltonian [33]. The nature and physical interpretation of the anharmonicity-induced ESQPT is different from the already known ESQPT associated with the ground-state QPT. Hence, in principle, there is no a-priori reason for both behaving in the same way. However, as we see if we compare the results for the critical ξcvalues in the different columns, the results obtained for the α=0 and the anharmonic cases are similar. The survival probability is oscillatory and regular except once ξ2is close to ξc, the critical value for the first or second ESQPT. In all cases, when ξ2=ξc, the quenched system reaches the critical energy and the survival probability suddenly drops down to zero and oscillates randomly with a small amplitude (red curves). This is similar to what happens in the α=0 case. Once ξ2is smaller than ξc, a periodic oscillatory decaying 064134-4 EXCITED-STATE QUANTUM PHASE TRANSITIONS IN … PHYSICAL REVIEW E 107, 064134 (2023) 0.25 0.5 0.75 1 F(t) ξ2=0.5 α=0, |ψ0>=|gs> ξ2=0.5 α= −0.6, |ψ0>=|gs> ξ2=0.6 α=−0.6, |ψ0>=|ψ∗> 0.25 0.5 0.75 1 F(t) ξ2=0.4 ξ2=0.3 ξ2=0.5 0.25 0.5 0.75 1 F(t) ξc=0.300 ξc=0.241 ξc=0.255 0 10203040 t 0.25 0.5 0.75 1 F(t) ξ2=0.1 010203040 t ξ2=0.1 010203040 50 t ξ2=0.1 FIG. 2. Survival probability F(t) as a function of time (t) for a system size N=300. The leftmost column includes F(t) results for α=0 and the middle and rightmost columns for α=−0.6. The initial state for the leftmost and middle columns is the ξ1=0.6 Hamiltonian ground state, |0=|gs, and the initial state for the rightmost column is the ξ1=0.7 Hamiltonian highest excited state, |0=|∗, to reach the second critical line of the energy spectrum. Different quantum quenches are shown for different values of ξ2. There are some critical values for ξ2,ξc, for which the system is settled in the critical energy of an ESQPT (third row) and the survival probability drops down to zero (with small random fluctuations). behavior is observed in F(t). As explained above, the quench from the ground state never reaches the second ESQPT line. For that purpose, one has to start from a different initial state. Thus, to explore how F(t) is affected by the second ESQPT, the quantum quench is performed using as an initial state the highest excited state of the system, |∗, for a given value of ξ1. In this way, the second critical line (yellow dashed line) for the ESQPT is accessible after the quench. In the panels of the right column of Fig. 2, the survival probability for decreasing ξ2values is plotted for an initial state equal to the highest excited state of the Hamiltonian with ξ1= 0.7 and α=−0.6. For this parameter selection, the second critical line is reached at ξc=0.255. In this column, again, results are very similar to the ones obtained in the preceding cases. The fidelity F(t) oscillates regularly while ξ2>ξ c2, but when the ξ2parameter gets close to the critical value, ξc2=0.255, the survival probability drops down to zero and randomly oscillates with a small amplitude. Once the critical line is crossed, F(t) recovers an oscillatory decaying periodic behavior, but at a certain time, this periodic oscillatory behavior becomes distorted. The reason for this phenomenon is that the tangent line to the highest excited-state curve at ξ1in the plane ε×ξremains very close to the critical line εc2=0.4 after the quench for lower values of ξ2up to 0. One should note that when starting from the highest excited state the first ESQPT critical line is not accessible after the quench (see Fig. 1). If we denote the eigenstates of ˆ H(ξi,α) with i=1,2as |ψj(ξi)with j=0,1,...,N/2, then we can write the initial state |0in the basis of ˆ H2=ˆ H(ξ2,α) eigenfunctions as |0=jCj|ψj(ξ2)and then F(t)=0|e−ıˆ H2t|0 2 = j |Cj|2e−iEjt 2 =dEe−iEt ρ0(E) 2 ,(7) where Ejis the energy of the jth ˆ H2eigenstate and ρ0(E)=j|Cj|2δ(E−Ej), called the strength function or local density of states (LDOS) [43,44], is the energy distribution of |0weighted by the Cjcomponents. From Eq. (7) it is clear that the fidelity F(t) is the absolute value of the LDOS Fourier transform squared and this quantity can provide some clues on the F(t) time dependence for the quench at the critical values ξc, denoted in red in the third row of Fig. 2.InFig.3we plot the LDOS for the same cases included in Fig. 2, hence in the first column, we show the LDOS for the ground state of ˆ H1=ˆ H(ξ1=0.6,α = 0) for different ˆ H2cases, all of them with α=0. In the 064134-5 J. KHALOUF-RIVERA et al. PHYSICAL REVIEW E 107, 064134 (2023) 0.00 0.25 0.50 0.0 0.1 0.2 0.3 2=0.6 0.0 0.2 0.4 0.00 0.05 0.10 0.15 2=0.5 0.0 0.2 0.00 0.01 0.02 0.03 2= 0.255 0.0 0.2 0.00 0.02 0.04 2=0.1 0.0 0.2 0.4 0.00 0.25 0.50 0.75 2=0.5 0.0 0.2 0.0 0.1 0.2 2=0.3 0.0 0.1 0.2 0.000 0.025 0.050 0.075 2= 0.241 0.0 0.2 Normalized energy 0.00 0.05 0.10 2=0.1 0.0 0.5 0.0 0.2 0.4 0.6 2=0.5 0.0 0.5 0.0 0.1 0.2 LDOS | j(2)| 0(1)| 2 2=0.4 0.0 0.5 0.00 0.02 0.04 0.06 2=0.3 0.0 0.5 0.00 0.05 0.10 2=0.1 FIG. 3. LDOS |ψj(ξ2)|ψ0|2as a function of the normalized excitation energy εj(arb. units) for systems with ξ1=0.6andα=0.0(left column), ξ1=0.6andα=−0.6 (middle column), and ξ1=0.7andα=−0.6 (right column) (N=300 in all cases). The chosen states are the ground state |ψ0=|gs(ξ1)(left and middle columns) and the most excited state with even parity |ψ0=|∗(ξ1)(right column), expressed in all cases in the basis of eigenstates of the Hamiltonian ˆ H(ξ2,α), being ξ2the quench parameter. The cases that correspond to a critical value of ξ2(third row) are plotted using red color. second and third columns we depict the LDOS for initial states that are the ground state of ˆ H1=ˆ H(ξ1=0.6,α =−0.6) and the most excited state of ˆ H1=ˆ H(ξ1=0.7,α =−0.6). The LDOS for the critical quench values are depicted in red. It can be clearly seen that, for all columns, in the critical control parameter cases the LDOS is nonzero at the ESQPT critical energy and has a clear local minimum at this energy value. The Fourier transform of such LDOS produces the particular time dependence shown in the panels of the Fig. 2third row. Another quantity of interest, inspired on Loschmidt’s objections to Boltzmann H theorem, is the Loschmidt echo [45,46]. This quantity, considered a probe to the sensibility of a system dynamics under perturbations, is used to benchmark the reliability of quantum processes [47]. It was shown to be a valid QPT detector [48] and, more recently, it has been used to check the influence of the ESQPT on the dynamics of the LMG model [49]. Consider an initial wave function, |ψ, which evolves a time tunder a Hamiltonian ˆ H1,|ψ(t)= e−iˆ H1t|ψ. We can reverse the time evolution with another Hamiltonian ˆ H2,eiˆ H2te−iˆ H1t|ψ. The squared overlap of the resultant state with the initial state |ψis the Loschmidt echo (LE), denoted as M(t)[46,50], M(t)=|ψ|eiˆ H2te−iˆ H1t|ψ|2.(8) 064134-6 EXCITED-STATE QUANTUM PHASE TRANSITIONS IN … PHYSICAL REVIEW E 107, 064134 (2023) Another physical interpretation of this quantity is possible, since Eq. (8) is the distance between the same initial state once it is evolved for a time twith two different Hamiltonian operators. One of the properties of ground-state and excited-state QPTs is that, near the critical region, states are quite sensitive to perturbations. A way to quantify this effect is computing the LE for the eigenstates of the system ˆ H1=ˆ H(ξ,α) with a time-reversal under ˆ H2=ˆ H(ξ+δ, α), Mj(t)=|ψj(ξ,α)|eiˆ H(ξ+δ,α)t|ψj(ξ,α)|2,(9) where |ψj(ξ,α)is the jth eigenstate of ˆ H1and δis a small perturbation. The LE, as well as its long-time average value, detects the ESQPT in the LMG model without anharmonicity [49]. In Fig. 4we plot Mj(t) for several eigenstates of a system with ξ=0.3 and α=−0.6. The total number of bosons is N=300, the system has been perturbed with δ=0.01, and only states with even parity are considered. Results are shown for j=0,20,48,82,103, and 120. The two states with energies closest to ESQPTs critical energies (j=48 and 103) are plotted in red. As expected, the ground state j=0 perform small oscillations with a single frequency around a value close to one, with a maximum value equal to one. Other states far from the critical region, as j=20,82, and 120, have a more complex oscillation pattern, not harmonic, with a larger amplitude and without reaching unity in the considered time range. However, in the case of eigenstates close to the critical energy, j=48 and 103, M(t) is only one for t=0 and the oscillations of the LE are of a much more irregular nature, something similar to what happens for the fidelity F(t) in Fig. 2. As shown in Ref. [49], the time averaged value of the LE for the jth eienstate, Mj, is a convenient probe to detect an ESQPT. This quantity is defined as Mj=lim T→∞ 1 TT 0 dt M(t)= kcδ jk 4,(10) where cδ jk are the coefficients of the jth eigenfunction of the Hamiltonian operator ˆ H(ξ+δ, α), expressed in the basis of eigenstates of ˆ H(ξ,α): |ψj(ξ+δ, α)=kcδ jk |ψk(ξ,α). The LE time averaged value is equal to the inverse of the participation ratio (PR) [51]of|ψj(ξ+δ, α)computed using the basis {|ψk(ξ,α)}.InFig.5we plot the time averaged LE versus the normalized excitation energy for all even parity states of the system studied in Fig. 4. The states included in Fig. 4have been marked using red pluses for critical ones (j=48 and 103) and blue crosses for others (j=20,82, and 120). Mjhas local maxima located for the eigenstates close to the critical energies, as it was observed in the LMG model without anharmonicity [49]. Hence, the time averaged LE detects the new ESQPT associated to the anharmonic term in the LMG Hamiltonian and confirms that this quantity is a good ESQPT probe. IV. ESQPTs AND OTOC Out-of-time-order correlators (OTOCs), that appeared for the first time in the 1960s in the context of superconductivity [52], are a four-point temporal correlation function able to measure the entanglement spread in a quantum system from FIG. 4. Loschmidt echoes for a system with ξ=0.3, α=−0.6, N=300, and a perturbation across the control parameter ξof δ=0.01. From top to bottom we display Mj(t) for the jth state with even parity: 0, 20, 48, 82, 103, and 120. The states closer to the ESQPTs are plotted in red. the degree of noncommutativity in time between operators. Since then, after a long period of relative inactivity, there has been a tremendous frenzy around this concept on various fronts [53]. They returned to the limelight with the proposal of OTOCs as a viable quantum chaos indicator, due to its exponential increase at early times in certain systems [54–57], and to diagnose the scrambling of quantum information [58–61]. Besides, OTOCs are sensitive probes for quantum phase transitions [62–68]. Despite the fact that the experimental access to out-of-time-order correlators is hindered by the unusual 064134-7 J. KHALOUF-RIVERA et al. PHYSICAL REVIEW E 107, 064134 (2023) FIG. 5. Time-averaged of M(t) versus the normalized excitation energy εfor the same system introduced in Fig. 4. The highlighted states (red plus symbols for transition states and blue crosses for others) correspond to the states studied in Fig. 4. time ordering of its constituents operators that precludes the measurement using local operators, several approaches using different experimental platforms have successfully provided OTOC results [23,69–75]. Given two operators, ˆ Wand ˆ V, it is possible to probe the spread of ˆ W(t) with ˆ Vthrough the expectation value of the square commutator Cw,v(t)=[ˆ W(t),ˆ V(0)]†[ˆ W(t),ˆ V(0)],(11) where ˆ W(t)=eıˆ Ht ˆ We −ıˆ Ht is the operator ˆ Win the Heisenberg’s representation [53,76–79]. The expectation value is usually computed in the canonical ensemble. However, in recent works, it has also been computed over given initial states or over the system eigenstates (microcanonical OTOC) [77,78]. The squared commutator Eq. (11) can be rewritten as Cw,v(t)=Aw,v(t)−2Fw,v(t). The first term is a two-point correlator, Aw,v(t)=ˆ W†(t)ˆ V†(0) ˆ V(0) ˆ W(t)+ ˆ V†(0) ˆ W†(t)ˆ W(t)ˆ V(0)and the out-of-time order appears in Fw,v(t), the real part of a four-point correlator, Fw,v(t)=[ˆ W†(t)ˆ V†(0) ˆ W(t)ˆ V(0)].(12) Without loss of generality, if we consider operators that are unitary, then Eq. (11) reads Cw,v(t)=2−2Fw,v(t). In a recent LMG model study, the ESQPT effects on the microcanonical OTOC and the OTOC following a quantum quench were explored for ˆ W=ˆ V=ˆ Sx/S[64]. The time evolution of the OTOC after a sudden quench was analyzed and it was concluded that the equilibrium value (the long time average value) of this observable can be used as a good marker for the ESQPT because it behaves as an order parameter, able to distinguish between the phases below and above the ESQPT, respectively. Our goal here is to analyze how the OTOC behaves once the ALMG system goes through the anharmonicity-induced ESQPT line. This study is of relevance since the physical nature of this ESQPT is different from the one of the already known ESQPT for the usual LMG model. Moreover, the possibility of using an OTOC as an order parameter for both ESQPTs is considered. We have used in our analysis the microcanonical OTOC [77,80], defined as Fn(t)=[n|ˆ W†(t)ˆ V†(0) ˆ W(t)ˆ V(0)|n],(13) where the state |nis the nth eigenstate of the Hamiltonian Eq. (2), whose energy is En. This state is computed for a given set of Hamiltonian parameters, ξand α. Following Ref. [64], we have first selected ˆ W=ˆ V=ˆ Sx/S as the OTOC operators. The reason behind this election is twofold. On the one hand, the expectation value of the ˆ Sx operator is known to be an order parameter for the QPT in the LMG model, and it has also been shown in previous works that it behaves as an order parameter for the ESQPT [81]. On the other hand, the ˆ Sxoperator is related with the breaking of parity symmetry in the spectrum eigenstates [43]. However, the obtained results (not shown) indicate that in this case the Fn(t) equilibrium value only detects the occurrence of the first ESQPT, independently of its nature, and not the second one. We decided to explore other possibilities such as ˆ W=ˆ Sy/S, ˆ V=ˆ Sx/Sor ˆ W=ˆ S+/S,ˆ V=ˆ S−/S. In both cases we obtain the expected results, with equilibrium values sensitive to the anharmonicity-induced ESQPT in the symmetric phase and to the two ESQPTs in the broken symmetry phase. Numerical solutions for the time evolution of the OTOC Eq. (13) with ˆ V=ˆ S−/Sand ˆ W=ˆ S+/Sare presented in Fig. 6. These are results for a selected set of positive parity states of a system with size N=300 that are obtained by the diagonalization of the Hamiltonian Eq. (2). The time evolution of the microcanonical OTOC is depicted for different initial states and ξ=0.5 with either α=0 (left-column panels) or α=−0.6 (right-column panels). Despite the different operators included in the OTOC, a quite similar phenomenology to that pointed out in Ref. [64] is observed. However, it is worth emphasizing that if we kept ˆ V=ˆ W=ˆ Sx/S, once the first critical energy is crossed, the time average value of the OTOC is zero as Fn(t) oscillates around zero. The behavior of the microcanonical OTOC, Fi(t), depends on the region of the spectrum in which the system is located. Particularly, Fi(t) develops a regular behavior, with small amplitude oscillations around a positive value. This value decreases until the Fi(t) oscillates around zero, when the critical energy value is reached. For the states close to the ESQPT critical energy (red color curves), not only Fi(t) oscillates around zero, but it also behaves in a highly irregular way, as in Ref. [64]. This is a feature shared by both columns in Fig. 6, though in the right column panels the second and fourth panel correspond to critical energies for the two ESQPTs that arise in this case. Let us now to discuss in more detail the left column (α=0). Recall that in this case there is just one ESQPT located in the mean-field limit at energy ε=ξ, its value for these plots is ε=ξ=0.5. We have selected the ground state and four other positive parity eigenstates, i=0, 30, 58, 100, and 140 in Figs. 6(a),6(c),6(e),6(g), and 6(i), respectively. The state with the closer energy to the critical ESQPT energy is i=58—Fig. 6(e)—where the ESQPT precursors are clearly manifested. In the cases with energies below the 064134-8 EXCITED-STATE QUANTUM PHASE TRANSITIONS IN … PHYSICAL REVIEW E 107, 064134 (2023) FIG. 6. Time evolution of the microcanonical OTOC, Fi(t), for selected positive parity eigenstates of an ALMG model with a system size N=300. In all panels ξ=0.5, the left column panels refers to Fi(t)forα=0 and the right column panels include results for α=−0.6. OTOCs for different initial states are shown: For the left column from top to bottom; (a) |i=0(ground state), (c) |i=30,(e)|i=58,(g) |i=100,and(e)|i=140. For the right column from top to bottom: (b) |i=0(ground state), (d) |i=74, (f) |i=95,(h)|i=115,(j) |i=140. There are some energies in which the eigenstate is settled at the critical energy of an ESQPT. These are the cases for panels (e) in the left column and (d) and (h) in the right column, highlighted using a red color. critical energy Fi(t)>0. However, as can be seen in Fig. 6(e), once the critical energy for the ESQPT is reached, the OTOC oscillates randomly around zero. For energies larger than the critical energy—left Figs. 6(g) and 6(i)—the Fn(t) display high and low frequency oscillations around a zero mean value. Therefore, the steady-state value of Fn(t) will be equal to zero for these states. As one goes up in energy in the spectrum, the same kind of oscillatory behavior is observed, with smaller amplitudes. There are some new features arising in the right column panels, that include Fi(t) results for the α=−0.6 anharmonic case. As previously mentioned, in this case there are two critical ESQPT lines that in the mean-field limit lie at ε=1+α=0.4 and ε=ξ=0.5. We show the results for the two eigenstates with local minimum PR values i=74 and 115 in Figs. 6(d) and 6(h).Again,theFi(t) OTOC oscillations at the critical lines are markedly irregular. These two states have been highlighted using red color. The other three values included in the right column of Fig. 6are i=0 (ground state), 95, and 140. In the region between the two critical lines, the envelope for Fi(t) has a sine-like oscillatory behavior around zero, so its steady-state value equals zero. Once the second critical line is crossed and the system energy increases, Fi(t) presents again an oscillatory behavior around positive values, as can be clearly seen in Fig. 6. It is worth pointing out that the characteristic times of the different microcanonical OTOCs span a wide range of frequencies. In particular, Fig. 6(f) exhibits a much longer period (smaller frequency) than the rest of the panels. The oscillatory frequency of the four-point correlator can be traced back to energy differences between pairs of states of different parity [68]. Therefore, whenever different parity eigenstates are degenerate, the stationary value of the OTOC has a nonzero contribution. This occurs at energies less than the critical energy of the first 064134-9