Changing the order of a dynamical phase transition through fluctuations in a quantum p-spin model

We study the non-equilibrium phase diagram of a fully-connected Ising $p$-spin model, for generic $p>2$, and investigate its robustness with respect to the inclusion of spin-wave fluctuations, resulting from a ferromagnetic, short-range spin interaction. In particular, we investigate the dynamics of the mean-field model after a quantum quench: we observe a new dynamical phase transition which is either first or second order depending on the even or odd parity of $p$, in stark contrast with its thermal counterpart which is first order for all $p$. The dynamical phase diagram is qualitatively modified by the fluctuations introduced by a short-range interaction which drive the system always towards various paramagnetic phases determined by the strength of time dependent fluctuations of the magnetization.


I. INTRODUCTION
Equilibrium phase transitions, either at zero or finite-temperature, are known to leave a substantial imprint in the non-equilibrium dynamics of a quantum many-body system [1].For example, even when a stationary state attained after a quantum quench does not reveal signatures of order as in low-dimensional systems [2][3][4], a linear ramp through a second order quantum critical point leaves universal signatures in the scaling of the number of excitations with the ramp speed [5][6][7][8][9][10][11], as confirmed extensively in a number of experiments [12][13][14][15][16][17][18].Analogous signatures are left when a first-order quantum phase transition is crossed [19][20][21] through the nucleation of resonant bubbles of the new phase close to the critical point [22][23][24] which leads to a modified Kibble-Zurek-like power-law scaling [25].
Among the signatures of criticality observed out-of-equilibrium, Dynamical Phase Transitions (DPT) occupy a special place.A dynamical quantum criticality can be observed as singular temporal behavior of the Loschmidt echo (LE) most notably after a quench across a quantum phase transition [26][27][28], even in situations where long-range order cannot be sustained in stationary states.In systems with long-range interactions, on the contrary, intertwined with the singular behavior of the LE [29][30][31][32][33][34], a standard Landau-type critical behavior based on the dependence of a time-averaged order parameter with respect to the quench parameters can be observed [35,36].Peculiar to the second-order dynamical transitions arising in these models is the fact that they are associated to critical trajectories with a divergent time scale in the dynamics separating revivals with a finite order parameter [29].In the presence of fluctuations critical trajectories become unstable and second order dynamical critical points widen up into chaotic dynamical phases [37][38][39][40].
While a great deal is known about second-order dynamical phase transitions, the dynamics of systems displaying equilibrium first order transitions is much less explored.The notion of dynamical criticality associated to the LE has been extended to include first order behavior [41] while first-order and dissipative phase transitions in infinite range p-spin systems coupled to an external bath have been studied in Ref. [42].However, dynamical transitions occurring in systems displaying first order equilibrium transitions are much less studied.
In this work we address this issue by studying dynamical phase transitions and their stability against fluctuations in a system displaying a first order equilibrium transition: a spin system with infinite range p-spin interactions in a transverse field.We show that, already at mean-field level where the dynamics is effectively classical, the system undergoes a DPT after a quench of the transverse field g, whose order depends non-trivially on p, despite its equilibrium counterpart being always of first order.In particular, we show that the order of the dynamical transition can be inferred entirely from the profile of the underlying energy landscape.We then perturb the model by a short-range two-body interaction tuning the strength of spin fluctuations [37,38,43].While for p = 2 a chaotic dynamical region opens up near mean-field criticality [37,38], we show that for p > 2 dynamical chaos is almost entirely replaced by a new prethermal regime, which we define as "dynamical paramagnetic phase", which appears for sufficiently large short-range coupling.This is due to the emission of energy in the form of spin-waves, which predominantly drive the system into a paramagnetic minimum even in the presence of other minima in the energy landscape.

II. MEAN-FIELD DYNAMICS
In this section, we study the dynamics of N spins 1/2 subject to all-to-all p-body and a global transverse field g.The corresponding Hamiltonian is given by: Here, the operators σα i denote the Pauli matrices at site i.The fully connected p-spin model in Eq. ( 1), which was originally introduced in the context of spin glasses [44,45], plays a central role for studies on quantum annealing [46,47].Its zero-temperature equilibrium phase diagram can be derived analytically [46][47][48][49], utilizing a Suzuki-Trotter decomposition in the thermodynamic limit [50].The phase diagram displays a quantum phase transition driven by g and detected by the magnetization S x = i ⟨σ x i ⟩ /N along the x-axis.The transition, located at some g = g c , separates a ferromagnetic state, defined by a non-vanishing S x , from a paramagnetic where S x = 0.The transition is continuous for p = 2, where S x vanishes with a square-root singularity [48], while it is of the first order for p > 2, where S x displays a discontinuity.In this section we will address the dynamics and dynamical phase transitions of this p-spin model.6), projected on the plane S y = 0.The projection is obtained from the black cut on the Northern hemisphere of the corresponding Bloch spheres, where S z = 1 − (S x ) 2 .Each panel illustrates the shape of the energy landscape corresponding to a different set of values of p: p = 2 (left), p = 2n + 1 (center) and p = 2n + 2 (right), respectively, with n ≥ 1 integer.For practical purposes, we fixed n = 1.Each profile in each panel corresponds to a different value of g: g = g dyn − 0.1 (a), g = g dyn + 0.1 (b) and g = g dyn (c).The value of g dyn is obtained from Eqs. ( 8) and (10), as discussed in the main text.For each energy profile, the corresponding trajectory evolving from ⃗ S(0) = x (black dot on the profile) is represented by a dashed horizontal line.In particular, the gray dashed horizontal line corresponds to a separatrix (white trajectory on the Bloch sphere).See also the orbits in Fig. 2 (a) and (d) for a comparison.

