Non-equilibrium criticality in quench dynamics of infinite-range spin models
Paraj Titum Mohammad F. Maghrebi







Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742, USA
Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742, USA
Johns Hopkins University Applied Physics Laboratory, Laurel, Maryland 20723, USA

Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA





Long-range interacting spin systems are ubiquitous in physics and exhibit a variety of ground state disorder-to-order phase transitions. We consider a prototype of infinite-range interacting models known as the Lipkin-Meshkov-Glick (LMG) model describing the collective interaction of N spins, and investigate the dynamical properties of fluctuations and correlations after a sudden quench of the Hamiltonian. Specifically, we focus on critical quenches, where the initial state and/or the quench Hamiltonian are critical. 
Depending on the type of quench, we identify three distinct behaviors where both the short-time dynamics and the stationary state at long times are effectively thermal, quantum, and genuinely non-equilibrium, characterized by distinct universality classes and  static and dynamical critical exponents. These behaviors can be identified by an infrared effective temperature that is finite, zero, and infinite (the latter scaling with the system size as N^1/3), respectively. The quench dynamics is studied through a combination of exact numerics and analytical calculations utilizing the non-equilibrium Keldysh field theory. Our results are amenable to realization in experiments with trapped-ion experiments where long-range interactions naturally arise.



Non-equilibrium criticality in quench dynamics of infinite-range spin models
    Mohammad F. Maghrebi
    
============================================================================




The dynamics of isolated quantum systems has intrigued physicists since the dawn of quantum mechanics <cit.>. Furthermore, this topic has been in the spotlight in the past twenty years thanks to the experimental advances in ultracold atoms <cit.> and trapped ions <cit.> among others <cit.>. These platforms are some of the  prominent candidates for the quantum simulation of quantum phases of matter, but they are also well suited to investigate the dynamics away from equilibrium. A typical 
experimental setting is one where a system parameter suddenly changes – a scenario commonly described as a quantum quench. 

There is mounting evidence, both theoretical and experimental, that generic non-integrable systems thermalize upon a quantum quench and local correlations are best described by a finite-temperature ensemble <cit.>. On the other hand, integrable systems defined by an extensive set of conserved quantities fail to thermalize and instead are described by generalized Gibbs ensembles that also take into account all conserved quantities <cit.>. But, even integrable systems often thermalize in a weaker sense of thermalization if their long-wavelength properties are described by a finite effective temperature. For example, such effective thermal behavior has been identified in one-dimensional condensates <cit.> and even observed in experiments <cit.>; similar behavior is predicted in integrable O(N→∞) models <cit.>. This weaker notion of thermalization (with obvious merits for critical properties) is one that we adopt in this work.
A natural question is then if, upon a quantum quench and depending on the initial state,
even integrable systems always thermalize at long wavelengths, or, alternatively, can they exhibit genuinely non-equilibrium (critical) behavior?

In this manuscript, we consider the quench dynamics  of the prototypical Lipkin-Meshkov-Glick (LMG) model, an integrable model of spins with collective interactions.
We focus on the role of critical fluctuations and their universal properties. We show that, depending on the nature of the initial state (disordered or critical), distinct universal behaviors emerge in the dynamics. In particular, we show that fluctuations within the stationary state at late times can be described by an effective temperature which drastically depends on the initial state, and may vanish or even diverge for a quench from a critical state. We shall see that the latter divergence is a signature of a genuinely non-equilibrium critical behavior. 



We first emphasize what distinguishes our results from previous work. Conventionally, quench dynamics in the LMG model has been investigated with the spins initially in the ordered phase with the dynamics governed by mean-field equations  <cit.>. In contrast, we consider quenches where the order parameter is zero and the dynamics of fluctuations cannot be captured by a mean-field treatment. Quench dynamics near critical points have also been explored extensively in the context of the Kibble-Zurek mechanism <cit.> as the system is ramped through a quantum critical point <cit.>. However, this mechanism requires the initial and final states to be close to each other and a critical point setting it apart from our work.

Model.—
We consider a  prototypical model of infinite-range interactions known as the Lipkin-Meshkov-Glick (LMG) model <cit.> which was originally introduced in nuclear physics, but has also been used to describe other physical systems such as Bose-Einstein condensates <cit.>, small ferromagnetic particles <cit.> as well as  trapped ions <cit.>. The LMG Hamiltonian is given by

    H=-J/N∑_i<jγ_xσ_i^xσ_j^x+γ_yσ_i^yσ_j^y-Δ∑_iσ_i^z,

