{'fraction_non_alphanumeric': 0.03977091950365892, 'fraction_numerical': 0.009619883588179147, 'mean_word_length': 4.243572129538763, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical transitions and quantum quenches in mean-field models

We develop a generic method to compute the dynamics induced by quenches in completely connected quantum systems. These models are expected to provide a mean-field description at least of the short time dynamics of finite dimensional system. We apply our method to the Bose-Hubbard model, to a generalized Jaynes-Cummings model, and to the Ising model in a transverse field. We find that the quantum evolution can be mapped onto a classical effective dynamics, which involves only a few intensive observables. For some special parameters of the quench, peculiar dynamical transitions occur. They result from singularities of the classical effective dynamics and are reminiscent of the transition recently found in the fermionic Hubbard model. Finally, we discuss the generality of our results and possible extensions.


Introduction
Thanks to the fast experimental progress on cold atoms [1,2,3], where direct control over the parameters of the effective Hamiltonian is available, many theoretical questions about out of equilibrium dynamics of quantum systems have been raised and started to be addressed. Among them, we cite the problem of thermalization for an isolated system (see e.g. [4,5]) the existence of long-lived out of equilibrium states (see e.g. [6,7,8,9]), dynamical transitions out of equilibrium [10,11,12,13], etc.
In this work we focus on the off-equilibrium dynamics induced by quantum quenches, i.e. the evolution of an isolated quantum system after a sudden change of a parameter of the Hamiltonian. This problem has been intensively studied these last years. It provides an useful idealization of the protocols followed in experiments, where the change of parameters takes instead place at finite rates. The literature on quantum quenches is already quite broad, see for example the reviews [14,15,16]. One dimensional systems have been intensively studied by exact analytical and numerical methods. Results for higher dimensional systems are instead scarcer; the previous methods cannot be applied and one has to resort to approximations of some kind. The ones that have been more used are the Bogoliubov method [17], path integral saddlepoint expansions [18,19], the large-N limit [20] and mean-field (or large dimension) approximations [21,10,11,12,22]. Interestingly, mean-field approaches revealed that quantum quenches may lead to out of equilibrium dynamical transitions. This was first found in the analysis of the Fermionic Hubbard model by time dependent Dynamical Mean Field Theory (t-DMFT) in [10]. Later, a confirmation and explanation of this result was found by using a Gutzwiller time-dependent Ansatz in [11]. In [12] we solved the Bose Hubbard model on a completely connected lattice and found a very similar dynamical transition at integer fillings. This was also noticed in the hard-core Bose-Hubbard model in a superlattice potential [23]. Finally, clues to argue that these dynamical transitions are generic even beyond mean-field were recently presented in [13] by mapping the problem to the one of classical phase transitions in films.
In this article we first describe a generic approach to solve the dynamics of quenches in completely connected quantum models. Then, we apply our method to the Bose-Hubbard model, the generalized Jaynes-Cummings model and the Ising model in a transverse field. We find that the dynamical transition, discovered for the Hubbard model, is a systematic dynamical effect present in completely connected systems characterized by a quantum phase transition in equilibrium. A shorter version of this work, that focused only on the Bose-Hubbard model, appeared in [12].

Summary of the method and results
In this article, we focus on the out of equilibrium dynamics in several completely connected models: the Bose Hubbard model (BHM), the transverse field Ising model (IM) and the Jaynes-Cumming model (JCM). The reason for focusing on systems defined on completely connected graphs is that this allows one to study a well defined model and-at the same time-to obtain an approximate solution for finite dimensional lattices. Indeed, such models are related to several approximations used in the literature. Among the most known ones, we cite the limit of infinite dimensions for models defined on a hyper-cubic d-dimensional lattice ‡ and the Gutzwiller Ansatz, a widespread method to obtain mean field equations. Actually, in the case of completely connected models, the Gutzwiller Ansatz can be shown to be exact by using a Hubbard-Stratonovich transformation.
We shall focus on out of equilibrium dynamics induced by a quantum quench. This procedure is defined as follows. LetĤ(λ) be the quantum Hamiltonian, which depends on a coupling λ. The system is prepared in its ground state |ψ(t < 0) = |GS(λ i ) at some coupling λ i . At t ≥ 0, the system is driven out of equilibrium in a controlled fashion, tuning the coupling λ(t) in time according to a predefined procedure. Except for quasistatic procedures, often called "adiabatic", the wave function becomes different from the ground state |ψ(t > 0) = |GS(λ(t)) . We call sudden quench the procedure in which the coupling is suddenly switched to a final value λ f : λ(t) = (λ f − λ i ) θ(t) + λ i .
As we shall show, the quantum dynamics in a completely connected system can be solved by mapping the unitary evolution onto an effective model undergoing Newtonian dynamics. This is a drastic simplification which is possible thanks to the symmetry of the completely connected Hamiltonian under any permutation of sites. We restrict our analysis to cases where |ψ(t = 0) is the ground state at some coupling λ i . Thus, |ψ(t = 0) is also symmetric under permutation of sites and since both the initial state and the Hamiltonian are symmetric, the generic unitary evolution takes place in the sector of symmetric states only. In this subspace, the states are parametrized using a few local macroscopic observables. The unitary evolution can then be written as a Schrödinger equation in the symmetric space, which involve an effective = V −1 , where V is the number of sites of the system. Thanks to this property, at the thermodynamic limit, the entire dynamics of the system can be encoded in the one of few macroscopic variables, which undergo an effective classical Hamiltonian evolution.
We first apply this approach to the Bose-Hubbard model, our main interest because of its applicability to cold atom experiments. In this case, it is the onsite repulsion U that plays the role of the coupling λ. In a first stage we make the additional assumption that the number of bosons per site is less or equal to n max b = 2. The average superfluid order after the quench is a non-monotonous function of |U f − U i |. Actually, it decays logarithmically to zero at some special values of U i and U f . For these special quenches, the superfluid order relaxes exponentially to zero. We call this peculiar feature a dynamical transition. Surprisingly, the microcanonical equilibrium characterized by same energy, towards which the system would relax on large times if it were able to thermalize §, has non-zero superfluid order. Thus, this transition is a purely dynamical phenomenon, which has nothing to do with the fact that the superfluid order vanishes in the high temperature phase. We extend our analysis to n max b ≥ 3. In this case the effective dynamics involves n max b − 1 degrees of freedom, and the classical trajectories can be either regular or chaotic. Apart from this difference, the dynamical transition is found to be qualitatively unchanged and the quantitative differences with n max b = 2 are small.
Finally, we focus on two other systems, for which some out of equilibrium properties have already been studied, a generalized Jaynes-Cummings model [25], and the Ising model in a transverse field [21]. We apply our formalism to both systems and show that also in these cases there is a dynamical transition for sudden quenches.
This paper is organized as follows: In section 3, we briefly describe the three considered models. In Section 4, we describe how the effective classical dynamics can be derived for an arbitrary Hamiltonian. Then in section 5, we derive the classical dynamics for the Bose-Hubbard model with truncation n max b = 2, and analyze quenches and the dynamical transition in section 6. The general case n max b ≥ 3 is considered in section 7. The effective dynamics and quench properties of the generalized Jaynes-Cummings model are described in section 8, and in section 9 for the Ising model. Section 10 and 11 contain discussions and possible extensions of our work. Appendix A is devoted to prove that the same effective dynamics is recovered with a Gutzwiller Ansatz wave function, and in Appendix B the WKB eigenstates of the completely connected model are discussed.

Bose-Hubbard model and the Mott-Superfluid transition
The Bose-Hubbard lattice model has regained a lot of interest recently since it can be realized and studied in experiments on cold atoms [26]. The Hamiltonian of its completely connected version, suited to describe the limit d → ∞ of the lattice model, reads where b † i , b i are the bosonic creation and annihilation operators, The first term is an on-site repulsion between bosons, the second is a tunneling term, of amplitude J, and rescaled by the volume V (number of sites). This rescaling is needed to obtain a well-defined thermodynamic limit V → ∞. The limit of infinite dimensions that is equivalent to this model is obtained by taking a coupling J 2d for a BHM defined on a d-dimensional hyper-cubic lattice.
This model undergoes a quantum phase transition at commensurate fillings: n i = n with n integer. See [27] for a detailed discussion of the phase diagram. Notice that to the Gibbs ensemble, as we shall discuss later. This implies that in the large dimension limit the timescale for equilibration diverges as a function of d.
since we consider the dynamical behavior, the number of particles is fixed because it corresponds to a quantity conserved by the dynamics. In consequence we shall discuss the phase diagram in terms of the density and not the chemical potential, as it is instead done usually. Let us focus first on cases where the number of particles per site is an integer. To grasp why there is a quantum phase transition, it is helpful to compare the ground state in the two limits U → ∞ (or J = 0) and U → 0. In the first case, the ground state is diagonal in the occupation number |ψ = ⊗ i |n i = n . It remains the same for all U < U c , and is called a Mott insulator. It is incompressible because adding or removing a boson requires an energy of the order of U. In the second case, U = 0, the ground state is akin to a product of coherent states |ψ = C i e −αb † i |0 . For U = 0, and also moderate values of U, the ground state is superfluid and compressible. The usual order parameter b i is zero because the density is fixed, hence the appropriate parameter is the off-diagonal long range order measured by |Ψ 0 | 2 = lim d ij →∞ b † jb i . In the completely connected model, all distances d ij between two different sites are equal to 1. This is the largest distance in the problem. In consequence one can define ¶: Because the system is completely connected, phonons are absent and the spectrum always has a gap. Beyond a critical coupling U c , the order parameter |Ψ 0 | 2 vanishes and the ground state becomes a Mott insulator. This is the Mott-superfluid quantum phase transition. For non-integer fillings, the system is instead always superfluid and compressible as it can be understood considering perturbation around the U → ∞ limit. At non zero temperature, there is a second order phase transition from the superfluid phase to a Bose gas at a critical temperature T c (U).

Generalized Jaynes-Cummings model and superradiance transition
The second model that we consider is the generalized Jaynes-Cummings model studied in [25], which describes N distinguishable two-level systems in interaction with a single quantized electromagnetic mode. It is useful in various contexts, and has been suggested as a qualitative description of the BEC/BCS crossover in fermionic condensates, of molecular magnetism and of the formation of dimers by pairing of two bosons. Recently, new experiments in the setting of cavity quantum electrodynamics [28] are realizations of the Dicke Hamiltonian [29], which can be mapped onto the Jaynes-Cummings Hamiltonian if the "rotating wave approximation" is made. The dynamics of sweeps, starting from an empty bosonic mode, was studied in [25] using quasiclassical approximations and the truncated Wigner approximation. Here, we describe the dynamics of sudden quenches from the broken symmetry phase.
Formally, this state is the true ground state only in the grand canonical ensemble. However, the canonical and grand canonical ensemble are equivalent up to 1/V corrections. ¶ In equilibrium, it is easy to check that the definition below gives back the usual one used in the grand-canonical ensemble, The model is defined as follows. There is one bosonic [b,b † ] = 1 degree of freedomthe electromagnetic mode-and one SU(2) spin S = N/2. The spin degree of freedom keeps track of the number of excited two-level systems, n exc = S z + N/2. Withn =b †b , the Hamiltonian iŝ It includes a potential energy for bosons (ω), for the spin (ω 0 ), and a coupling g between the two. Each time a boson is lost, a two-level system is excited, and reciprocally: the quantityQ =Ŝ z +b †b is conserved by the dynamics. It is thus possible to define 2λ = ω 0 − ω such that up to a constant term one finds: In the following we take g = 1 and consider λ in units of g, and the time in units of /g. We consider the regime Q ≥ N/2 (which can be reduced to Q = N/2), for which there is a quantum phase transition. At T = 0, the system is in the normal ground state n = 0 if λ < λ c = −1 and is in the super-radiant ground state for λ > λ c where n = 0. This is the super-radiance quantum phase transition. This transition is the consequence of the competition between the two terms in the Hamiltonian: the first one favors removing bosons (if λ < 0), whereas the second plays the role of a kinetic energy and it can lower the energy if the bosons modes are filled, see section 8. As for the Bose-Hubbard model, there is a finite temperature phase transition at T c (λ), above which the super-radiant phase disappears. Notice that this phase transition is usually studied in the canonical/microcanonical ensemble, whereas here it is more natural to focus on given value of Q since this quantity conserved by the dynamics.

Ising model in a transverse field and the ferromagnetic transition
The Ising model in a transverse field is a paradigm of quantum phase transitions [30]. The Hamiltonian of its completely connected version readŝ where J is the ferromagnetic coupling and Γ the transverse field. For simplicity, we set J = 1, measure Γ in units of J and time in units of /J. The Hamiltonian can be written using a single spin Ŝ = i Ŝ i and readŝ The large N limit corresponds to the large spin S = N/2 limit, which is also the classical limit. The ground state is found, minimizing the corresponding classical Hamiltonian + Ŝ → S = S{sin θ cos φ, sin θ sin φ, cos θ} [21]. The minimum is at φ = 0, and for Γ < 1/2, sin θ = 2Γ, which is a ferromagnetic ground state, in which the spins align toward + The kinetic term of the classical limit is nontrivial to obtain, but we do not need it for the moment.
the z-axis. For Γ > 1/2, θ = π and the ground state is called quantum paramagnet, because it is unoriented in the z-basis |± i , |ψ = 1 The quantum phase transition from z-ferromagnet to quantum paramagnet takes place at Γ c = 1/2. Starting from the ferromagnetic ground state and increasing the temperature the system undergoes a second-order phase transition at a finite temperature T c (Γ).
4. Quantum quenches in completely connected models: mapping to an effective classical dynamics

Site permutation symmetry and classical Hamiltonian dynamics
In the following, we show how the dynamics induced by quantum quenches, in arbitrary completely connected models, can be mapped onto an effective classical Hamiltonian dynamics. As already stated, any model defined on a completely connected lattice has a Hamiltonian which is symmetric under any permutation of sites. Thus, one expects that the ground state of the model is also site permutation symmetric * . We shall study quenches starting from the ground state of the system. Since the site permutation symmetry is conserved by the unitary evolution, the dynamics only takes place into the subspace of symmetric wave functions. This is a drastic simplification, because symmetric states can be described using a few variables only.
To clarify this point in a concrete but simple case, let us see what these variables are for the Bose-Hubbard model, with the additional constraint that there are n max b = 2 or less bosons per site. First, remark that given a particular Fock state |{n i } , the linear combination of every possible site permuted Fock state is a site permutation symmetric state |S({n i }) ∼ P |{n P (i) } . This new state is completely characterized by the fraction x 0 , x 1 , x 2 of sites with 0, 1 and 2 bosons respectively, and we call it |x 0 , x 1 , x 2 = |x where x is a shorthand for all variables. Since the Fock states form a basis of generic states, the {|x } states form a basis of the symmetric sector. In order to express the Schrödinger evolution in this basis, we have to compute the transitions elements . .} a vector of integers. For example in the Bose-Hubbard model, the matrix elements of the Hamiltonian connect states that differs by one boson jump. Thus x 0 , x 1 and x 2 can only differ by 1/V or 2/V . For example, after a jump Thanks to the "locality" of transition elements between symmetric states, the Schrödinger equation in the symmetric basis |ψ = x ψ x (t)|x takes a simple form.
The Schrödinger equation projected on x| reads x|i∂ t |ψ = x|H|ψ . Denoting the diagonal transition element where we made use of ψ x+m/V = exp(m i ∂ x i /V )ψ x with an implicit summation over the index i. Strikingly, the Schrödinger equation (7) involves an effective = 1/V , thus the regime of interest V → ∞ corresponds to the classical regime. Furthermore, for the ground state, the initial wave function is a narrow wave packet of width 1/ √ V (this general property is discussed in more detail for the Bose-Hubbard model in section 5.2). Thus, by the Heisenberg uncertainty principle, we expect that momentum fluctuations are proportional to √ ∝ 1/V . Therefore the evolution of the wave-packet can be fully described in the thermodynamic limit by its average (or center) x(t) = x and average momentum p(t) = p . Both quantities evolve following a classical Hamiltonian dynamics obtained from the quantum one by replacingp The classical Hamiltonian evolution of the variablesẋ i (t) = ∂H/∂p i andṗ i (t) = −∂H/∂x i yields the evolution of the wave packet after the quantum quench, and give access to all observables as a function of time. A more careful analysis is performed in section 5.3 and fully supports this picture. As a conclusion, the analysis of quench dynamics in connected models is tractable, and takes the form of an effective classical Hamiltonian evolution. The initial condition is provided by the ground state values of x i , p i at the initial coupling U i . The sudden quantum quench dynamics is then described by the classical dynamics of x(t), p(t) induced by the effective Hamiltonian H(U f ), with initial conditions x i , p i . In the following, we study the Bose-Hubbard model with two or less bosons per site. We derive explicitly the effective classical Hamiltonian in section 5.1. The nature of the ground state, in particular the peaking of the wave packet when V → ∞ and the quantum phase transition, are discussed in section 5.2. A more thorough derivation of the effective Hamiltonian dynamics is presented in section 5.3. Finally, we conclude by describing the phase space of effective trajectories in section 5.4.