A. Semi-classical theory for the post-quench dynamics
Our goal is to study the dynamics of the average magnetization, ⃗ S(t) = ⟨ j ⃗ σ j (t)⟩ /N , after a quench in the transverse field g in Eq. (1), where the field g abruptly changes from a value g 0 = 0 to a value g > 0. Specifically, we prepare the system in the fully polarized state, corresponding to a ground state Ĥ0 at g 0 = 0.The dynamics of ⃗ S(t) is obtained by averaging over |ψ 0 ⟩ the Heisenberg equations generated by the Hamiltonian Ĥ0 .However, due to the all-to-all interacting nature of Ĥ0 , this dynamics becomes effectively classical in the thermodynamic limit N → ∞.This classical behavior is derived from the general theory outlined in Ref. [35], which we briefly review in the following.
To begin with, we observe that the Hamiltonian in Eq. ( 1), for the post-quench transverse field g = g, can be rewritten in terms of the collective spin operators, Ŝα = j σα j /N , as follows: The collective spin components satisfy the commutation relations where ϵ αβγ is the Levi-Civita symbol.The commutators, controlled by an effective Planck constant 1/N , vanish in the thermodynamic limit.As a consequence, the dynamics of the average magnetization becomes effectively classical for N → ∞ and is governed by the Hamilton equation (see also [42,51]): The right-hand side of Eq. ( 5) is derived by substituting the rescaled average commutators, N ⟨[ Ĥ0 , Ŝα ]⟩ /i, with the corresponding Poisson brackets, {H cl ( ⃗ S), S α }.
Here, the effective classical potential H cl ( ⃗ S, g) is obtained as the thermodynamic limit of ⟨ Ĥ0 ⟩ /N and is given by The initial condition for Eq. ( 5) is ⃗ S(0) = x, corresponding to the fully polarized state in Eq. ( 2).Here and throughout this paper, we use the notation {x, y, z} to denote the unit vectors along the corresponding axes.We also observe that, as the modulus of the magnetization is a constant of motion, the classical dynamics from Eq. ( 5) takes place on the Bloch sphere | ⃗ S| 2 = 1.Thus, the system under consideration is always in a non-equilibrium state, as all the microscopic spins perform a coherent, undamped precession.
The dynamics from Eq. ( 5) is strongly influenced by the shape of the effective Hamiltonian H cl ( ⃗ S, g).Depending on the value of p and for sufficiently small values of g, the profile H cl ( ⃗ S, g) on the Bloch sphere exhibits various topologies, characterized by the number and positions of its maxima and minima.To determine their locations, we parameterize the magnetization with the spherical angles (θ, ϕ) ∈ [0, π/2]×[0, 2π], as ⃗ S = sin θ cos ϕ, sin θ sin ϕ, cos θ .The stationary points of H cl ( ⃗ S, g) are then defined by the equations: In the following, consider only the stationary points falling in the in the Northern hemisphere (S z > 0) of the Bloch sphere, where the dynamics is confined for g > 0. One possible solution of the system of Eqs. ( 7) is given by θ = 0 (the North Pole ⃗ S = z of the Bloch sphere), being a maximum for p = 2 and a minimum for p > 2. All the other solutions are obtained by solving the system: From the first equation, we get that stationary points lie in the plane S y = 0 for every p, while the number and the precise location of the solutions depend on the value of p. Specifically, for p = 2, we observe two symmetric minima separated by a maximum at the North pole, ⃗ S = z.For p ≥ 3 odd, the same topology persists, but the profile becomes asymmetric with respect to the North pole.Conversely, for p ≥ 4 even, the potential features three minima: one at the North pole and two symmetric minima with respect to it.To summarize, the projection of H cl ( ⃗ S, g) in the plane S y = 0 exhibits three distinct shapes, corresponding to different values of p: -A symmetric double-well for p = 2.
-An asymmetric double-well for p ≥ 3 (odd), with one paramagnetic and one ferromagnetic minimum.
-A symmetric triple-well for p ≥ 4 (even), with one paramagnetic minimum and two opposite ferromagnetic minima.
It is worth noting that these three profiles undergo a qualitative transformation beyond the spinodal point [47], defined as More precisely, for each value of p the profile of H cl ( ⃗ S, g) becomes a single well, centered around its only minimum at ⃗ S = z, when g > g sp .For these values of g, the second of Eqs.(8) has no solution.