with the XY-like anisotropic interaction characterized by γ_x,y and the transverse field Δ. The collective nature of the model allows us to write the Hamiltonian in terms of the total spin operators S_a=1/2∑_i σ^a_i with a=x,y,z. Note that the Hamiltonian commutes with ^2=S_x^2+S_y^2+S_z^2, making it block diagonal in a basis defined by the total spin S. In fact, this model is integrable <cit.> and exactly solvable using Bethe ansatz <cit.>. The model of interacting spins is  schematically shown in <ref>(a).

The ground-state phase diagram of the LMG model [<ref>(b)] is well known, exhibiting a transition from a disordered paramagnet to an ordered ferromagnet  <cit.>. The in-plane magnetization serves as an order parameter where 1/N⟨ S_x⟩≠ 0 and/or 1/N⟨ S_y⟩≠ 0 in the ordered phase and 1/N⟨ S_x,y⟩=0 in the disordered phase, in the thermodynamic limit N→∞. The two phases are separated by a continuous transition along Δ=J   max{γ_x,γ_y}. Given the infinite-range interactions, the phase diagram can be obtained from a mean-field analysis. However, mean-field theory is insufficient where the order parameter is zero (outside the ordered phase) and particularly fails to capture the divergent fluctuations at the critical point. These fluctuations scale with N as <cit.>

    1/N⟨ S_x^2⟩∼ N^1/3,     1/N⟨ S_y^2⟩∼ N^-1/3,

along the critical line γ_x>γ_y; here, the prefactor 1/N factors out the trivial (square-root-volume) scaling away from criticality. Notice that the normalized spin fluctuations diverge only along the direction with the larger interaction strength. 
The exponent characterizing this divergence (1/3) is a distinct signature of the quantum phase transition <cit.>. 
For comparison, the fluctuations at a finite-temperature phase transition for this model diverge with a different critical exponent of 1/2 <cit.>. Therefore, the critical exponents distinguish between the quantum and thermal phase transitions.

Quench dynamics.—It is widely believed that the dynamics following a sudden quench leads to thermalization. The LMG model being integrable does not fully thermalize; nonetheless, we will see that a generic quench from the disordered phase gives rise to 
an effectively thermal behavior (including critical exponents) describing low-frequency modes. The question, however, remains if effective thermalization can be evaded at all, and, specifically, if a new, non-thermal scaling could emerge? Remarkably, the answer is in the affirmative. To show this, we study the dynamics for different types of quenches and initial states. To expose the critical behavior, the post-quench Hamiltonian is considered to be one at a critical point; without loss of generality, we take {γ_x=1,γ_y=0,Δ=J}. We consider three different initial states, each corresponding to the ground state of the LMG Hamiltonian but at different parameters: (i) : Initial state deep in the disordered phase, (ii) : Critical initial state on the critical line Δ=γ_x J, (iii) : Critical initial state on the critical line Δ=γ_y J; see <ref>(b).
Types II and III are distinguished by their initial divergent fluctuations in S_x and S_y, respectively.

It is instructive to first discuss a quench within  the disordered phase. In this regime, we can simplify the dynamics by using the Holstein-Primakoff approximation <cit.>, S_z=N/2-a^†a and S_x-iS_y ≈√(N) a. Defining a=1/√(2)(x+ip), we can write S_x≈√(N/2) x and S_y≈-√(N/2) p. The LMG Hamiltonian (<ref>) can be then cast as a harmonic oscillator with the frequency and the mass defined as Ω^2≡ 4Δ^2(1-Jγ_x/Δ)(1-Jγ_y/Δ) and 2m≡1/(Δ-γ_yJ), respectively; see the Supplemental Material (SM) <cit.>. In this picture, the quench can be viewed as a sudden change of the mass and frequency of the oscillator, {m_0,Ω_0}→{m,Ω}. It is straightforward to characterize the  fluctuations at long times in a quench to the vicinity of the critical point (mΩ≪ m_0Ω_0) <cit.>, 

    1/N⟨ S_x^2(t)⟩_t→∞=1/2⟨ x^2(t)⟩_t→∞≈m_0Ω_0/8m^2Ω^2.