Schrödinger equation in the site permutation symmetric basis and effective Hamiltonian
In this section, the method sketched above is applied to the Bose-Hubbard model with truncation n b ≤ n max b = 2. This restriction is well-suited to derive analytical expressions, and preserves the main features of the model, such as the Mott insulator to superfluid quantum phase transition. The case of a larger number of bosons per site is considered in section 7. The results remain qualitatively the same.
The symmetric states are written under the form |x 0 , x 1 , x 2 . Since the total number of sites V = V (x 0 +x 1 +x 2 ) is fixed, and the overall density of bosons n = N/V = x 1 +2x 2 is conserved by the dynamics, the symmetric states can be labeled by one variable only. Choosing x 1 as such a variable, the symmetric normalized wave functions are denoted and ′ means the sum over Fock states of fixed fraction x 1 . In the following for simplicity of notation we drop the subindex one and use x instead of x 1 . Note that contrary to the previous paragraph, x is now just a number and not a vector.
To compute the transition rates, we proceed as follows.
The repulsion term Using the previous notations, the only possible transitions correspond to m = 2, m = −2, m = 0 representing respectively x + 2/V |Ĥ|x , x − 2/V |Ĥ|x and x|Ĥ|x . The first amplitude reads There are V !/(N 0 !N 1 !N 2 !) factors on the ket side. For each of these factors, there are N 0 N 2 transition to a state with N ′ 1 = N 1 + 2. Therefore using that and after a partial cancellation of the normalization factors, we obtain: In the following, only dominant contributions of order V are kept, the sub-leading ones are dropped since they do not matter for dynamics taking place on times not diverging with the system size. Rewriting x 0 and x 2 as functions of x, the transition rates are The Schrödinger evolution on ψ x (t) in the site permutation symmetrical basis is Therefore we have obtained that the Hamiltonian in the subspace of symmetric states readsĤ Note that the dimension of the Hilbert space has been reduced from e V to V in this case. The case n max b = 2 is especially convenient because the effective classical motion takes place in one dimension, therefore the equation of motion are integrable and the evolution easy to understand and describe.