B. The dynamical transition
The presence of multiple local minima in the profile has a strong impact on the dynamics of ⃗ S(t), which evolves according to Eq. ( 5) from the initial condition ⃗ S(0) = x, located in the rightmost ferromagnetic well.We denote the position of the local maximum that separates this well from the rest of the landscape as ⃗ S m (g).The dynamics exhibits qualitatively different orbits, depending on the value of the post-quench transverse field g [52], as qualitatively depicted in Fig. 1: 1. Below a certain threshold, g < g dyn (trajectory (a)), the dynamics starts from a ferromagnetic well, and the initial energy is insufficient to surmount the energy barrier in correspondence of nearest maximum in ⃗ S = ⃗ S m (g).Consequently, the magnetization oscillates in the ferromagnetic well, with S x (t) > S x m (g) at every time t.
2. For g > g dyn (trajectory (b)), the post-quench energy is larger then the energy barrier and the corresponding orbit encompasses all the minima in the landscape.
3. At precisely g = g dyn (trajectory (c)), the system has sufficient energy to reach the top of the barrier in ⃗ S m (g), but it is unable to surpass it.Here the period of the oscillations, T cl (g), diverges and the magnetization approaches ⃗ S m (g) infinitely slowly.The resulting orbit forms a separatrix.
The three cases listed above correspond to three different possible topologies for the underlying orbits at large times.We classify each distinct topology as a dynamical phase [53].Consequently, the singular dynamics at g = g dyn leads to a dynamical phase transition (DPT).The corresponding dynamical critical point g dyn is obtained by equating the energy of the initial configuration, ⃗ S(0) = x, with that of the local maximum in ⃗ S c = ⃗ S m (g dyn ).In terms of the variables θ and ϕ, the equation reads: Then, the simultaneous solutions of Eqs. ( 8) and (10) determines both g dyn and the spherical coordinates of the maximum in ⃗ S c .From a practical point of view, the DPT can be studied also in terms of a "dynamical order parameter", such as the time-averaged longitudinal magnetization As discussed also in Appendix A, the two indicators are equivalent to each other: as g approaches g dyn , the change of topology in the orbits is signalled by a divergence of their period T cl (g), in turn creating a non-analiticity in the function S x (g).In particular, both T cl (g) and S x (g) display a log-singularity while approaching the dynamical transition, a feature already known for other mean-field models driven away from equilibrium [35,54].It is important to note that such a DPT can be observed in the N → ∞ limit only: for a finite size, quantum fluctuations are restored and induce dephasing between local spins, leading to the relaxation of time-averaged observables to their thermodynamic expectation values [38].In the following, we show that the nature of the mean-field DPT is determined by the topology of the landscape, which depends solely on the parity of p (if p > 2).As a consequence, we will focus on studying Eq. ( 5) for p = 3 and p = 4 only.These cases, alongside with p = 2, are paradigmatic and encapsulate the three possible landscape shapes discussed in the previous section, respectively.
First, it is crucial to note that the order of the DPT is not necessarily identical to that of the thermal phase transition for the p-spin model at the same p value.The case of p = 2 is special, exhibiting both a second-order thermal phase (Color online) From left to right: orbits of the magnetization on the Bloch sphere ((a) and (d)), time-averaged magnetization S x calculated from Eq. ( 11) as a function of g ((b) and (e)), time evolution of S x (t) near the critical point ((c) and (f)).The plots in the upper row corresponds to the case of p = 3, the ones below to p = 4. (Left) Orbits on the Bloch sphere, resulting from the integration Eq. ( 5), from the initial condition ⃗ S = x.Each trajectory corresponds to one of 20 values of g,uniformly sampled in the interval [0.2, 3.6], where g dyn falls for both p = 3 and p = 4.For both values of p, the system evolves along trajectories that are either confined in a ferromagnetic well (green trajectories), for g < g dyn , or encompass all the local minima (orange trajectory), for g > g dyn .The two regions are divided by a separatrix (white), corresponding to g = g dyn .g dyn is calculated from Eqs. ( 8) and (10),as discussed in the main text.(Center) Plot of the dynamical order parameter S x , from Eq. ( 11), as a function of g.The red dot correspond to the dynamical critical point (g dyn , S x c ), obtained by simultaneously solving simultaneously solving Eqs. ( 8) and (10).At the dynamical critical point, S x is continuous for p = 3 only.(Right) Time evolution of the longitudinal magnetization S x (t), for a transverse field g either equal to g dyn − 10 −7 (green plots) or to g dyn + 10 −7 (orange plots).For each value of p, the dashed lines correspond to the local maxima of the landscape when g = g dyn .For each value of p, g dyn and S x c are obtained as in the central panels.
transition [55] and a second-order dynamical phase transition [29,35,48].For p = 2, the dynamical transition occurs at g = λ and is identified through S x , which is positive when g < λ (indicative of a dynamics confined within a ferromagnetic well) and vanishing for g > λ (due to oscillations between the two symmetric wells).We classify these two behaviours as dynamical ferromagnetic and dynamical paramagnetic phases, respectively.The continuity of this DPT is not just a consequence of the symmetry of H cl ( ⃗ S, g) under the reflection S x → −S x .Instead, it emerges as a general property linked to the topology of the phase space, which is characterized by only one maximum.To understand why, we notice that, as g approaches g dyn either from above or below, the energy of the orbit gets close to the one of the separatrix.The separatrix has a crossing at the only local maximum of H cl ( ⃗ S, g dyn ).Thus, every dynamics starting on the separatrix approaches asymptotically this maximum in ⃗ S c ( ⃗ The orbits asymptotically close to the separatrix, approached as g → g ± dyn , exhibit a plateau at the only local maximum in ⃗ S c , whose length diverges with the period T cl (g).As a consequence, we have that lim and the dynamical transition is continuous.
For p > 2, the results illustrated Fig. 2 consistently demonstrate the occurrence of a DPT at some g = g dyn , whose order depends on the parity of p, unlike its static counterpart.For p = 3 the topology of the phase space is the one of a double well, akin to the case of p = 2 but without Z 2 symmetry.Consequently, the transition is continuous, following the same argument made for p = 2 and as illustrated in Fig. 2-(b).Specifically, a dynamical ferromagnetic phase re-emerges for g < g dyn , as all the orbits are confined in a ferromagnetic well (green trajectories, Fig. 2-(a)).Conversely, for g > g dyn , the oscillations of the orbits between the two asymmetric wells (orange trajectories, Fig. 2-(a)) lead to a dynamical bistable phase.In this phase, S x is a weighted average of the two minima, that varies continuously with g and exhibits a cusp at g = g dyn .The connection between the continuity of S x to the presence of a single local maximum is quantitatively validated by the results presented in Fig. 2-(c).This figure illustrates the plateaus developed by the longitudinal magnetization near g dyn , both for g ≲ g dyn (green plot) and g ≳ g dyn (orange plot).
For p = 4 the DPT becomes discontinuous, as shown in Fig. 2-(e).The discontinuity results from a drastic change in the topology of the energy landscape, which now exhibits one paramagnetic and two opposite ferromagnetic minima.Specifically, for g < g dyn , the system is still in a dynamical ferromagnetic phase, where the orbits (green plots in Fig. 2-(d)) are confined in the rightmost well and S x > 0.Then, approaching g dyn from below, S x still tends to the value determined by the position S x c of the nearest, non-vanishing local maximum.However, upon moving even slightly above g dyn , the orbits (orange plots, Fig. 2-(d)) become symmetric with respect to the y − z plane, resulting in S x = 0 and marking a discontinuous dynamical transition.From a broader perspective, the discontinuity arises from the existence of two maxima of H cl ( ⃗ S, g dyn ), having the same energy.These maxima are at located at S x = ±S x c , the two points where the separatrix (illustrated as a white trajectory in Fig. 2-(d)) has crossings.As a consequence, S x (t) displays two distinct plateaus for g ≳ g dyn , around S x = ±S x c (see green and orange plots in Fig. 2 (f)): when g approaches g dyn from above, S x converges to the vanishing average of the two plateaus, creating a discontinuity at g dyn .

III. NON-EQUILIBRIUM SPIN-WAVE DYNAMICS A. Non-equilibrium spin-wave theory
The dynamical transitions investigated in Sec.II-(b) originate from the coherent dynamics of the local spins in the N → ∞ limit.In this limit, the spins are effectively decoupled from each other due to the all-to-all interaction among them.However, this coherence is expected to be unstable with respect to the inclusion of short-range interactions, which introduce fluctuations in the dynamics and drive the system towards eventual thermalization.However, if the short-range coupling is small, the amplitude of these fluctuations is expected to be small for a parametrically long time, leaving possible instances of dynamical phases in the prethermal stage of the dynamics [56][57][58].
This scenario was investigated in Ref. [37], using the non-equilibrium spin-wave theory (NEQSWT).For p = 2, it was shown that the dynamical critical point retrieved at the mean-field level is melted into an entire chaotic region of the phase diagram when short-range fluctuations are included.In this work, we extend this analysis to the more general case of p > 2, employing the NEQSWT to study the same quench dynamics discussed in Sec.II, under the influence of an extra, short-range term in the Hamiltonian.For the purposes of our study, we assume that our system is on a one-dimensional lattice with periodic boundary conditions, so that the short-range interaction term is expressed as [37]: To provide a clearer interpretation of Eq. ( 14), we introduce Fourier modes σα k = N j=1 e −ikj σα j , with k = 2πn/N and n = 0, . . ., N − 1.Here, N denotes the system size.Then, the full Hamiltonian is given by Ĥ depends on the zero-momentum mode components σα 0 , related to the magnetization through S α (t) = ⟨σ α 0 ⟩ /N , solely via Ĥ0 .On the other hand, the perturbation Û from Eq. ( 14) comprises only k ̸ = 0 contributions, which are In particular, the yellow and the orange regions correspond to dynamical ferromagnetic and paramagnetic phases, respectively.(Right) The corresponding orbits are either confined in the ferromagnetic well (plot T1) or encircle both the two minima (plot T2).In the blue region, instead, is the dynamical paramagnetic region where the magnetization revolves around the minimum on the North pole (plot T3); here S x is closer to 0 than in any other phase, though non-vanishing, and increases discontinuously when moving across the border with the orange zone.The crossover region is instead a narrow chaotic phase, reminiscent of the one found in [37], where collective spin can localize at large times in either the paramagnetic or the ferromagnetic minima (plot T4).
the ones expected to induce dynamical fluctuations in the magnetization.The expression in Eq. ( 14) is then the simplest perturbation that breaks the permutation symmetry in Eq. (1).Nevertheless, as we shall discuss later, our findings are expected to be independent of the range of interaction and of the dimensionality of the lattice.Hereafter, we outline the technical steps needed for implementing the non-equilibrium spin-wave theory (NEQSWT), for the Hamiltonian presented in Eq. (15).
The fundamental hypothesis underlying NEQSWT is that, at least as long as J is sufficiently small, the net effect of the term from Eq. ( 14) is to give rise to small spin-wave excitations on top of the classical magnetization ⃗ S(t).In particular, the magnetization length is still close to its maximal value, | ⃗ S(t)| ≃ 1, when the short-range fluctuations are sufficiently small.This allows the dynamics to be still effectively described by trajectories near the Bloch sphere, perturbed by fluctuations induced by the finite k degrees of freedom.The NEQSWT is then implemented by studying the dynamics in a time-dependent, rotating reference frame R, identified by the time-dependent Cartesian vector basis {X(t), Y(t), Z(t)}.The frame R is constructed such that the magnetization is aligned with the Z-axis at any time t.The base vectors of R can be identified by their components in the original frame {x, y, z}, given by: We implement this change of frame in the Hilbert space, through the time-dependent rotation operator V (t) = exp −iϕ(t) j σz j /2 exp −iθ(t) j σy j /2 , so that the spin operators transform accordingly:  25) and ( 29) for the same default parameters of N ,λ and the maximum time T listed in the caption of Fig. 3.For both the plots, each color corresponds to a value of the observable we plot, as specified in the legends reported on the right of each diagram.(Left) From the plot of S x , we identify two main regions: one corresponds to a dynamical ferromagnetic phase(yellow), where S x > 0 (plot T1), while the other is a paramagnetic phase (orange) where S x = 0.In between, some chaotic spots (blue and yellow spots) are found occasionally, where the magnetization eventually falls in one of the other two symmetric, ferromagnetic minima (plot T4).(Right) Looking at the average fluctuations (δS x ) 2 , we clearly see that the in dynamical paramagnetic phase, coinciding with the non-blue region of the phase diagram, can be split in a sub-phase 1 (yellow), where the dynamical trajectories either surround symmetrically all the three minima of the landscape (plot T2) and a sub-phase 2 (blue), where the magnetization localizes (predominantly) in a paramagnetic well (plot T3).The border between the two identifies the transition line J = J dyn (g).
In the new frame R, the Heisenberg equation of motion for the components σα j , for α = X, Y, Z, can be written as Here, the modified Hamiltonian includes an additional term i V (t)∂ t V † (t) = −⃗ ω(t) • j ⃗ σ j (t)/2, with ⃗ ω(t) = − sin θ(t) φ(t), − θ(t), cos θ(t) φ(t) , plays a role analogous to the one of apparent forces in classical mechanics.The new Hamiltonian H(t), describing the dynamics in the rotating frame, can be explictily written as: Here and in the following, we will often omit the time dependence of the operators and of the basis vectors to keep the notation compact.
H is the Hamiltonian describing the Heisenberg dynamics in the rotating frame R, where the magnetization ⃗ S(t) is along Z(t).The time evolution of the angles θ(t) and ϕ(t) is determined by self-consistently imposing that the transverse components of ⃗ S(t) vanish in the new frame, that is: Solving Eqs. ( 20) is in general a formidable task.However, as long as the fluctuations transverse to the classical magnetization are small, we reasonably assume that the dynamics is dominated by terms at the lowest nontrivial order in the Holstein-Primakoff (HP) expansion [59]: where s = 1/2 for the current case.As in the case of Ref. [37], we retain perturbative terms from the HP expansion which are quadratic in the spin-wave modes qk = j e −ikj qj / √ N and pk = j e −ikj pj / √ N , i.e. the Fourier transforms of qi and pi .This is equivalent to keep the following terms in H: In Eq. ( 22), we also defined the spin-wave number operator nSW = i (q 2 i + p2 i )/2.Specifically, the terms Û (0) 2 and Û2 represent the quadratic terms from the HP expansion of the mean-field term Ĥ0 and the short-range perturbation Û , respectively.Within this quadratic approximation, imposing the Eqs.( 20) is equivalent to requiring that the average of each coefficient, which appears either in front of q0 or p0 in Eq. ( 22), vanishes self-consistently.After some algebra, this requirement leads to the following equations of motion: s φ =pλ(sin θ) p−2 (cos ϕ) p cos θ 1 − (p − 1)ϵ(t) − g − 2Jδ qq (t) cos θ cos 2 ϕ + 2Jδ qp (t) sin ϕ cos ϕ s θ =pλ(sin θ cos ϕ) p−1 sin ϕ{1 − (p − 1)ϵ(t)} − 2Jδ pp (t) sin θ sin ϕ cos ϕ + 2Jδ qp (t) sin θ cos θ cos 2 ϕ .(25) In Eqs.(25), we defined the "quantum feedback" variables for α, β ∈ {p, q}.These variables couple the classical spin to the corresponding spin-wave correlation functions, defined by defined for each value of k = 2πn/N , where n = 1, . . ., N − 1.We also defined the spin-wave density, The equations of motion for the spin-wave correlators are derived from the Heisenberg equations generated by the sum of the quadratic Hamiltonians Û 0 2 + Û2 and read as follows: We refer the reader to Refs.[37,38] for a more detailed calculation.We observe that, in the limit of J → 0, the equations ( 25) decouple from the quantum feedback and consistently reduce to a representation of Eq. ( 5) in the spherical coordinates θ(t) and ϕ(t).
The NEQSWT is expected to be valid as long as the density of spin-wave excitations is small, that is ϵ(t) ≪ 1.In this case, the modulus of the magnetization | ⃗ S(t)| = 1 − ϵ(t) is close to one and the dynamics can still be described in terms of classical trajectories.In this regime, spin waves behave as free bosonic excitations, which interact with the macroscopic collective spin only, corresponding to the k = 0 mode.Higher-order terms appearing in Eq. ( 1), which account for nonlinear scattering among the spin waves, can be neglected: they are expected to contribute significantly to the dynamics only at longer times and to drive the system away from the prethermal regime relevant for the DPT discussed here [38].