This expression is reminiscent of a high-temperature harmonic oscillator (T≫Ω) where the equipartition theorem dictates 1/2 m Ω^2⟨ x^2⟩≈1/2 T which hints at the emergence of an effective temperature T_ eff=m_0Ω_0/4m <cit.>, or, equivalently,  

    T_ eff=Δ/2√(Δ_0-Jγ_x0/Δ_0-Jγ_y0).

Further insight can be obtained by examining the behavior of T_ eff for the different quenches. For the Type-I quench, the initial state is disordered (Δ_0> Jγ_x0,Jγ_y0) giving rise to a finite effective temperature. For a Type-II quench, the initial state is critical (Δ_0=Jγ_x0>Jγ_y0), resulting in a vanishing effective temperature. Most surprisingly, for the Type-III quench, the critical initial state  (Δ_0=Jγ_y0>Jγ_x0) leads to a divergent effective temperature.
This simple analysis hints at qualitatively different behaviors in Type I, II and III, which we will identify with an effective thermal, quantum, and non-equilibrium critical behavior, respectively. To this end, we shall go beyond the Holstein-Primakoff approximation both numerically using exact diagonalization and analytically via the Keldysh field theory. 
 
First, we introduce universal scaling functions that capture the dynamics of the correlations and fluctuations,


    1/N⟨ S_x^2(t)⟩   = N^α f(t/N^ζ),  
    
    C=1/2N⟨ [S_x(t_2),S_x(t_1)]_+⟩_ st.   = N^α C̃(t_2-t_1/N^ζ), 
    χ=1/2iN⟨ [S_x(t_2),S_x(t_1)]_-⟩_ st.   = N^ζ χ̃(t_2-t_1/N^ζ),


where f, C̃, χ̃ are scaling functions. 
The two-time correlators C and χ denote the correlation and response functions, respectively, which in equilibrium are related via the fluctuation-dissipation theorem <cit.>. The subscript (st) indicates the long-time limit (t_1,t_2≫ |t_1-t_2|) when a stationary state is approached. We introduce two critical exponents: ζ defines a dynamical exponent for the scaling of the characteristic time scale of the dynamics with system size, and α characterizing the scaling of fluctuations. Remarkably, we shall see that the same exponents describe the entire dynamics both at short times as well as the stationary state at long times.

Numerical results.—
Let us discuss the numerical results for the quench dynamics in the LMG model. The total-spin conservation allows us to simulate using exact diagonalization for large system sizes up to N=9000. For the numerical simulations, we set J=1, and restrict ourselves to the largest spin sector S=N/2. 
Beside the post-quench Hamiltonian with {γ_x=1,γ_y=0, Δ=1}, the initial states correspond to the ground state of the Hamiltonian with the following parameters: (i) Type-I: {γ_x=1,γ_y=0, Δ=4}, (ii) Type-II: {γ_x=1,γ_y=0.5, Δ=1}, (iii) Type-III: {γ_x=0,γ_y=1, Δ=1}. 



The fluctuations in the order parameter are shown in <ref> where each row corresponds to a given quench type. In the case of the Type-I quench [Figs. <ref>(a) and (b)], the initial state is in the disordered phase with small (i.e., non-critical) fluctuations. Fluctuations grow initially (t≲1/J) independent of the system size, but peak at longer times that increase with the system size. 
The scaling collapse of the different curves [Fig. <ref>(b)] indicates that fluctuations diverge as N^0.5 and furthermore evolve with a characteristic time scale ∼ N^0.25 before reaching the stationary state; hence, we identify the exponents ζ=0.25 and α=0.5. Indeed, the same exponents govern the two-time correlators in the stationary state consistent with the scaling in Eq. (<ref>); see the SM <cit.>. Interestingly, these exponents are identical to those governing a thermal critical point <cit.>. This might be surprising because a true thermal phase transition only occurs at Δ/J<1 <cit.> in contrast with  Δ/J=1 chosen above, hence the stationary state cannot be described as a Gibbs state <cit.>.
However, the fact that the critical behavior is consistent with a thermal phase transition hints at an effective thermalization at low frequencies. 

For the Type-II quench, the initial state is critical with divergent fluctuations in S_x. As shown in Fig. <ref>(c), fluctuations do not significantly grow over time. This observation indicates that the sudden quench has only slightly disturbed the system. Indeed, the scaling collapse shown in Fig. <ref>(d) reveals the exponents ζ=0.33 and α=0.33, consistent with a quantum critical behavior already present in the initial state. The emergence of quantum criticality in a quench from a critical state is also observed in the O(N) model at large N <cit.>.