Nature of the ground state and quantum phase transition
As in section 4, the Schrödinger equation (13) involves an effective = 1/V , and its classical limit is obtained throughx → x(t) andp → p(t).
The possibility of describing the quantum dynamics after a quench in terms of classical dynamics relies on the fact that the initial wave function is initially of small width (∼ 1/ √ V ). In order to check this, let's first expand the eigenvalue equation at lowest order in 1/V : Since the kinetic term has a factor 1/V 2 in front, the ground state corresponds to the minimum of D x − 2W x and small quantum fluctuations around it. Indeed, around the minimum, x GS , a quadratic expansion of W x and D x maps this problem onto a quantum harmonic oscillator with Thus, we find that the ground state wavefunction is centered around the absolute minimum x GS of the potential D(x) − 2W (x) and has a width σ ∼ /mω, which is of the order of 1/ √ V . An appropriate expression for any eigenstate is provided by the WKB expansion is given in section Appendix B.
The quantum phase transition taking place as a function of U can be recovered within this formalism as we now show. We recall that for the Bose-Hubbard model, the quantum phase transition and the Mott phase are present at commensurate fillings only, so we have to focus on n = 1 for n max b = 2. The ground state of the quantum Hamiltonian corresponds to the global minimum of the effective Hamiltonian (15), given by ∂H ∂x = ∂H ∂p = 0. These equations lead to p = 0 and x equal to the value x GS , which verifies ∂(D(x)−2W (x)) ∂x = 0. It is easy to check that where U c = 3 + 2 √ 2 ≃ 5.82843 (note that x ≤ 1). In the thermodynamic limit, the Mott insulator ground state is |n 0 = 1, n 1 = 1 . . . = |x = 1 . Instead the superfluid ground state is |x = x 0 , with 1/2 < x 0 < 1. Therefore, in the superfluid state, a fraction of sites proportional to the intensive kinetic energy, and can be computed using the identity It is the order parameter for the transition: it is positive in the superfluid phase and vanishes in the Mott insulating one.