B. Modified non-equilibrium phase diagram
The dynamics of the magnetization, as described by Eqs.(25) , is still equivalent to the one of a particle moving into a multi-well shaped energy profile [60], as discussed in Sec.II A. However, in this case, an additional damping effect arises from the exchange of energy with the spin-wave degrees of freedom, essentially acting as a self-generated bath.Further details on this mechanism are discussed in Sec.III C. The strength of the damping is controlled by the coupling J in Eqs.(25).For p = 2, this mechanism is responsible for the melting of the dynamical critical point into a chaotic crossover phase [37,38].In this phase, the magnetization asymptotically localizes into one of the two ferromagnetic wells, although the dynamics is strongly sensitive to both the initial conditions and the integration parameters.It is natural to ask if the spin-wave emission is going to have the same effect in the first-order case p > 2.
To understand the effect of fluctuations on the mean-field dynamical transition, we study the post-quench dynamics of the system at J > 0, by simultaneously integrating equations ( 25) and ( 29), for several values of the couplings g and J. Before the quench, we prepare again the system in the fully polarized state from Eq. ( 2), which is equivalent to the initial conditions ϕ(0) = 0 and θ(0) = π/2.In this initial state, there are no spin-wave excitations.For each choice of g and J and within our maximum simulation time, we verify the spin wave density ϵ(t) from Eq. ( 28) always settles around a small value, so that the dynamics is consistently in a prethermal regime.In particular, as shown in Fig. 5, ϵ(t) typically grows from zero and saturates to a small value.Consequently, the magnetization is asymptotically damped to a trajectory whose energy, H cl ( ⃗ S(t), g), is slightly lower than the one at t = 0. We define the new dynamical phases in terms of the topology of these trajectories, asymptotically reached after the damping.It is crucial to note once again that these new phases are prethermal : they are expected to be disappear at longer times, as soon as the non-linear interaction among the spin-wave degrees of freedom is taken into account, leading to thermal relaxation of the system.Below we will reconstruct the dynamical phase diagrams by looking at S x as a function of g and J and at individual trajectories.Similar to Sec.II, our focus will be on the cases of p = 3 and p = 4, which are paradigmatic for p > 2 odd and p > 2 even, respectively.If spin waves were emitted at a constant rate in time, the orbits encompassing all the minima (i.e.orange trajectories in Fig. 2-(a) and (d) ) would localize down one of the wells of H cl ( ⃗ S, g) with approximately equal probability.The phase diagram for S x , depicted for p = 3 in Fig. 3, reveals a more nuanced behavior.While the dynamical ferromagnetic phase is of course robust against fluctuations, the dynamical bistable phase loses stability for sufficiently large values of J and its trajectories predominantly localize around the paramagnetic minimum, asymptotically converging to stable orbits.These asymptotic orbits identify a third, new dynamical paramagnetic phase on the phase diagram.Quantitatively, each dynamical phase corresponds to a narrow interval of values of S x , with the greatest values in the dynamical ferromagnetic phase and the smallest in the paramagnetic one (although non-vanishing due to the asymmetry of the energy profile).In our discussion, we use the notation J dyn (g) to identify the transition line between the dynamical bistable and dynamical paramagnetic phases.As detailed in Appendix B, S x exhibits a discontinuity when crossing J dyn (g), giving rise to a new first-order DPT driven by J.The predominance of localization around S x = 0 is softened close to the mean-field critical point g = g dyn , where localization in the ferromagnetic basin becomes more frequent.The line separating the dynamical ferromagnetic phase from other phases is consequently melted in a narrow chaotic crossover region, akin to the one observed for p = 2 [37].
As shown in the phase diagrams from Fig. 4, a similar phenomenon is observed for p = 4. Here, as soon as J is moved above a critical threshold J dyn (g) (generically distinct from the one retrieved for p = 3), the trajectories from the mean-field dynamical paramagnetic phase, initially encompassing all the three minima of the landscape,  28) and computed by within the NEQSW described in the main text, by fixing either p = 3, J = 0.25, g = 1.15 (red plot) or p = 4, J = 0.2, g = 1.12 (blue plot).In both plots we fix N = 100 and λ = 1.
become unstable and localize in the paramagnetic well.However, this discontinuity in the time-evolution can not be observed from S x , displayed in Fig. 4 (left).Specifically, S x is zero for orbits either surrounding all the minima or asymptotically localizing in the paramagnetic well, while S x > 0 in the ferromagnetic basin.Thus, we also examine the behaviour of the time-averaged fluctuation [42], defined as From the phase diagram in Fig. 4 (right), it is clear that (δS x ) 2 is discontinuous across the transition line J dyn (g), where the topology of the asymptotic trajectories abruptly changes.As a consequence, the dynamical paramagnetic phase, corresponding to S x = 0, can be divided in two sub-phases: -A dynamical paramagnetic phase 1, where the asymptotic orbits surround all the minima, like in the corresponding mean-field phase.
-A dynamical paramagnetic phase 2, where the trajectories eventually localize in the paramagnetic well.
Like in the case of p = 3, these phases are identified by different narrow intervals of values of (δS x ) 2 .The smaller values correspond to the dynamical paramagnetic phase 2 and to eventual localization in the paramagnetic basin.

