Quantum quench influenced by an excited-state phase transition

PHYSICAL REVIEW A 83, 033802 (2011) Quantum quench influenced by an excited-state phase transition P. P ´ erez-Fern´ andez,1P. Cejnar,2J. M. Arias,1J. Dukelsky,3J. E. Garc´ ıa-Ramos,4and A. Rela˜ no5 1Departamento de F´ ısica At´ omica, Molecular y Nuclear, Facultad de F´ ısica, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain 2Institute of Particle and Nuclear Physics, Faculty of Mathematics and Physics, Charles University, V Holeˇ sovi ˇ ck´ ach 2, Prague, 18000, Czech Republic 3Instituto de Estructura de la Materia, CSIC, Serrano 123, E-28006 Madrid, Spain 4Departamento de F´ ısica Aplicada, Universidad de Huelva, E-21071 Huelva, Spain 5Grupo de F´ ısica Nuclear, Departamento de F´ ısica At´ omica, Molecular y Nuclear, Universidad Complutense de Madrid, Avenida Complutense s/n, E-28040 Madrid, Spain (Received 20 September 2010; published 3 March 2011) We analyze excited-state quantum phase transitions (ESQPTs) in three schematic (integrable and nonintegrable) models describing a single-mode bosonic field coupled to a collection of atoms. It is shown that the presence of the ESQPT in these models affects the quantum relaxation processes following an abrupt quench in the control parameter. Clear-cut evidence of the ESQPT effects is presented in integrable models, while in a nonintegrable model the evidence is blurred due to chaotic behavior of the system in the region around the critical energy. DOI: 10.1103/PhysRevA.83.033802 PACS number(s): 42.50.Nn, 05.45.Mt, 64.70.Tg, 05.70.Fh I. INTRODUCTION Diverse quantum effects in systems depending on external control parameters represent an interesting field of theoretical and experimental investigation. A lot of recent attention in this field has been focused on two different types of dynamical phenomena, namely, on quantum phase transitions and socalled quantum quenches. These phenomena and their mutual relation are addressed in the present work. A quantum phase transition (QPT) is a sudden change of the ground-state structure at a certain critical value of the control parameter. It can be observed as a nonanalytic evolution of the system’s energy and wave function induced by an adiabatic variation of the control parameter λacross the quantum critical point λc0at zero temperature. First discussed in the 1970s [1–3], the QPT phenomena become very important in the context of solid-state physics [4–6] as well as in nuclear and many-body physics—see, e.g., reviews [7,8]. Recent experimental realization of a QPT by means of laser-driven Bose-Einstein condensates [9] makes the concept highly actual also in the framework of quantum-optical models, which are studied in the present work. A quantum quench (QQ) represents an abrupt, diabatic change λ1→λ2of the control parameter followed by unitary quantum evolution associated with the Hamiltonian at λ2.The dynamics after the quench is determined by the fragmentation of the initial state (typically an eigenstate of the Hamiltonian at λ1) in the Hamiltonian eigenstates at λ2and depends on the system-specific features of the energy spectrum and eigenfunctions in the relevant range of energy and control parameter. Pioneering theoretical works in this field appeared already in the late 1960s [10], but a really rapid growth of interest was triggered at the beginning of this millennium by experimental studies using Bose-Einstein condensates confined in optical lattices [11]. For an extensive list of QQ-related references, see Ref. [12]. Not surprisingly, the QPT and QQ effects are mutually related. The response of the system to a quantum quench is most complex if the initial state is prepared near a quantum critical point λc0[12–15]. Since the QPT is connected with a sharp variation of the system’s dynamical properties, the QQ-induced time evolution in these cases depends sensitively on the size of the parameter change. In particular, the quenches leading to the other side of the quantum critical point are expected to differ in their responses from the quenches ending on the same side. However, it is clear that a more detailed picture connecting the two effects is needed. This can be seen from the fact that the QQ-induced dynamics depends on the features of a number of excited states (those populated by the diabatic parameter change) while the QPT is by its definition related merely to the ground state. This problem is addressed in the present paper. In particular, we discuss the QQ-induced dynamics in connection with a novel concept related to quantum criticality—a so-called excited-state quantum phase transition (ESQPT) [16–18]. This phenomenon represents a nonanalytic evolution of individual excited states in the system with a variable control parameter λ. A chain of ESQPTs affecting energy levels Eiwith increasing excitation index iusually originates from the ground-state QPT affecting E0at λc0. So far, such effects have been studied mostly in integrable systems with one effective degree of freedom (one-dimensional configuration spaces) showing a singularity in their classical dynamics at a certain energy [19–26], but they seem to exist in a much richer variety of incarnations. The ESQPTs can be viewed as close analogs of thermal phase transitions formulated in the microcanonical language [27]. An important step is the scaling of energy and other observables by a suitably defined size parameter ℵ.Ifthe thermal fluctuation of the scaled energy E≡E/ℵvanishes in the thermodynamic limit ℵ→∞, a thermal phase transition becomes localized at sharp values of E. This leads to anomalous behavior of the level density ρ(E,λ) along the critical energy curve Ec(λ) and affects also the λ-dependent “flow” of individual energy levels Ei(λ). 033802-1 1050-2947/2011/83(3)/033802(14) ©2011 American Physical Society P. P ´ EREZ-FERN ´ ANDEZ et al. PHYSICAL REVIEW A 83, 033802 (2011) In standard thermodynamic systems, the size parameter ℵ is directly proportional to the number of degrees of freedom so that the ℵ→∞limit leads to an infinite dimension of the configuration space. However, for “finite models” including those investigated in the present work the dimension of the configuration space is fixed (the degrees of freedom in these models are usually some collective ones and their number does not depend on the number of constituents). In this situation, the singularities in ρ(E,λ) and Ei(λ) may be of a different type than those observed in standard thermodynamic systems. For instance, the level density can form an infinite peak at E=Ec(λ) (in contrast to thermodynamic systems whose level density is always an increasing function of energy). The peak in ρ(E,λ) is accompanied by vanishing slope and diverging curvature of individual energy levels crossing the critical boundary Ec(λ)[16–21]. One can anticipate that singularities in ρ(E,λ) and Ei(λ) may have a major impact on the dynamics following a quantum quench. Indeed, if the energy distribution of the initial state in the new Hamiltonian after the quench is spread across the critical energy Ec(λ2), the form of this distribution must be changed. In this case, the induced time evolution will be different from that in cases in which the distribution does not interfere with the critical energy. If the effect is sufficiently large, the QQ-induced dynamics would constitute a possible detection method for various types of ESQPT. To study the above-formulated conjecture, we use three simple models. They describe a single-mode bosonic field interacting with an algebraic subsystem, which is based either on the SU(1,1) or on the SU(2) dynamical algebra. The SU(1,1) model may serve as a toy for the description of formation and dissociation of diatomic molecules and bosonic atoms [28]. The SU(2) Hamiltonian represents either the well-known Dicke [29] or Jaynes-Cummings [30](TavisCummings [31]) model of quantum optics, or may alternatively describe an interacting mixture of diatomic molecules and fermionic atoms [28]. In this work we will adopt the former interpretation. While the SU(1,1) model and the SU(2)-based Jaynes-Cummings model are integrable, the SU(2)-based Dicke model is not. All the three models show a rather similar phase structure. First, considering the ground-state properties, it turns out that an increasing strength of interaction between the bosonic field and the algebraic subsystem drives the entire system to a quantum critical point where the ground state abruptly changes its form. We then show that this QPT is in all three models followed by a chain of ESQPTs and demonstrate that these have a strong impact on the character of quantum dynamics after some fine-tuned quenches—namely, those leading the system to a narrow region around the critical excitation energy. The type of ESQPT and its QQ signatures depend on the dimensionality of the model: they are of the strongest type for the SU(1,1) and the SU(2) Jaynes-Cummings Hamiltonians, and of a softer type for the SU(2) Dicke Hamiltonian (see below for a clarification of this concept). We will also see that breaking of integrability in the latter model blurs the effects of criticality in the quench dynamics. We anticipate the same general trend also in more complex situations. The plan of the paper is the following: In Sec. II we describe the models, and analyze in Sec. III their classical and phase-transitional properties, particularly those related to excited states. Section IV collects the results on quantum quenches. We introduce the general concept of the critical quench, driving the system into the ESQPT region, and continue to more specific numerical results for the models employed. Section Vcontains a brief summary and outlook. II. MODELS A. Algebraic structure Below we will investigate quantum quenches in three simple models formulated within a common algebraic framework. Before analyzing specific features of these models, we first outline the underlying algebraic structure. The models describe a composite system consisting of two interacting subsystems: (i) a single bosonic mode given by creation and annihilation operators b†and b, therefore described by the HeisenbergWeyl algebra HW(1), and (ii) a subsystem represented by pseudospin operators J±=Jx±iJyand J0=Jzsatisfying commutation relations of the SU(2) algebra, [J0,J±]=±J±,[J+,J−]=2J0,(1) or by analogous operators K±=Kx±iKyand K0=Kz satisfying commutation relations of the SU(1,1) algebra [K0,K±]=±K±,[K+,K−]=−2K0.(2) The SU(2) or SU(1,1) algebras will be realized more specifically in terms of fermionic or bosonic operators. The complete dynamical algebra is HW(1) ⊗SU(•), where the bullet stands for a specification of the respective special unitary algebra, but in the following we will use just the abbreviated names SU(2) and SU(1,1) for the two models. A schematic representation of both models is given in Fig. 1. The Hilbert space of the coupled system is identified with the tensor product H=H(i) ⊗H(ii), where H(i) is the space of HW(1) spanned by the set of basis vectors |Nb, with Nb=0,1,... denoting the number of bbosons, and H(ii) coincides with the space associated with one of the irreducible representations (irreps) of the group SU(2) or SU(1,1). The irreps are classified by the eigenvalues c(2) SU(•)of the respective second-order Casimir invariants, C(2) SU(2) =J2 x+J2 y+J2 z,(3) C(2) SU(1,1) =K2 x+K2 y−K2 z.(4) These are parametrized as c(2) SU(2) =j(j+1) (with jinteger or half integer) and c(2) SU(1,1) =−k(k−1) (with k>0 known as the Bergmann index). The irreps are finite dimensional in the SU(2) case (compact group) and infinite dimensional in the SU(1,1) case (noncompact group). The respective basis states |j,m(with m=−j, −j+1,...,+jbeing the eigenvalues of J0) and |k,n(with n=0,1,... enumerating the eigenvalues k+nof K0) are generated from the lowest state |j, −jand |k,0by consecutive actions of the raising operators J+and K+, respectively. 033802-2 QUANTUM QUENCH INFLUENCED BY AN EXCITED-STATE ... PHYSICAL REVIEW A 83, 033802 (2011) FIG. 1. (Color online) A schematic representation of the models used. The SU(1,1) model (the upper panel) may describe the coexistence of two-atom molecules (lower level) with dissociated atoms (upper level). The SU(2) model (the lower panel) describes the interaction of a single-mode radiation field with an array of two-level atoms. B. SU(1,1) model Since the group generated by the SU(1,1) algebra is noncompact, its irreps are infinite dimensional, the generators being expressible through creation and annihilation operators a†,a of another type of boson. For instance, a single-boson-pair realization reads as K+=1 2(a†)2,K −=1 2a2,K 0=1 2a†a+1 2.(5) Alternatively, one can use some other boson-pair realizations (e.g., with two kinds of boson), which together with Eq. (5) constitute various forms of the Schwinger representation of the SU(1,1) algebra. Note that in Sec. III A we will also introduce the Holstein-Primakoff bosonic representation. To construct the Hamiltonian, we assume the simplest realization (5). In this case, there are just two irreps, one with k=1 4and the other with k=3 4. Their respective Hilbert spaces are spanned by vectors |Nacontaining even and odd numbers of abosons. The interaction between aand bbosons is considered such that the creation of one bboson leads to the destruction of a pair of abosons, and vice versa. Such a model can schematically describe, for example, the formation and dissociation of two-atom molecules [28]. The total Hamiltonian reads as H(1) =ω0K0+ωb†b+λ √M(1) [bK++b†K−],(6) where λ/√M(1) ⩾0 is a scaled coupling parameter (the meaning of M(1) will be explained below) and ω,ω0stand for single-particle energies (we set ¯h=1). For each of the SU(1,1) irreps (classified by the quantum number k), there are two commuting operators (quantum degrees of freedom) which determine the basis in the whole Hilbert space of physical states: one is associated with the number of bbosons, Nb=b†b, the other with K0,or equivalently with the number of abosons, Na=a†a.Atthe same time, there exist two different integrals of motions: one is the energy Hand the other one can be written in the form M(1) =2Nb+Na−4k−1 2=2(Nb+K0−k).(7) The value of M(1) ⩾0 is always even, M(1)/2 counting the number of bbosons plus the number of a-boson pairs. The conservation of M(1) implies that the Hamiltonian (6) represents an integrable system, which for each fixed value of M(1) can be associated with an effective one-dimensional configuration space (one quantum degree of freedom). In the following we will assume ω0>ωfor the SU(1,1) model. This means that the λ=0 ground state can be identified with a molecular condensate with no pair of atoms. It has the form |Nb=M(1)/2⊗|k,0, where the first term represents a state with a maximal number of bbosons and the second one stands for the lowest-weight SU(1,1) state with the minimal value of Na=0(Na=1) for k=1 4(k=3 4). However, for sufficiently large values of the coupling parameter λthe interaction between the molecules and atomic pairs supports a more balanced distribution of the expectation values Na and 2Nb. With an increasing size of the system the crossover between the two types of ground-state structure gets sharper and in the infinite-size limit, M(1) →∞, it becomes a phase transition. The calculation of the critical value of the interaction strength, λ(1) c0, will be presented in Sec. III B. C. SU(2) models The SU(2) algebra yields a compact group with finitedimensional irreps. Its generators can therefore be constructed from fermionic operators. For instance, they can be associated with an array of spin-1 2particles (or two-level atoms) located on 2jsites: J+= 2j  i=1 a† ↑ia↓i,J −= 2j  i=1 a† ↓ia↑i, (8) J0=1 2 2j  i=1 (a† ↑ia↑i−a† ↓ia↓i). Here, a† ↑ior a↑iand a† ↓ior a↓icreate or annihilate spinup and spin-down states of the fermion on site iand the ladder operators J±describe spin flips along the array. Alternatively, one can use fermion-pair realizations of the SU(2) algebra (with J±creating and annihilating a pair of fermions), or the Schwinger or Holstein-Primakoff types of bosonic realization with truncated Hilbert spaces (the latter bosonic realization will be discussed in Sec. III A). Depending on the specific realization, the model can receive different physical interpretations. Below we will implicitly consider the realization (8), which may schematically describe interactions of single-frequency photons with two-level atoms in maserlike systems. 033802-3 P. P ´ EREZ-FERN ´ ANDEZ et al. PHYSICAL REVIEW A 83, 033802 (2011) The Hamiltonian is taken in either of the following forms: H(2) =ω0J0+ωb†b+λ √M(2) [bJ++b†J−],(9) H(3) =ω0J0+ωb†b+λ √M(3) [(b+b†)(J−+J+)],(10) where λ/√M(2) or λ/√M(3) is a properly scaled coupling parameter (λ⩾0) and ω,ω0two single-particle energies. The Hamiltonian H(2) is known as the Jaynes-Cummings [30]or Tavis-Cummings model [31], while the Hamiltonian H(3) is referred to as the Dicke model [29]. Note that a so-called rotating-wave approximation of H(3) leads to the simpler Hamiltonian H(2). The Jaynes-Cummings Hamiltonian (9)isverysimilarto that of Eq. (6). It conserves the quantity M(2) =2(Nb+J0+j),(11) analogous to Eq. (7), and therefore corresponds to an integrable system described effectively by a one-dimensional configuration space. (The full model has again two degrees of freedom, associated with commuting operators Nband J0.) The Dicke model violates the conservation of M(2), but it still conserves the parity =(−1)M(2)/2labeling individual eigenstates. In this case, the size parameter is taken as M(3) =4j, which is the total number of fermionic states (twice the number of sites). As in the SU(1,1) case, the ground states of both the SU(2) Hamiltonians change their nature suddenly as the coupling strength λincreases above a certain value, the transition having a critical character in the infinite-size limit, M(2),M(3) →∞. For the Hamiltonian H(2) from Eq. (9), we will assume ω>ω 0, identifying the λ=0 ground state with a photon vacuum, Nb=0, combined with a maximally excited state of the atom array: J0=1 2M(2) −j(below we set M(2) =4jso that J0= +jat λ=0). At the critical coupling strength, this structure eventually changes into a state with Nb>0, in which a part of the energy is transferred from atoms to the photon field. It should be stressed that here we are running the model in a nonstandard regime, taking into account only a finite set of states with a single fixed value of M(2). This is in contrast to the rotating-wave approximation of the Dicke model, for which one usually considers the infinite spectrum with all values of M(2). For the Dicke Hamiltonian H(3) in Eq. (10), we set ω0=ω, which corresponds to the resonance absorption and emission of photons by the atoms. The λ=0 ground state has the form |Nb=0⊗|j, −j, describing a photon vacuum and an unexcited array of atoms (recall that for the Dicke model Nb+J0is not conserved). For a sufficiently strong interaction between matter and light, the ground state flips to a form with J0>−jand Nb>0, showing a macroscopic excitation of both subsystems [32,33]. This may be considered as a toy example of the maser phase transition. D. Numerical solution The model Hamiltonians from Eqs. (6), (9), and (10) can be diagonalized numerically in an appropriate basis. Due to the tensor product structure of the Hilbert space, the basis is naturally chosen in the form |Nb⊗|k,nfor the SU(1,1) model, and |Nb⊗|j,mfor both SU(2)-based models, where |Nb∈H(i) stands for a state with a given number of bbosons while |k,n,|j,m∈H(ii) are basis vectors of the respective SU(1,1) or SU(2) irreps. Recall that in the SU(1,1) case, the link of |k,nwith the states |Na, counting the number of abosons, is achieved via setting k=1 4or 3 4for the evenor odd-Nairreps, respectively, and n=1 2(Na−4k−1 2). Matrix elements of individual Hamiltonian terms in these bases can be easily calculated from the known action of the b†and b operators on the vectors |Nband the action of {K+,K−,K0} or {J+,J−,J0}operators on vectors |j,mor |k,n. For both SU(1,1) and SU(2) integrable models, the basis includes a finite set of vectors. These are determined by the chosen values of the size parameters M(1) or M(2), which permit only a finite number of combinations of (Nb,n)or(Nb,m) satisfying Eqs. (7) and (11), respectively. Thus the corresponding Hamiltonian matrices are finite and the diagonalization is just a routine problem. The nonintegrable Dicke model, on the other hand, has no conservation-dictated constraint on the allowed combinations of basis vectors. Its basis is therefore infinite and must be numerically truncated for Nb>N trunc, making the convergence tests of the diagonalization outputs an important issue. Speaking qualitatively, the truncation with a fixed Ntrunc can be safely done in the low-energy part of the spectrum only for a sufficiently small interaction strength λbetween the bbosons (photons) and atoms. Indeed, we know that for λ=0, the Dicke Hamiltonian is diagonal in the |Nb⊗|j,mbasis, the states with increasing Nbbeing associated with increasingly high excitations. Therefore, for moderate values of λthe states with high photon numbers are only weakly admixed to the low-lying energy eigenstates of the system. As λincreases, the cutoff parameter Ntrunc must increase accordingly for a given low-energy portion of the spectrum to be well reproduced. We stress that the results presented below were tested for stability against the change of Ntrunc (see Sec. IV D). The setting of the model parameters, as used in the calculations below, and some important model-specific values are summarized in Table I(some of the symbols will be explained later). TABLE I. Summary of parameter setting for the three models employed: frequencies ω,ω0, the size parameter M(n), a parameter R(n)in the potential energy, the QPT critical parameter λ(n) c0,andthe ESQPT critical scaled energy E(n) c. SU(1,1) SU(2) n=1n=2n=3 Naeven Naodd Jaynes-Cummings Dicke ω,ω0ω0−ω=1=ωω−ω0=1=ω0ω=1=ω0 M(n)2Nb+Na2Nb+Na−14j4j R(n)1 2M(1) 3 2M(1) 1 2 1 2 λ(n) c0 1 √2 1 √2 1 √2 E(n) c 1 2 1 4−1 4 033802-4 QUANTUM QUENCH INFLUENCED BY AN EXCITED-STATE ... PHYSICAL REVIEW A 83, 033802 (2011) III. PHASE TRANSITIONS A. Classical limit The above Hamiltonians are specimens of a rather large general class of systems—namely, those described by finite algebraic models [34]. For these models, the relevant observables are constructed in terms of a finite set of generators Giclosing a dynamical algebra [Gi,Gj]=kcij k Gkwith structure constants cij k . The corresponding systems have a finite number of degrees of freedom and their thermodynamic (infinite-size) limit coincides with the classical limit ¯h→0. To see this, recall that the thermodynamic limit is generally achieved for asymptotic values of a properly defined size parameter ℵsuch that thermal fluctuations of a scaled Hamiltonian H=H/ℵvanish with ℵ→∞. In algebraic systems, this parameter needs to be introduced on the level of individual generators, via scaled generators Gi≡Gi/ℵκ(with κ>0) whose substitution into the Hamiltonian H(Gi) should be consistent with the definition of H; thus H(Gi/ℵκ)= H(Gi)/ℵ. In fact, this is why the size parameter ℵ≡M(n) was included into the effective coupling constant λ/√ℵ of the above Hamiltonians H(n)(n=1,2,3). The known commutation relations for the bare generators Githen ensure that the scaled generators Giyield vanishing commutators in the ℵ→∞limit, [Gi,Gj]→0, which constitutes the classical behavior of the correctly scaled observables. A general method for approaching the classical limit in finite algebraic models is based on coherent states [35,36]. These for the above-described composite systems are naturally considered in the form of a tensor product |ζ⊗|ξ, where |ζ∝eζb†|0(with ζ∈C) is the HW(1) coherent state of the subsystem (i), and |ξis a yet unspecified SU(1,1) or SU(2) coherent state of the subsystem (ii) [37,38]. The latter states can be taken in several alternative forms, depending on the concrete realization of the two algebras. One possibility is to use |ξ∝eξK+|k,0or |ξ∝eξJ+|j, −j(with ξ∈C) and associate the classical limit with j→∞or k→∞.The corresponding phase space of subsystem (ii) is then identified with a two-dimensional (2D) surface of constant positive or negative curvature, which is the sphere j2 x+j2 y+j2 z=1in the SU(2) case (ji=Ji/j ) or a two-sheet hyperboloid k2 x+ k2 y−k2 z=−1 in the SU(1,1) case (ki=Ki/k). We do not directly follow this path, partly because in our case the value of the SU(1,1) invariant is fixed to k=1 4or 3 4, so we would not be able to keep the same treatment in both SU(1,1) and SU(2) systems. Instead, we employ the HolsteinPrimakoff transformation of both algebras onto a bosonic field c†,c, which reads as K+=c†(2k+c†c)1/2,K −=(2k+c†c)1/2c, K0=c†c+k(12) in the SU(1,1) case (with c†c≡Nc⩾0) and J+=c†(2j−c†c)1/2,J −=(2j−c†c)1/2c, J0=c†c−j(13) in the SU(2) case (with 0 ⩽Nc⩽2j). In this way, we obtain a mapping of the original HW(1) ⊗SU(•) dynamical algebra onto a new algebra associated with both band ctypes of boson, which is then analyzed with the aid of bosonic coherent states |ζ,ξ∝eζb†+ξc†|0, where |0is a common vacuum of both b and cbosons. For Hamiltonians (6) and (9), the new dynamical algebra can be identified with the algebra U(2) ≡{b†b,c†c,b†c,c†b}, since in both these cases the total number of bosons N=Nb+ Ncis conserved. The size parameters M(1) and M(2) introduced above both coincide with the value 2N. On the other hand, for the nonintegrable Hamiltonian (10) the new dynamical algebra can be identified with HW(2), the Heisenberg-Weyl algebra of band cbosons. The total number of bosons is not conserved and therefore the only sensible size parameter is the value of j(for consistency reasons we have chosen M(3) =4j), which measures the size of the subsystem (ii). The method proceeds via evaluating the expectation value ζ,ξ|H(n)|ζ,ξ/M(n)for any of the above Hamiltonians H(n) (with n=1,2,3) in the coherent states of band cbosons. The result can be gained directly by substituting c √M(n)=x+ip √2,b √M(n)=y+iq √2(14) and the Hermitian conjugate expressions for c†and b†into the scaled Hamiltonian H(n)=H(n)/M(n). Note that in Eq. (14) we define coordinates x,y and the associated momenta p,q, respectively, which satisfy canonical commutation relations [x,p]=[y,q]=i/M(n)with 1/M(n)playing the role of the Planck constant. For M(n)→∞, the coordinate and momentum operators can be treated as commuting variables. We obtain a general form H(n)=H(n) 0+λH(n)(15) for the scaled classical Hamiltonian. The first term, H(n) 0=−R(n)ω0 2+ω0 2(p2+x2)+ω 2(q2+y2),(16) which describes the system with λ=0, is common to all three models n=1,2,3, with the additive constant expressed as R(1) =2k/M(1),R(2) =2j/M(2), and R(3) =2j/M(3) =0.5 (cf. Table I). The second term corresponds to the interaction and has a model-dependent form, H(1) =2R(1) +(p2+x2)(xy +pq)/√2, H(2) =2R(2) −(p2+x2)(xy +pq)/√2,(17) H(3) =2R(3) −(p2+x2)√2xy. The constraints (7) and (11) on the conservation of M(1) and M(2) read as p2+q2+x2+y2=1,(18) which makes it possible to completely eliminate one degree of freedom in both integrable models. To do this, we set one of the momenta to zero, in our case q=0, and use Eq. (18) to fix the corresponding coordinate: y=± 1−p2−x2.On the level of coherent states, the choice q=0 is achieved by considering only a relative phase between band cbosons, setting the overall phase factor to unity. Note that this choice is dynamically consistent since the elimination of yensures that ˙ q=∂H/∂y =0. 033802-5 P. P ´ EREZ-FERN ´ ANDEZ et al. PHYSICAL REVIEW A 83, 033802 (2011) B. Ground-state phase transitions To analyze the ℵ→∞properties of the ground state of the above three models, we set both momenta pand qin Eqs. (16) and (17) to zero, yielding a potential V(n)=H(n)|p=q=0 (n=1,2,3). Indeed, any increase of por qtakes one away from the minimum of the Hamiltonian H(n)(p,q,x,y) and therefore corresponds to an excitation of the system above the ground state. The energy obtained by the minimization of the potential V(n)(x,y) represents an estimate of the scaled ground-state energy E0=E0/ℵ. The problem can be further simplified for both integrable systems, i.e., for Hamiltonians H(1) and H(2), where the constraint (18) with p=q=0 restricts the ground-state solution to the unit circle x2+y2=1. In the SU(1,1) case (n=1) with ω0>ω,wetake x=sin ϑ, y =cos ϑ, (19) yielding the potential V(1) =V(1) 0+ω 2sin2ϑ+λ √23sin(2ϑ)2R(1) +sin2ϑ. (20) For the integrable SU(2) Hamiltonian (n=2) with ω>ω 0, we redefine the angle ϑso that x=cos ϑ, y =sin ϑ, (21) and V(2) =V(2) 0+ω 2sin2ϑ+λ √23sin(2ϑ)2R(2) −cos2ϑ. (22) In both cases ω ≡|ω0−ω|>0, while V(1) 0=(ω+ R(1)ω0)/2 and V(2) 0=(1 −R(2))ω0/2. For n=3, the nonintegrable SU(2) Hamiltonian H(3) of the Dicke model with ω0=ω, the constraint (18) is not applicable. We take polar coordinates x=rcos ϑ, y =rsin ϑ, (23) arriving at V(3) =V(3) 0+ω 2r2+λ √2r3sin(2ϑ)2R(3) −cos2ϑ, (24) with V(3) 0=−R(3)ω/2. It is clear that the M(1) →∞limit of the SU(1,1) model with k=1 4,3 4gives R(1) →0. On the other hand, in the integrable version of the SU(2) model it is natural to take M(2) =4j, that is, R(2) =0.5. Indeed, with this choice the total number of bosons Ncan be arbitrarily partitioned into Nband Nc, including the extremal choice (Nb,Nc)=(0,N), which corresponds to the photon vacuum combined with a fully excited array of atoms. With these settings, both expressions (20) and (22) become identical (except for the additive constants). The potential energy surface V(1) alias V(2) is shown in Fig. 2. The minimum of both potentials V(1) and V(2) for λ=0 is at ϑ=ϑ0(0) =0 (we may equivalently choose ϑ0=π, which would have no influence on the conclusions below). For FIG. 2. The potential energies V(1)(ϑ)andV(2)(ϑ) of the integrable SU(1,1) and SU(2) models as functions of a rescaled coupling parameter g=λ/(√2ω). The thick curve demarcates the trajectory of the potential minimum; the dashed line indicates the saddle point position. increasing λ, the minimum ϑ0(λ) remains at the same place until the critical value λ(1) c0=λ(2) c0=ω √2(25) is reached. Here, ϑ=0 becomes a saddle point and the minimum ϑ0(λ) deviates to negative values, following a trajectory sin2ϑ0=12g2−1−12g2+1 18g2,g≡λ √2ω ⩾1 2.(26) This nonanalytic evolution represents a second-order quantum phase transition. From Eq. (26)forg>g c0=0.5 we get ϑ0(g)∼√g−gc0, which means that the critical exponent for the order parameter ϑ0(or x0=sin ϑ0) is equal to 1 2. For the Dicke Hamiltonian H(3) the dimension of the system is not reduced, so the properties of potential (24)mustbe analyzed in the plane (r,ϑ)≡(x,y). Nevertheless, the physics is similar to that described above. For λgrowing from 0 to a critical value λ(3) c0=ωω0 2(27) (λ(3) c0=ω/√2forω0=ω), the potential minimum is located at (x,y)=(0,0), which corresponds to a separable state of unexcited atoms and the field vacuum. At the critical point (27), the determinant of the Hessian matrix (composed from second derivatives of V(3) with respect to both variables) evaluated at the minimum becomes negative, which means that (x,y)=(0,0) becomes a saddle point of the potential. Starting at this point, two degenerate minima deviate symmetrically to the quadrants with xy < 0 (for λ>0). These minima correspond to a (nearly) degenerate parity doublet of the 033802-6 QUANTUM QUENCH INFLUENCED BY AN EXCITED-STATE ... PHYSICAL REVIEW A 83, 033802 (2011) FIG. 3. The potential energy surface V(3)(x,y) of the Dicke model with ω0=ω. Values of a rescaled coupling strength g≡λ/(√2ω) are given in each panel. ground-state solutions involving excitations of both atomic and field subsystems. The distance r0of the minima from the origin increases with g≡λ/√2ωω0as √g−gc0above the critical point gc0=0.5, so we have again a second-order QPT with the order parameter r0characterized by the critical exponent 1 2. Various stages of evolution of the potential V(3) are shown in Fig. 3. C. Excited-state phase transitions Any excited-state phase transition can be recognized in the dependence of a quantum level density ρ(E,λ)onthe scaled energy E. At the ESQPT point Ec(λ) this dependence shows a nonanalyticity whose type enables one to classify the critical behavior in analogy with the standard typology of thermal phase transitions. The nonanalyticity of ρ(E,λ) is reflected by specific discontinuous features in the flow of scaled energy levels Ei(λ) through the critical boundary Ec(λ)[17,18]. The semiclassical theory of the level density [39] leads to the decomposition ρ(E,λ)=¯ρ(E,λ)+˜ρ(E,λ),(28) where ¯ρand ˜ρrepresent smooth and oscillatory components, respectively. The oscillatory component can be expressed as a sum over periodic orbits in the general form ˜ρ= kAkcos(¯h−1Sk+φk), where Akand φkstand for the amplitude and the phase shift of the kth-orbit contribution, while Sk=p·dxrepresents the action over this orbit. In the limit ¯h→0 this part leads to infinitely rapid oscillations which cancel out if the level density is integrated over an arbitrary narrow interval of energy. In this limit, only the smooth component of Eq. (28) is relevant. It is expressed via orbits of zero length, yielding the formula ¯ρ(E,λ)=(2π¯h)−fδ(E−H(p,x,λ))dpd x   (E,λ) ,(29) where His the classical Hamiltonian depending in general on f-dimensional vectors of coordinates xand momenta p.The quantity (E,λ)dErepresents a 2f-dimensional volume of the available phase space for the interval of energy (E,E+dE). As follows from these considerations, in systems with synonymous thermodynamic (ℵ→∞) and classical (¯h→0) limits, any kind of nonanalyticity in the energy dependence of the classical phase space volume generates an ESQPT on the quantum level. Such nonanalyticities most commonly follow from the presence of the Hamiltonian stationary points [40,41]. The relation between such points and thermodynamic phase transitions was recently investigated [27] for systems with asymptotically increasing numbers of degrees of freedom. The situation is similar also in the present type of system, described by finite algebraic models, in which the number of degrees of freedom is fixed (independent of the increasing size parameter). A moderate dimension of the phase space allows for nonanalyticities of different types than those typically observed in the infinite-dimensional systems. Numerous examples of such effects can be found in the literature; see, e.g., Refs. [16–20,42]. Because the quantum microcanonical entropy is proportional to a logarithm of ρ(E,λ), the dependence of the level density on the scaled energy can be used as a key for the ESQPT classification. For instance, a jump of (E,λ)atEc(λ) may be seen as an analog of a first-order phase transition. It is known that in systems with f=1 there exist ESQPT effects even stronger than those of the first-order type. These effects are associated with an infinite peak of (E,λ), which shows up the corresponding dependence ρ(E,λ). The origin of this behavior is often found in a local maximum of the one-dimensional potential at energy Ec(λ). On the other hand, a softer type of nonanalyticity in (E,λ), like a discontinuous or infinite derivative, causes a continuous phase transition. This typically happens in systems with more than one degree of freedom [18]. In case of a discontinuous (n−1)th derivative (n⩾2) the transition can be declared to be of nth order. The ESQPT effects are present in the spectra of all three models described above. We saw in Sec. III B that the ground-state QPTs are located at the critical points λ(n) c0from Eqs. (25) and (27). It turns out that for λ>λ (n) c0the singularity propagates into the excited spectrum. Let us first consider the integrable Hamiltonians H(1) and H(2), both corresponding to an effectively one-dimensional configuration space. When the critical point is reached in these systems, the global minimum of the potential V(1) or V(2) changes into a saddle point, which remains present for all values λ>λ (1) 0cor λ(2) 0c.The saddle point represents a singularity of (E,λ), causing the strongest type of ESQPT characterized by the infinite peak in the semiclassical level density. To see this, recall that in systems with f=1 the integral in Eq. (29) is equal to the period τof the single (uniquely determined) classical orbit at energy E; hence (E,λ)=τ(E,λ). If Ecoincides with the 033802-7 P. P ´ EREZ-FERN ´ ANDEZ et al. PHYSICAL REVIEW A 83, 033802 (2011) FIG. 4. Level dynamics for the SU(1,1) model (left) and for the SU(2) integrable model (right) with M(1) =M(2) =100 and ω =1. The scaled energies were obtained by an exact diagonalization. The ESQPT above the ground-state critical point λc0=0.707 is apparent in the bunching of levels around critical energies E(1) c=0.5andE(2) c= 0.25, respectively. energy of the x=0(ϑ=0) saddle point, the period becomes infinite because x=p=0 is a stationary point ( ˙ x=˙ p=0) of both Hamiltonians H(1)(p,x) and H(2)(p,x). Let us note that the same type of ESQPT is observed in systems with one quantum degree of freedom showing a local maximum of the potential, for instance in the Lipkin model [19,22,23] and many others [17,18]. A local increase of the level density at the saddle-point energy E(n) c=V(n) 0in both integrable models, i.e., the SU(1,1) and Jaynes-Cummings models (n=1 and 2, respectively), is demonstrated in Fig. 4. The two panels capture the evolution of quantum spectra for both models with the interaction parameter λ, showing clear indications of the ground-state QPT and its extension into the ESQPT on the right-hand side of the critical point, which for ω =1isatλ(1) c0=λ(2) c0=0.707 (see Table I). The calculation was done in a finite-size case, but it shows well-pronounced precursors of the phase transitional behavior. Additional ESQPT signatures are depicted in Fig. 5, which shows expectation values of operators proportional to K0= Nc+kand J0=Nc−jin individual excited states as a function of scaled energy Efor a fixed value of λ=1.5. Specifically, we consider the operators Na/M(1) =(2K0− 1 2)/M(1) for the SU(1,1) model and J0/j for the JaynesCummings model. Note that these two observables act as order parameters of the respective standard QPTs: their ground-state expectation values change from Na0=0toNa0>0, and from J00=jto J00<j,asλcrosses the critical point λc0. In both panels of Fig. 5, the respective system is well above the QPT critical point. The energy dependences of the respective expectation values NaEand J0Eshows a cusplike shape with a singularity localized at the ESQPT energy E(n) c (cf. Fig. 4). The two shapes are mutually reversed: while for the SU(1,1) model, the expectation value drops sharply to the lowest value at the critical energy, for the SU(2) model it has FIG. 5. Expectation values of Na(left) and Jz(right) for individual states across the spectrum at λ=1.5. The left and right panels, respectively, correspond to the SU(1,1) and the SU(2) integrable models with M(1) =M(2) =2000. The ESQPT is indicated by needlelike singularities located at the critical energies. a needle-shaped maximum. This is connected with a singular localization of the semiclassical wave function for E=E(n) cat the saddle point of the potential, i.e., at x=0 for the SU(1,1) and x=1 for SU(2) model (in both cases ϑ=0). This implies NcEc=0 for the SU(1,1) case (hence NaEc=0or1for even or odd systems, respectively) and NcEc=1 2M(2) for the SU(2) case (so J0Ec=1 2M(2) −j). An analogous effect (explained by infinite dwell times of a classical particle at the stationary point) is known from one-dimensional systems with a local maximum of the potential [40]. The ESQPT critical energies E(1) cand E(2) cdrop to the ground-state energy as λ decreases to the corresponding critical points λ(1) c0and λ(2) c0, and so do both cusp singularities in Fig. 5. Below the critical point the singularities disappear. For the nonintegrable Dicke model with the Hamiltonian H(3), the two-dimensional potential (24) has a saddle point at (x,y)=(0,0). This is connected with a nonanalytic dependence of the phase space volume (and the level density) on the scaled energy, although of a softer type than in the previous case, as follows from a higher dimensionality of the phase space for the Dicke model. Specifically, for λ>λ (3) c0 the level density exhibits an anomalous growth with an infinite derivative (singular tangent) at E=E(3) c, which coincides with the saddle-point energy V(3) cof the potential [18]. The resulting ESQPT is continuous, although its finite-size precursors very much resemble those of a first-order phase transition (the level density is close to a steplike function). The steplike increase of the level density in the Dicke model can be seen in Fig. 6, where the level dynamics with variable λis shown for j=2. Even for such a moderate value of the angular momentum, a sharp precursor of the ESQPT effect at absolute energy E(3) c=M(3)E(3) c=−1 is well visible in the spectrum above λ(3) c0=0.707 (for ω=ω0=1) as the lower interface between the horizontal and sloped level contours. Note that the effects of the Hilbert-space truncation (the cutoff for the number of photons; see Sec. II D) become relevant for the high-energy part of the spectrum. 033802-8 QUANTUM QUENCH INFLUENCED BY AN EXCITED-STATE ... PHYSICAL REVIEW A 83, 033802 (2011) FIG. 6. Level dynamics for the SU(2) nonintegrable model H(3) with j=2, obtained by a numerical diagonalization with Ntrunc ≈40. We show absolute energies in units of ω=ω0. A steep growth of the level density at the energy E=−2 [corresponding to the saddle point of the potential (24)] indicates a continuous ESQPT in the infinite-size limit. IV. QUENCH DYNAMICS A. Survival probability and energy distribution The Hamiltonians introduced in Sec. II have the common form (15), that is, H(λ)=H0+λHif the model-specifying superscript nis omitted. Here H0and Hrepresent the free and interaction terms, respectively, and λis a dimensionless control parameter. Let us stress that here we are working with the scaled Hamiltonian H=H/ℵ, but consider a finite-ℵcase, so that in general [H0,H]= 0. As seen from the expression H(λ2)=H(λ1)+Hwith =λ2−λ1, the above Hamiltonian allows one to apply perturbation techniques with the same perturbation Hfor all initial points λ1. Suppose that the system is initially prepared in one of the eigenstates |ψi(λ1)≡|ψ1of H(λ1)≡H1with energy E1(λ1)/ℵ≡E1. Below we will consider the initial state |ψ1coinciding with the ground state |ψ0(λ1),butthe formalism can be very easily developed for the general case. At time t=0, the value of the control parameter is abruptly changed from λ1to λ2=λ1+. The state |ψ1is no longer an eigenstate of the new Hamiltonian H(λ2)≡H2and starts evolving. The evolution after the quench can be monitored by the survival probability p1(t)=|a1(t)|2, where a1(t)=ψ1|e−iH2t|ψ1=|E2|ψ1|2   ω1(E2) e−iE2tdE2(30) is the amplitude describing the decay and recurrence of the initial state |ψ1for t>0. A formula of this form captures in general all quantum decay processes and has been studied in many different contexts (e.g., in analyses of the fidelity or Loschmidt echo [43]). Note that the use of the scaled Hamiltonian H2in Eq. (30) is equivalent to the t→t/ℵtransformation of time in the expression with unscaled Hamiltonian H2. Expanding the initial state |ψ1in the eigenbasis |E2i≡|Ei(λ2)of the Hamiltonian H2(with i=1,2,...enumerating discrete eigenvalues E2i), |ψ1= iE2i|ψ1   ci |E2i,(31) the survival probability reads as p1(t)= i|ci|4+2 i>j |ci|2|cj|2cos[(E2i−E2j)t].(32) As indicated in Eq. (30), the survival amplitude a1(t) can be written as the Fourier transform of the energy distribution ω1(E2)≡|E2|ψ1|2of the initial state in the eigenbasis of H2. The precise energy distribution is given by ω1(E2)≡ i|ci|2δ(E2−E2i).(33) Since both functions p1(t)inEq.(32) and ω1(E2)inEq.(33) are expressed in terms of the discrete energies E2iand the corresponding occupation probabilities |ci|2, they comprise fully equivalent information on the quench-induced relaxation process. The discrete form (33) of the energy distribution ω1(E2) can be approximated by its smoothed form ¯ω1(E2), obtained 033802-9