Effective classical evolution of wave packets
In the following we show in detail how the classical Hamiltonian dynamics emerges from the quantum evolution of symmetric states in the thermodynamic limit. For simplicity we consider the situation of section 5, where n max b ≤ 2, but the argument is more general. The Schrödinger equation reads A ground state wave function is a wave packet characterized by a small width of the order of 1/ √ V as shown in section 5.2. The width broadens after a long time, which diverges with V . On times that do not diverge in the thermodynamic limit, the ground state wave function and its subsequent evolution can be written as with f (x, t) = g(x, t) − iθ(x, t), g and θ being respectively the envelope and the phase of the wave-packet. The envelope g(x, t) has a maximum at x(t), and expanding g(x, t) around x(t) shows that indeed this wave function describes a packet of width 1/ √ V . After these preliminary considerations, we proceed and evaluate the evolution of the wave function, plugging the expression (18) of ψ(x, t) into (17): For the amplitude and the phase, this leads to: The position of the peak x(t) corresponds to the maximum of the amplitude, thus it satisfies the implicit equation ∂ x g| x(t) = 0. To obtain its evolution, we differentiate it with respect to time: The term ∂ t ∂ x g above can be computed from (20). At the point where ∂ x g| x(t) = 0, (21) becomes: This equation involves only one unknown quantity, ∂ x θ| x(t) , which, as we shall show, is akin to a momentum. A self-consistent equation for this quantity can be found taking its derivative with respect to time: In the previous expression, ∂ t ∂ x θ| x(t) is the space derivative of (20). Using (22) the two terms in ∂ 2 x θ| x(t) cancel out, and (23) reads: This, together with eq. (22), provides a set of closed differential equations for x(t) and ∂ x θ| x(t) . The last thing we have to show is that the quantity ∂ x θ| x(t) is the average value of the momentum operatorp, which reads: Because of the form (18) of ψ x (t) these integrals can be performed by the saddle point method (in the limit V → ∞). At the saddle point, ∂ x g(x(t), t) = 0 and ∂ x θ(x(t), t) has a non zero value, thus p = ∂ x θ| x(t) +O(1/V ). The same is true for x = x(t)+O(1/V ). Thus, rewriting (22) and (24), one finds that the phase p(t) and the position x(t) obey the equations: As a conclusion, the average values x(t) and p(t) of the position and momentum operator, which describe the global behavior of the wave packet, have the same time evolution as classical canonical variables with Hamiltonian H[x, p] = D x − 2W x cos(2p). This property fits into the general picture of the propagation of wave packets in the semiclassical regime (see for example [31] for a mathematically oriented review) and, hence, was expected on general grounds as discussed previously. The generalization to more than one variable is straightforward. The width of the packet can also be evaluated.
It is related to the separation in time of two classical trajectories initially separated by the width of the packet. When the effective dynamics is one dimensional, the classical Hamiltonian is integrable and periodic orbits separate linearly in time. In this case, since trajectories start from an initial distance 1/ √ V , the typical time of separation is t ∼ √ V . We check numerically that our analysis of the limit of large V is correct. In figure 1, the exact quantum evolution, obtained by diagonalization of the discrete equation (7), is compared to the classical evolution for short times. We find an excellent agreement.