C. The mechanism behind the localization of the magnetization
Our results are very different from the ones found in Ref. [37] for the case of p = 2, where the same short-range perturbation from Eq. ( 14) led to dynamical chaos.In this section we show that, despite the apparent difference, the dynamical phases retrieved both for p = 2 and p > 2 share a common origin, which can be related to the predominant emission of spin waves when the classical trajectory visits the paramagnetic well.
Intuitively, the inhomogeneous emission is due to the specific form of the short-range perturbation in Eq. ( 14), which induces fluctuations in the collective dynamics when the local spins are not aligned with the x-axis: as the maximal misalignment in the dynamics is reached when the magnetization crosses the plane S x = 0, spin waves are expected to be mostly emitted there.The previous argument can be made quantitative by investigating the time evolution of the spin-wave density ϵ(t), which quantifies the degree of spin-wave excitations in the system.By differentiating both sides of Eq. ( 28), we obtain the spin-wave emission rate: Then, substituting Eqs. ( 29) on the right-hand side we obtain the following evolution equation for ϵ(t): The quantum feedback terms δ αβ (t) are the ones defined in Eq. ( 26).
From Eq. ( 32), it is evident that the emission rate dϵ/dt depends explicitly on the position of the magnetization on the Bloch sphere.This dependency is encapsulated in the coefficients: We observe that the coefficients C 1 (θ, ϕ) and C 2 (θ, ϕ) do not depend on p. Instead, they are derived from a combination of the coefficients appearing in the Hamiltonian Û2 from Eq. ( 24), obtained by a quadratic expansion of the perturbation Û , from Eq. ( 14), in the moving frame R.This observation aligns with the physical expectation that Û is the only term in the Hamiltonian which generates spin-waves excitations.Both C 1 (θ, ϕ) and C 2 (θ, ϕ) vanish when the magnetization is along the x-axis, i.e. θ = π/2 and ϕ = 0.In Fig. 6, we plot the time evolution of the amplitudes |C 1 (θ, ϕ)| and |C 2 (θ, ϕ)|, along two sample trajectories from the dynamical paramagnetic phases retrieved for p = 3 and p = 4, respectively.The results shown therein confirm that these amplitudes are maximised when the magnetization crosses the plane S x = 0, located in the paramagnetic well.The previous observations confirm our intuition that the maximum spin-wave emission occurs concomitantly with the maximal misalignment between ⃗ S(t) and x.
As the predominant emission in the paramagnetic well is determined only by the short-range perturbation, the different phenomena observed for p = 2 and p > 2 respectively can be addressed to the different stability properties of paramagnetic the stationary point ⃗ S = z of the energy landscape in Eq. (6).For p = 2, ⃗ S = z is unstable and symmetric fluctuations in the two wells [61] induce dynamical chaos [37,38].However, for p > 2, ⃗ S = z becomes a stable minimum so that the magnetization moves on stationary orbits after being damped, giving birth to the dynamical paramagnetic regions shown in Fig. 3 (left) and Fig. 4 (right).This phenomenon is reminiscent of Hopf bifurcations [62] occurring in classical dynamical systems.
It is also worth noticing that the mechanism by which the fluctuations induce the localization of the collective spin is slightly more subtle than a simple dissipation mechanism.In particular, we observe that when ϵ(t) > 0, the magnetization length is | ⃗ S(t)| = 1 − ϵ(t) decreases, so that the dynamics in Eq. ( 25) takes place in the time-dependent modified potential As shown in the animated plots (see Ancillary Files [63]), the profile of the H ϵ(t) (θ, ϕ) is squeezed towards zero energy when ϵ(t) grows: this eventually leads the magnetization to be trapped in the paramagnetic region, where ϵ(t) exhibits large spikes, while the ϵ(t) is nearly stationary across the ferromagnetic wells.
All the results presented in this section are expected to be independent of the range of interaction of the perturbation in Eq. ( 14) and of the dimensionality of the lattice: replacing the k-dependent couplings J cos(k) with generic Jk [64] leaves Eq. ( 32) invariant [65] and spin-waves are still expected to be emitted around the plane S x = 0, as explicitly shown for the case of p = 2 [38].On the other hand, if the short-range interaction was along a direction not coinciding with the x-axis, the dependency of the coefficients from Eq. ( 33) on the angles θ and ϕ could be different: in this case, spin-wave emission may occur in different regions on the Bloch sphere, leading to a different non-equilibrium phase diagram.