Most interestingly, the Type-III quench exhibits novel non-equilibrium behavior that is neither thermal nor quantum critical. As shown in <ref>(e), fluctuations grow faster than the other quenches. While the dynamics might seem similar to the Type-I quench, the exponents are markedly different: (ζ,α)=(0.17,0.67). This indicates that fluctuations diverge with the system size even more strongly than those in Type-I or, equivalently, at the thermal critical point. Indeed, as we shall see shortly, the effective temperature in this case diverges with the system size. We also compute the two-time correlators for this quench as shown in <ref>(a). Again, we find that the correlation and response functions obey the scaling forms in <ref> with approximately the same critical exponents.



Inspired by the fluctuation-dissipation theorem, we can also identify an effective temperature describing the stationary state. In equilibrium and at low frequencies (relative to temperature), this theorem relates fluctuations and dissipation as C(ω)=(2T/ω) iχ(ω) <cit.>. This equation can also be recast in the time domain,  χ(t)=1/2T∂_t C(t) with t the time difference. Away from equilibrium, this relation is commonly used to identify a time-dependent effective temperature <cit.>. Indeed, for the Type-III quench, the two-time correlators in <ref>(a) satisfy the latter equation for an approximately constant T(t); see <ref>(b). In <ref>(c), it is clear that the average temperature, T_ av=1/t∫_0^t dt̃ T(t̃), scales with the system size as T_ av∼ N^1/3, and indeed diverges with the system size. 
This should be contrasted with the constant effective temperature for the Type-I quench <cit.>.We note that T_ av is technically different from T_ eff(ω=0); the latter is evaluated from the IR limit of the ratio of correlation and response functions in the frequency domain which is difficult to access numerically. Nevertheless, we expect the scaling behavior to be the same. 

Scaling analysis.—
The scaling behavior observed in the numerics hints towards an underlying scaling theory. In fact, this behavior and the scaling relations postulated in <ref> can be obtained from an effective low-frequency theory. At (or near) the critical point of the post-quench Hamiltonian, the relevant degree of freedom can be identified from the Holstein-Primakoff bosons as x=(a+a^†)/√(2), which also serves as the order parameter within the ordered phase. The non-equilibrium nature of the model necessitates the Keldysh formalism at the expense of doubling the fields x→ x_c/q that define the Keldysh basis. The non-equilibrium partition function is given by <cit.>

    Z=∫ [dx_c,q] W(x_c0, Kẋ_c0)   e^iS_K.

Here, W(x_c,p_c) is the Wigner function describing the initial state (hence the variables with the subscript 0), and S_K is the Keldysh action that is given by S_K= -∫_0^∞ dt [K x_q ẍ_c+ r x_q x_c +u/N(x_c^2+x_q^2)x_cx_q]. The coefficients in the action are given in terms of the microscopic parameters as 
K^-1≡ 2(Δ-Jγ_y), r≡ 2( Δ-Jγ_x) and u≡ J γ_x; see the SM <cit.>. Notice that r= 0 defines the critical point after the quench. The Wigner function describing the initial state can be written as W(x,p)= W(x^2N^-α_0,p^2N^α_0) where α_0=0, 1/3,-1/3 for Type I, II and III, respectively. This correctly produces ⟨ x^2⟩∼ N^α_0 (similarly for p) in the initial state; see <ref>. The precise form of the Wigner function as well as the coefficients in the action are not required for the scaling analysis; for convenience, we set K=1.