Different effective trajectories in the phase space
We now study the effective classical dynamics by focusing on the phase space properties.
Since the Hamiltonian (15) is one-dimensional, and energy is conserved there are lots of constraints on the motion. In particular all the trajectories are integrable, and can only be periodic in x(t) except for separatrix trajectories. We find that the phase space is divided in two regions by a separatrix. To see this, one can characterize trajectories in terms of their turning points, and considering whether the momentum p is bounded or not. There are three types of turning points: the derivative dx dt = 4W (x) sin(2p) vanishes if either p = 0, p = π/2 or W (2x ′ ) = 0. This last case is called "absorbing" for a reason that will be explained later.
Let us compare three trajectories at different energies for U < U c , which are representative of all the cases encountered. For this purpose, two representations are shown: the evolution of x(t), p(t) in figure 2, the phase space in figure 3 (right panel). The three types of trajectories are: • (A) has two p = 0 turning points, and thus its momentum p is bounded.
• (B) is a separatrix trajectory in the classical mechanics sense. The left turning point is at p = 0, and the right turning point at x = 0 is "absorbing", the time taken to reach it (or escape from it) is infinite. In the equation of the motion, at the point W (x) = 0, the effective mass tend to infinity.
• (C) has one turning point type at p = 0 and one at p = π/2. p is growing infinitely large with time, which is not pathological because only p modulo π enters the equations.
The diagram in figure 3 (left panel) is similar in the spirit to the figures of effective potential shown to explain central motions in textbooks. It allows one to understand   In the previous discussion we considered a specific form of the phase space and of D − 2W , D + 2W which are valid for certain values of U only. In the next section a careful investigation of the U dependence will be presented. However, before that, we want to stress that the effective potentials D − 2W , D + 2W can take three different qualitative form, depending on the coupling U. We now discuss the dynamical evolution following a quantum quench and its dependence on the final and initial value of U. At t < 0 the system is in the ground state at the coupling U i . At t = 0 the coupling is switched to U f and the quench dynamics is computed for t > 0. In the following we call a quench "from superfluid" when U i < U c and "from Mott" when U i > U c , and "to superfluid" or "to Mott" when U f > U c and U f < U c . The results of the following sections are summarized by the dynamical phase diagram shown in 6 (left panel).

Mott to superfluid
Starting from the Mott ground state x 0 = 1, the trajectory is stuck at x(t) ∼ 1 even for large times. In order to check this, let us linearize the equation of motion around x = 1. We use thatẋ = ∂H/∂p = 4W x sin(2p) and that sin(2p) can be extracted from x . Thus, at dominant order in ǫ we obtain the equation for ǫ = 1 − x: The trajectory ǫ(t) = 0 is unstable, and since the wave function has a width 1/ √ V , its typical evolution is given by ǫ(t) = 1/ √ V e t/τ . Therefore, in the effective picture, the trajectory is stuck at the Mott ground state on times of the order of log(V ). This result is a peculiarity, actually a pathology, of mean field models. Indeed, in a real finite dimensional system, spatial fluctuations drive the system away from the Mott state in a finite time.

Superfluid to Mott
In a superfluid to Mott quench, the initial condition is given by the ground state packet characterized by {x = U/Uc+1 2 , p = 0}. The trajectory after the quench is of the type (C), see quench Q1 of figure 5. The superfluid order parameter |Ψ 2 0 | (equation (2)) oscillates.

Superfluid to superfluid
Depending on the value of U i and U f , there are three different types of dynamical evolutions in this case, that we call Q2, Q3 and Q4 in figure 5. Q2 is of type (C) and Q4 is of type (A). Q3, reached from a special U d i (dependent on U f ), corresponds to a singular separatrix of type (B) and of energy E = 0. In this case the packet relaxes exponentially to the Mott insulator ground state |x = 1 , and concomitantly the superfluid order |Ψ 0 | 2 vanishes exponentially in time. This may be seen from the linearization (27), which becomesǫ = −ǫ/τ for a trajectory relaxing to x = 1. Approaching the transition oscillations take place on a time scale that diverges as −τ ln(|U d f − U f |). In consequence we find that a dynamical singularity, or transition, takes place at U d f . The values of U d f depends on U i and thus defines a dynamical transition line U d f = (U c + U d i )/2, in the U i , U f plane. In figure 6 (right panel), as an example of singular behavior, we show the time average |Ψ 0 | 2 as a function of U f for quenches starting from the non interacting case U i = 0. It is interesting to compare to |Ψ 0 | 2 its equilibrium counterpart obtained from the microcanonical average corresponding to the same energy. From this we clearly see that the system is not thermalized. Actually, large quenches (Q2) have non vanishing oscillations of the superfluid order contrary to the equilibrium value which is instead zero (because it corresponds to an effective high temperature). Moreover, we find that the dynamical transition takes place at an energy at which the system, if it were relaxed, would be superfluid. Thus, the exponential relaxation of the superfluid order is a purely dynamical effect.