IV. CONCLUSIONS
In conclusion, in this work we have studied the the post-quench dynamics of a fully-connected p-spin model (for p > 2) perturbed by a short-ranged interaction, controlled by the coupling J, generalizing to arbitrary values of p the system studied in previous work [29,35,37].In the mean-field limit of J = 0, the dynamics is equivalent to the one of a particle in a classical energy landscape.By identifying the topology of each orbit as a dynamical phase, we observe a dynamical phase transition (DPT) driven by the transverse field g.The order of the DPT depends on the parity of p, for p > 2, which determines the qualitative shape of the effective potential.In particular, we find a second order dynamical transition for odd values of p (where the potential is an asymmetric double-well) and first order one for an even p (where the profile is made by three basins).The nature of this transition is modified by the presence of a short-range perturbation, treated in the framework of the non-equilibrium spin-wave theory.The latter generates an effective damping in the dynamics which is maximal when the magnetization visits the paramagnetic well.The damping induces a prethermal stage of the dynamics and changes the nature of the dynamical phases, now being defined by the asymptotic behaviour of the magnetization.For p ≥ 3 odd, the short-range coupling drives a new first-order transition on the phase diagram, while a more subtle transition appears for p ≥ 4 even, detected in the order parameter fluctuations and being of the first-order like the one obtained in the mean-field limit.Our analysis can be straightforwardly generalized to a wider class of fully-connected spin models with generic integrability breaking terms: the profile of the energy landscape and the axis where the integrability breaking interaction takes place are the only two ingredients which fully determine the features of the non-equilibrium phase diagram.FIG. 8. Plot of the singularity of the time-average magnetization S x , respectively for p = 3 (blue) and p = 4 (red).S x is plotted as a function of g in the dynamical ferromagnetic phase, where g < g dyn , and S x c is defined as an unstable stationary point of the potential in Eq. ( 6) of the main text at g dyn .We refer the reader to Appendix A for further details.
with a prefactor which is different above and below g dyn .
The divergence of T cl (g) at g dyn is connected to change of topology of the underlying trajectories, being confined in a single ferromagnetic well for g < g dyn and exploring the whole landscape for g > g dyn .In particular, at g = g dyn the trajectory is a separatrix, that is a singular, non-periodic orbit.Any dynamics evolving on the separatrix asymptotically converges to the nearest local maximum ⃗ S c , separating the well where the motion takes place from the others.For any value of p, this implies that S x (g dyn ) = S x c .Then, the divergence of T cl (g) also generates a singularity in S x , for g → g − dyn .In this case, the orbits evolving into the rightmost ferromagnetic well develop a plateau of diverging length near S x c (see Fig. 2-(c) and (f) of the main text), so that S x can be estimated as Eq. (A8) holds under the reasonable assumption that the integral T cl (g) 0 dt S x (t) − S x c is bounded and converges to a positive constant c > 0 at the transition point.Eq. (A8) implies that when approaching the transition from below.The log-singularity is quantitatively confirmed by the results shown in Fig. 8, obtained by computing S x numerically (as in Sec.II of the main text) and S x c from Eqs. ( 8) and ( 10) of the main text.We remark that, while the log-singularity in the time-average magnetization S x are intimately connected to the one of the periods, its continuity is determined by the topology of the effective potential, as discussed in the main text.smooth J above and below the discontinuity point, because of the noise induced by the spin-wave emission.In this regime, it is possible to make a semi-analytical estimation of the point J dyn (g) where the transition happens, as discussed in the following.First we remind that, as discussed in Sec.III C, the effective potential determining the dynamics is squeezed by the emission of spin waves and is expressed by Eq. ( 34) of the main text.From this expression, it is evident that the dynamics is driven by an effective transverse field Here, g ϵ depends on time through the spin-wave density ϵ(t), defined in Eq. ( 28).However, Eq. ( 9) also indicates the existence of a critical value of g, denoted as g sp , such that for g > g sp the energy landscape becomes a single well.These observations leads us to the conclusion that it exist a threshold value ϵ sp for the spin-wave density, defined by the equation g (1 − ϵ) p−1 = g sp , (B2) such that the squeezed potential from Eq. ( 34) displays a single paramagnetic well whenever ϵ > ϵ sp .Then, a sufficient condition for the localization of the magnetization in the paramagnetic well corresponds to ϵ(t) asymptotically exceeding ϵ sp .Quantitatively, this corresponds to: Eq. (B3) condition implies that, for large times, the squeezed potential from Eq. ( 34) displays a single, paramagnetic well, where the magnetization is by definition localized.The plots in Fig. 9 (left) show that, fixing a sufficiently large g, the values of J where ϵ overcomes the threshold ϵ sp matches the transition point J dyn (g): this provides a semianalytical argument for the prediction of J dyn (g) at large g, as the threshold ϵ sp can be predicted analytically, but there is no explicit formula relating J dyn (g) to ϵ.Our argument fails for smaller g, where the localization mechanism becomes more subtle (as discussed in Sec.III C the main text) and the first-order transition line is smeared into a chaotic crossover close to the mean-field critical point g dyn .A similar fate happens to the second-order critical line, extending from the point (g, J) = (g dyn , 0), to finite values of J, so that the chaotic crossover prevents the possibility of a precise estimation of the tricritical point where the two lines meet.However, we can roughly identify the transition point as the center of the finite-width area where the two transition lines merge completely with the chaotic region: in the inset in Fig. 9 (bottom-right), we show that this happens approximately around (g, J) ≃ (1.026, 0.17).For p = 4, the results in Fig. 10 (right) show a discontinuous transition driven by J, this time detected by the time-averaged fluctuations (δS x ) 2 , even though the semi-analytical estimation of J dyn (g) fails in this case.Moreover, here both the mean-field dynamical transition, driven by g, as well as the one driven by J, are of the first-order, so that we do not retrieve the tricritical behaviour discussed for p = 3.