We shall focus on Type I and III where the system is significantly disturbed upon quench. 
At the critical point, r=0, we expect (confirmed by numerics) that the fluctuations diverge, that is, 1/N⟨ S_x^2⟩∼⟨ x_c^2 ⟩∼ N^α for an exponent α(>0) to be determined; we then identify the scaling dimension [x_c]=α/2. This divergence can be scaled away by defining X_c=x_c N^-α/2. Anticipating (again confirmed by numerics) a critical slowdown, we also introduce a rescaled time T =t N^-ζ with a dynamical exponent ζ. By appropriately choosing these exponents, the action together with the Wigner function can be made scale invariant, i.e., independent of N. Writing the Wigner function in terms of the rescaled variables, we have W(x_c0,  ẋ_c0)= W(X_c0^2 N^α-α_0,  X_c0'^2 N^α-2ζ+α_0) with X'=dX/dT. 
Now, the function W is only significant when its arguments are of order 1. Therefore, fluctuations of X_c0 are greatly suppressed for Type-I and III quenches (α_0≤ 0), thus imposing X_c0≈ 0. On the other hand, fluctuations of X_c0' can be made scale invariant by setting α-2ζ+α_0=0. Similar analysis can be made for the terms in the action. The kinetic term (∼∫_t x_q ẍ_c) can be made scale-invariant by introducing the new rescaled variable  X_q=N^-α/2+ζ x_q, hence the scaling dimension [x_q]=-α/2+ζ. Next, we turn to the interaction term, (-1/N)∫_t u_c x_q x_c^3+u_q x_q^3 x_c, where “classical” (∼ u_c) and “quantum” (∼ u_q) vertices are distinguished although u_c=u_q at the microscopic level. Let us first cast the classical vertex in terms of the rescaled variables, upon which u_c/N → u_c N^2ζ+α-1. Making the latter scale invariant yields 2ζ+α-1=0. Combining the above relations, the critical exponents are determined as ζ=1+α_0/4 and α=1-α_0/2. 
One can then confirm that, upon the scaling transformation, the quantum vertex is suppressed by a power of N and thus can be neglected; in the language of renormalization group theory, the quantum vertex is irrelevant. This itself is a consequence of classical and quantum fields having different scaling dimensions ([x_c]=1-α_0/4 vs [x_q]=α_0/2). Our scaling analysis yields the exponents reported in  Table <ref> consistent with numerics, and also reproduces the behavior of the effective temperature (see the SM <cit.>). We remark that, in analogy with boundary critical phenomena <cit.>, our analysis has relied on scaling both “boundary” and “bulk” terms. 

While we have focused on finite-size scaling, we can just as well identify the scaling behavior in the thermodynamic limit (N→∞) away from the critical point (r>0). Interestingly, we find that the fluctuations at late times diverge as 1/r, 1/√(r) and 1/r^2 for Type I, II and III, respectively; see the SM <cit.>. While the first two scalings can be identified with the characteristic quantum and thermal behavior, the distinct scaling for Type III is indicative of genuinely non-equilibrium critical behavior.

The non-equilibrium dynamics reported in this work is accessible in a variety of experimental platforms, particularly in the context of trapped ions <cit.>. These systems are described by spin models with long- or even infinite-range interactions. 
A challenge is to prepare an initial critical state for Type-II and III quench dynamics. However, based on a quantum approximate optimization protocol <cit.>, variational quantum algorithms have been recently proposed <cit.> and implemented <cit.> to efficiently prepare quantum critical states.

Conclusion and outlook—
In this work, we  have studied quench dynamics in the context of infinite-range spin models. We have identified a rich dynamical universal behavior depending on the critical nature of the initial state. Specifically, we have identified a scenario (Type III) where a genuinely non-equilibrium critical behavior emerges both in the short-time dynamics and long-time stationary state. Our findings complement the distinct signatures of correlated initial conditions on the dynamics of classical systems <cit.>. 
An important future direction is to identify efficient routes to preparing critical initial states in such models. An intriguing alternative is to first perform a Type-I quench, whose initial state can be easily prepared, that leads to a (thermal-like) critical state at late times. Upon a second quench of Type III, the system then approaches another critical, yet non-equilibrium state at late times. Another particularly relevant direction is to extend our results to long-range spin models, V(r)∼ 1/r^p, specifically with p<1 where long-lived prethermal states emerge <cit.>. More generally, extending our results to other integrable models or those exhibiting long-lived prethermalization is worthwhile. 
Finally, it would be interesting to identify any ageing behavior <cit.> in such models. 


Acknowledgments.—
We thank Marcos Rigol and Daniel Paz for useful discussions. P.T. acknowledges support from the NIST NRC Research Postdoctoral Associateship Award. M.M. acknowledges support
from NSF under Grant No. DMR-1912799 and start-up funding from Michigan State University. This work is supported in  part  by  the  U.S.  Department  of  Energy (DOE),  Office  of  Science,  Office  of  Advanced  Scientific Computing Research (ASCR) Quantum Computing Application Teams program, under fieldwork proposal number ERKJ347. The authors acknowledge the University of Maryland supercomputing resources (http://hpcc.umd.edu) made available for conducting the research reported in this Letter. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.


apsrev4-1


in 1,...,7

[pages=]draftLMG-supp 