Superfluid to superfluid U f < U d
In this regime, there is a qualitatively new family of quenches Q6 (comparable to Q2). There is also a new special quench Q5 to the modified separatrix state at energy E d , which gives rise to another dynamical transition. However, unlike the other quench cases, these effects are artifacts of the truncation n max b = 2 and disappear for n max b ≥ 3. Therefore, they are not indicated in figure 6 (left panel), where the outcome of all possible sudden quenches is summarized.    (9): Specifically, when n max b = 3, if one chooses x 1 and x 2 as free variables, the Hamiltonian is Notice that because 0 < x i < 1, there are constraints on the possible values of x i , such as x 2 < (n − x 1 )/2. For all n max b , there is a Mott insulator to superfluid quantum phase transition at some coupling U c if the density n is an integer (lower than n max b ). Above the critical coupling U > U c the ground state is a Mott insulator x n = 1, x i =n = 0 (all sites have n bosons), and below U < U c the ground state is superfluid with all x i = 0.

Regularity of trajectories after a quench
In the previous case n max b = 2, the effective dynamics was one-dimensional, and thus integrable. For n max b > 2, the classical dynamics takes place in two or more dimensions, and the trajectories may be either regular or chaotic. In order to characterize them we study their regularity properties. For chaotic trajectories, neighboring trajectories separate exponentially in time in the phase space y = {x i , p i }, like δy(t) ∼ exp(λt)δy(0). The rate of separation λ is the largest Lyapunov exponent [32]. We computed λ for n max b = 3 at density n = 1 for several trajectories, using a simplified version of the traditional Gram-Schmidt orthonormalization of the Lyapunov vectors. Roughly speaking, if we write the Hamiltonian evolution under the formẏ i = f i (y), the deviation satisfiesδ y i = ∂ j f i (y)δy j . This can be integrated numerically and normalized at each step to avoid an overflow (the orthogonalization step is dedicated to find all Lyapunov exponents, whereas here we only need the largest one). The Lyapunov exponents for trajectories with different initial conditions {x i 1 , x i 2 , p 1 = p 2 = 0} are plotted in figure  7 in the superfluid phase U = 2.86 < U c . Trajectories are regular in some regions of the space (periodic or quasi periodic) with λ = 0 within the error bar, whereas some other regions are chaotic with λ > 0.1. The quench from E = 0 is exactly at the dynamical transition, the corresponding trajectory is chaotic. For U > U c , when the ground state is a Mott insulator, all trajectories are regular. We notice that the regularity of a trajectory affects the time of spreading of the packet (determined by the time of separation of two neighboring trajectories). Actually, the time of separation is typically polynomial t ∼ V 1/α for a regular motion but only t ∼ ln V for a chaotic one.

Dynamical transition
In quenches from the superfluid phase U i < U c , like in the previous case n max b = 2, a dynamical transition occurs at the special coupling U d f where the final energy equals the energy E = 0 of the unstable Mott trajectory x n = 1, x i =n = 0. Direct evidence of this transition is given by the singularity of the time averaged superfluid order |Ψ 0 | 2 as a function of U f for a given U i . In figure 8b, we compare this singularity in |Ψ 0 | 2 for n max can safely assume that this phase diagram is quantitatively representative of the phase diagram without truncation.
Even though the existence of the dynamical transition for any n max b is beyond doubt, because the singularity is numerically manifest, more precise results on this transition, even for n max b = 3, are hard to provide. Some features of the case n max b = 2 persist ; for example, at the dynamical transition, the momentum 2p 1 − p 2 becomes unbounded, see Fig. 8a. The fractions x i (t) are also oscillating, but without definite period. They are either quasi-periodic for low λ regions, or chaotic. Approaching the dynamical transition, the typical time of return to x 1 ∼ 1 is increasing, possibly logarithmically divergent in U f − U d f as suggested by the numerical divergence in figure  8b. Unfortunately, the analysis of the singularity for n max b ≥ 3 turns out to be out of reach. The numerical integration of classical equations of motions is exponentially sensible to numerical errors and thus not reliable. One could imagine that the trajectory at E = 0 is an unstable manifold, which would support the existence of a singularity for trajectory arbitrarily close to it. However, it is hard to decide whether the trajectory at the dynamical transition goes arbitrarily close to the point x 1 = 1 (ground state of the Mott insulator), at which the trajectory is exponentially slowed down. A possible scenario, in which the surface E = 0 is ergodic, does not seem to be validated by numerical integration of trajectories.