FIG. 1 .
FIG. 1. (Color online) Energy profiles described by Eq. (6), projected on the plane S y = 0.The projection is obtained from the black cut on the Northern hemisphere of the corresponding Bloch spheres, where S z = 1 − (S x ) 2 .Each panel illustrates the shape of the energy landscape corresponding to a different set of values of p: p = 2 (left), p = 2n + 1 (center) and p = 2n + 2 (right), respectively, with n ≥ 1 integer.For practical purposes, we fixed n = 1.Each profile in each panel corresponds to a different value of g: g = g dyn − 0.1 (a), g = g dyn + 0.1 (b) and g = g dyn (c).The value of g dyn is obtained from Eqs. (8) and (10), as discussed in the main text.For each energy profile, the corresponding trajectory evolving from ⃗ S(0) = x (black dot on the profile) is represented by a dashed horizontal line.In particular, the gray dashed horizontal line corresponds to a separatrix (white trajectory on the Bloch sphere).See also the orbits in Fig.2(a) and (d) for a comparison.

FIG. 2 .
FIG. 2.(Color online) From left to right: orbits of the magnetization on the Bloch sphere ((a) and (d)), time-averaged magnetization S x calculated from Eq. (11) as a function of g ((b) and (e)), time evolution of S x (t) near the critical point ((c) and (f)).The plots in the upper row corresponds to the case of p = 3, the ones below to p = 4. (Left) Orbits on the Bloch sphere, resulting from the integration Eq. (5), from the initial condition ⃗ S = x.Each trajectory corresponds to one of 20 values of g,uniformly sampled in the interval [0.2, 3.6], where g dyn falls for both p = 3 and p = 4.For both values of p, the system evolves along trajectories that are either confined in a ferromagnetic well (green trajectories), for g < g dyn , or encompass all the local minima (orange trajectory), for g > g dyn .The two regions are divided by a separatrix (white), corresponding to g = g dyn .g dyn is calculated from Eqs. (8) and(10),as discussed in the main text.(Center) Plot of the dynamical order parameter S x , from Eq. (11), as a function of g.The red dot correspond to the dynamical critical point (g dyn , S x c ), obtained by simultaneously solving simultaneously solving Eqs.(8) and(10).At the dynamical critical point, S x is continuous for p = 3 only.(Right) Time evolution of the longitudinal magnetization S x (t), for a transverse field g either equal to g dyn − 10 −7 (green plots) or to g dyn + 10 −7 (orange plots).For each value of p, the dashed lines correspond to the local maxima of the landscape when g = g dyn .For each value of p, g dyn and S x c are obtained as in the central panels.