Super-radiance quantum phase transition
We now consider the generalized Jaynes-Cummings model (4) and derive its quench dynamics. For simplicity and without loss of generality we focus onQ = N/2. Thanks to the conservation ofQ =Ŝ z +b †b , the spin degrees of freedom can be accounted for in terms of the bosonic ones. The states can be parametrized by the density of bosons n only, in the Fock basis |n bosons ⊗ |S; m spin : The action on |n of some useful operators can be readily computed at the dominant order in NŜ z |n From these equations, the transition elements D n = 1 N n|Ĥ|n and W n = − 1 N n ± 1/N|Ĥ|n are easy to compute, and yield an effective Hamiltonian on n and φ its conjugate momentum H[n, φ] = D n − 2W n cos(φ) which is the semi-classical approximation of the original Hamiltonian obtained in [25], up to a shift in the phase φ → φ + π. In [25], the authors studied sweeps from the empty state n = 0, and the link to the Landau-Zener problem. Here, we focus on sudden quenches and the related dynamical transition. The ground state of this effective system is at φ GS = π and The super-radiance quantum phase transition is the result of a competition between the potential energy term and the kinetic energy term, which can be understood if the quantum Hamiltonian is written in the |n basis: Qualitatively, the potential term −2λb †b encourages filling for λ ≥ 0 and discourages it for λ ≤ 0. The "kinetic" contribution ∼ g(b † +b) can lower the energy. As a matter of fact, the sign of g is almost irrelevant. Changing the sign of g leads to W n → −W n . However, the new Hamiltonian can be mapped into the old one by the translation p → p + π. In consequence, if the ground state of a given Hamiltonian is at p = 0, Around n = 1 − ǫ, the time is finite T ∼ 0 dǫ 1/ √ ǫ, whereas around n = ǫ, the time T ∼ 0 dǫ 1/ √ ǫ 2 diverges. In other words, the trajectories touching the turning point n = 1 have a finite period, whereas the trajectories touching n = 0 have infinite period. The situation is anologous to the one already studied for the BHM: the absorbing state n c = 0 plays the same role of the Mott state, and gives rise to the dynamical transition.
To identify the critical value λ d f at which the transition takes place, one can compute the energy after a quench λ i → λ f and compare it to the energy of the system at rest in n = 0. One can check that a trajectory n(t) starting with the same energy, actually energy zero, is indeed driven to the absorbing point n c = 0. This dynamical transition occurs on the line λ d f = − 1 − n GS (λ i ). A linearization similar to (27) shows that the relaxation is exponential at the transition with a characteristic time τ −1 = 2 √ 1 − λ 2 . Around the transition, the period of a trajectory diverges as τ ln(|λ f − λ d f |). The dynamical phase diagram is shown in figure 9b. For λ i > λ c , the initial state is in the broken symmetry phase n(t = 0) = 0. The situation is very similar to the BHM one: the filling number n(t) oscillates in time after the quench, and there is one region with bounded effective momentum φ, and two regions with unbounded φ. Notice that there are two separatrices between the three regions, and that one is not singular, and does not give rise to a dynamical transition, whereas the other is. The singular one is the one for which the absorbing state n = 0 is met, corresponding to the restored symmetry state. For quenches with λ i < λ c , the initial state is empty n(t = 0) = 0, and either drifts away from zero on large times (of the order of log N) if λ f < 1, and otherwise, the state is definitely stuck at n = 0. The bounded or unbounded nature of the final state is also indicated.
Finally, we remark that eigenstates of the completely connected model can be written within a WKB expansion, this is done for completeness in Appendix B.

Quenches in the Ising model in a transverse field
In our last example, let us consider the dynamics due to quantum quenches in the transverse field Ising model. To derive the classical effective dynamics, one can use the fact that the large spin limit is also the classical limit, see e.g. [21]. To do so, one should make the substitution Ŝ → S = S{sin θ cos φ, sin θ sin φ, cos θ}. The kinetic term (needed to describe classical dynamics) in the classical Hamiltonian of the large spin limit can be derived using path integral for spins [33]. Here we instead use our generic method based on site permutation symmetry. We denote spin symmetric states of total momentum S = N/2 with the notation |s . They are such thatŜ z |s = Ns|s with s ∈ [−1/2, 1/2]. At dominant order in N,Ŝ x |s = N 1 4 − s 2 |s + 1 N + |s − 1 N /2. One can proceed as usual to get the effective Hamiltonian: As expected, the Hamiltonian is very reminiscent of the classical Hamiltonian for a single rotor once the change of variables s → cos θ is made. Indeed, the equations of motion obtained here are the same as is [21]. In this paper, the authors studied the quench from the paramagnetic phase to the ferromagnetic phase and AC dynamics.
Proceeding as for the BHM, we analyze the phase space of symmetric trajectories. This is different in the ferromagnetic and paramagnetic phases, see figure 10. In the latter, obtained for Γ < Γ c , there is a separatrix. For sudden quenches from the ferromagnetic phase, where s = 0, to final values of the transverse field Γ f < Γ c , there is a dynamical transition at Γ d f = 1 4 + Γ 2 i . At the transition, the trajectory is a separatrix, and s decreases exponentially in time to s = 0, which is the quantum paramagnetic ground state. The symmetry is dynamically restored at this point but also for larger quenches: contrary to the two previous examples, for quenches beyond the dynamical transition Γ f > Γ d f , the trajectories are symmetric in s → −s, and the magnetization is oscillating around zero. Averaging over time, the order parameter is zero Ŝz = 0.

Discussion
Before concluding, two points need to be addressed further: (1) the possible link between dynamical transitions and equilibrium quantum phase transitions and (2) the pro and cons of our approach concerning the physics of finite dimensional systems.
• In three examples, we found that dynamical transitions occur in systems characterized by a quantum equilibrium phase transition. Therefore, despite the fact that the dynamical transition is not directly related to the equilibrium one, as we stressed before, it is natural to wonder whether instead its existence is related to it. In support of this, if one studies the completely connected Bose-Hubbard model in a parameter regime (non unitary fillings) where there is no equilibrium quantum phase transition then one finds that there is no dynamical transition either. Actually, we found the very same phenomenon in all models we studied, and this was also noticed for the Fermionic Hubbard model [11]. We certainly cannot address this issue in general, also because it is still far from clear what "the dynamical transition" is beyond mean-field theory. We shall instead endeavor to provide a partial answer at the mean-field level. Generically, we observe that in all models we studied, there is a regime in which the symmetric state (the Mott insulator, the paramagnet, etc.) is not the ground state but it plays the role of an excited metastable state that can trap the system for an infinite time. This gives rise to the dynamical transition. If the effective motion is one dimensional, this state is absorbing and belongs to a separatrix. If the effective motion has two or more degrees of freedom, this state seems to belong to an unstable manifold. Thus, it appears that the singularity at the dynamical transition is a result of this singular trajectory, which in turn is a consequence of the existence of a quantum phase transition. A physical argument to support this idea is the following. For a quench from the unbroken symmetry phase to the broken symmetry phase in completely connected model, the dynamical evolution is expected to be very slow. The reason is that the symmetry must break up and the order must develop starting from this point. This, however, is a local phenomenon, for example due to domain growth in finite dimensions, and therefore generically absent in completely connected models, except on times diverging with the system size: indeed, in our examples, we find that the typical time scale to depart from the initial value x(t = 0) is of the order of log V . The point is that this slow trajectory is precisely the time reversal of the singular trajectory that is responsible for the dynamical transition. To summarize, the existence of a quantum phase transition implies the existence of slow trajectories. They correspond to the escape from the symmetry unbroken state. Their time reversed counterparts are the trajectories responsible for the dynamical transition. Thus, the existence of a phase transition indeed appears to be intertwined with the existence of a dynamical transition within meanfield.
• Let us now discuss the range of applicability of our results for finite-dimensional systems. First, our analysis can be put on the same footing as the Gutzwiller Ansatz and functional integral saddle point approximations, since the corresponding dynamical equations of motion coincide (see Appendix A for the former equivalence and also [22]). Clearly, our mean-field approximation misses several essential physical effects. Relaxation and thermalization, for which spatial and temporal fluctuations must be taken into account, are not present. Furthermore, the dynamics of a quench from the unbroken symmetry phase to the broken symmetry phase is not properly described for the reasons mentioned above. Inhomogeneities, topological defects, domain growth are out of reach. On the contrary, our mean field approximation is expected to capture well the evolution of local quantities at least on short times, as long as they are homogeneous across the system. Contrary to many other approaches, it is not perturbative in parameters (U, J, Γ . . . ) of the Hamiltonian. Unlike pertubative approaches, it can describe the existence of dynamical transitions, and phenomena related to global observables.

Summary and outlook
The general purpose of this article was twofold. First, we described a method to study the quantum quench dynamics of generic completely connected models, providing a mapping to effective classical dynamics. There are two reasons for considering completely connected models. In some cases, such as in the generalized Jaynes-Cummings model, and the Dicke model, they provide the correct physical description. In others, such as the transverse field Ising model or the Bose-Hubbard model, they lead to an approximative treatment of finite dimensional systems. In the latter case, the range of validity of the approximation is possibly limited to short times only. Our second aim was to study and to reveal the existence of out of equilibrium dynamical transitions induced by quantum quenches. In agreement with other works [11,13] we showed that within the mean-field approximation dynamical transitions occur in quenches from the broken symmetry phase to other regions of the broken symmetry phase and that this is a quite generic phenomenon for systems that display a quantum phase transition at equilibrium. Clearly, a main and pressing question is to understand what this transition is really, beyond mean-field theory. Does it become a cross-over in finite dimensions? If it remains a bona fide transition, what are its critical properties? In the d = 1 BHM, a nonmonotonic behavior of the propagation velocity as a function of U f has been found in [34] by t-DMRG. This could well be a signature of a cross-over related to the dynamical transition found in mean-field. An exact diagonalization study of hard core bosons in d = 2 suggests the existence of a cross-over or a singularity of the revival time after a quantum quench [23]. From the analytic point of view going beyond mean-field is a difficult task. An argument for the occurrence of dynamical transitions beyond mean field is given in [13], within an imaginary time path integral formalism. Clearly, it is crucial to take into account fluctuations not captured by mean-field theory. This was started to be done in [35]. This work indeed suggests that these fluctuations could alter substantially the mean-field results. Promising ways to capture fluctuations are the projector operator formalism [36], 1/z expansions [37] (z is the connectivity of the lattice). Field-theoretic methods are also being developed, as in [19] for the Bose-Hubbard model in the grand canonical ensemble; the use of two particle irreducible actions also seems promising [38].
As a consequence, the Gutzwiller Ansatz wave function is indeed equivalent to an effective classical Hamiltonian evolution of the variables x l and p l . The effective Hamiltonian H[x l , p l ] = ψ|Ĥ|ψ [x l , p l ] must be computed for the model at hand, and is indeed the same as the one discussed in section 4. For instance, for the Bose-Hubbard model with truncation n max b = 2, one recovers (15).

Appendix B. Eigenstates in completely connected models
Using a Wentzel-Kramers-Brillouin (WKB) approximation [39], eigenstates of the Schrödinger equation (13) can be found in the limit of large V . To do so, we write the eigenstates φ E (x) within the WKB approximation. The derivation is standard, one looks for stationary solutions of the form ψ(x) = A(x)e iV S(x) of the equation (13), at dominant order in V . One finds where C is a normalization constant, E is the energy, p(x) satisfies the implicit equation E = D x − 2W x cos(2p(x)), and x 0 is the left turning point of the classical trajectory of energy E.
Using this expression, we can draw a parallel between the eigenstate φ E (x) and the effective classical trajectories of energy E. For an observable f (x), the average reads In the second line, we refer to the effective classical motion dx dt = ∂H ∂p , and the average f (x) is the average over the effective classical trajectory of f (x) over one peridod. In words, we found that in the semi classical regime, the average over a quantum stationary state of energy E is given by the average over one period of the classical trajectory of same energy E. Although this result is not often mentioned, this is a direct consequence of the WKB expression of the wave function.