FIG. 3 .
FIG. 3. (Color online) (Left) Non-equilibrium phase diagram for the time-average magnetization S x , from Eq.(11).Each point of the phase diagram is obtained corresponds to a value of S § , obtained by simultaneously Eqs.(25) and (29), for a given choice of the couplings g and J.The initial condition is fixed at ⃗ S = x.The results shown in this figure refer to the case of p = 3.We fix the maximum simulation time as T = 2000.Other simulation parameters are fixed at N = 100 and λ = 1.The color of each point on the phase diagram corresponds to a value of the time-averaged magnetization S x , as specified in the interval shown on the right of the diagram.In particular, the yellow and the orange regions correspond to dynamical ferromagnetic and paramagnetic phases, respectively.(Right) The corresponding orbits are either confined in the ferromagnetic well (plot T1) or encircle both the two minima (plot T2).In the blue region, instead, is the dynamical paramagnetic region where the magnetization revolves around the minimum on the North pole (plot T3); here S x is closer to 0 than in any other phase, though non-vanishing, and increases discontinuously when moving across the border with the orange zone.The crossover region is instead a narrow chaotic phase, reminiscent of the one found in[37], where collective spin can localize at large times in either the paramagnetic or the ferromagnetic minima (plot T4).

FIG. 4 .
FIG. 4. (Color online) Simultaneous plot of the time-averaged magnetization S x (left) and time-averaged fluctuations (δS x ) 2 (right), for p = 4.These are obtained integrating equations Eq.(25) and (29) for the same default parameters of N ,λ and the maximum time T listed in the caption of Fig.3.For both the plots, each color corresponds to a value of the observable we plot, as specified in the legends reported on the right of each diagram.(Left) From the plot of S x , we identify two main regions: one corresponds to a dynamical ferromagnetic phase(yellow), where S x > 0 (plot T1), while the other is a paramagnetic phase (orange) where S x = 0.In between, some chaotic spots (blue and yellow spots) are found occasionally, where the magnetization eventually falls in one of the other two symmetric, ferromagnetic minima (plot T4).(Right) Looking at the average fluctuations (δS x ) 2 , we clearly see that the in dynamical paramagnetic phase, coinciding with the non-blue region of the phase diagram, can be split in a sub-phase 1 (yellow), where the dynamical trajectories either surround symmetrically all the three minima of the landscape (plot T2) and a sub-phase 2 (blue), where the magnetization localizes (predominantly) in a paramagnetic well (plot T3).The border between the two identifies the transition line J = J dyn (g).

FIG. 9 .
FIG. 9. (Color online) Non-equilibrium phase diagram for the p-spin model in eq.(1), for p = 3. (Left) Over the same phase diagram shown in Fig. 3 (left) of the main text, we plot the threshold line from Eq. (B1) (red), above which ϵ > ϵsp (see details in Appendix B).The black vertical line and the green box indicate the values of g and J investigated in the plots on the right.(Top-right) Plot of the time-averaged magnetization S x , as function of J and at fixed g ≃ 1.142 (black line on the phase diagram on the left).S x is discontinuous around J ≃ 0.194.(Bottom-right) Inset from the phase diagram on the left (in the green box), around the chaotic region where the second-order and first-order transition line merge.

FIG. 10 .
FIG. 10. (Color online) Non-equilibrium phase diagram for the p-spin model in eq.(1), for p = 4. (Left) Same phase diagram shown in Fig. 4 (left).The black vertical line and the green box indicate the values of g and J investigated in the plots on the right.(Right) Plot of the time-averaged fluctuations (δS x ) 2 , as function of J and at fixed g ≃ 1.133 (black line on the phase diagram on the left).(δS x ) 2 is discontinuous around J ≃ 0.131.