Time delayed control of excited state quantum phase transitions in the Lipkin–Meshkov–Glick model

We theoretically investigate the role of dissipation in excited state quantum phase transitions (ESQPT) within the Lipkin–Meshkov–Glick model. Signatures of the ESQPT are directly visible in the complex spectrum of an effective Hamiltonian, whereas they get smeared out in the time-dependence of system observables. In the latter case, we show how delayed feedback control can be used to restore the visibility of the ESQPT signals.


I. INTRODUCTION
The properties of a quantum phase transition (QPT) [1] and an excited state QPT (ESQPT) [2,3] driven by quantum fluctuations in many body quantum systems at zero temperature feeds the interest since many decades.At a quantum level, the ESQPT is hidden in the system's density spectrum as discontinuities or divergence [2,4].For some famous models like Dicke superradiance [5] or Lipkin-Meshkov-Glick (LMG) [6], the main properties of both QPTs can be explained and understood at a semiclassical level in the thermodynamic limit, where a QPT corresponds to a bifurcation [7,8] and a ESQPT is connected to a saddle point in the semi-classical energy potential [2].Furthermore, the density of states can then be calculated analytically [4,9].For both models, the ESQPT manifests in the observable average as a peak at a certain energy [4,10,11] and was already observed in molecular systems [12] or microwave Dirac billiards [13].
Feedback control is a promising tool to change the system dynamics in a desired way.Influence of the laser statistics [14], neurosystems [15] or even control at a quantum level [16][17][18][19][20] are only some examples of its powerfulness.One usually distinguishes between the closed and open control loops, in the last one the feedback depends on the state of the system, where time delayed Pyragas control [21] has an important niche.For instance, it was used to speed up the convergence to a steady state in the dissipative quantum system [22].
In our previous works, we applied time delayed Pyragas control to a Dicke model [5] to create new nonequilibrium phases [23] and suggested a new method to extract the ESQPT-signal from the time evolution in the closed LMG system [6,24].Furthermore, closed loop control was already applied to the LMG model to induce new phases or to modify the divergence in density of states [25][26][27].Moreover, Oberthaler et al. implemented the LMG model using the existing ingenious experimental setup based on 87 Rb Bose-Einstein condensates [29], which offer a high degree of freedom for system preparation, providing a possibility to test new insights about the (ES)QTP.Surprisingly, the effects of a dissipative environment on the ESQPT seems not to be well studied in contrast to the QPT [30], though they are always present in experimental realisations.
Inspired by this, we study the effects of dissipation on the ESQPT signal in the LMG model and apply a time delayed feedback scheme to cancel them by the creation of new phases.In our model, we condition the atomic coupling on the difference of a spin observable average at two different times and perform the calculation at a mean-field-level.On top we show, that in dissipative systems the ESQPT is directly visible from a complex system spectrum.
Our work is organised as follows.In section II we introduce the dissipative LMG model and the feedback.In section III we study the dissipative effects on a ESQPT at a quantum level evaluating the complex spectrum of an effective Hamiltonian.In section IV we show the smoothing effects of the ESQPT signal due to dissipation at a mean-field level and show in section V how to compensate them using time delayed Pyragas control.In the last VI section we discuss the results.

II. LMG MODEL WITH DISSIPATION AND FEEDBACK
The LMG model [6] describes an interaction between N spins and represents a special case of a Heisenberg Model.In general, the LMG Hamiltonian reads where Ĵi = N j=1 σj i , i ∈ {x, y, z} are collective angular momentum operators, h is an effective parameter for external magnetic field in z-direction and γ x or γ y describes the spin-spin interaction strength, which is the same for all spins in the LMG model.In the following, we will always use the isotropic γ y = 0 case, if we do not explicitly point to γ y = 0 case.Especially in the experimental realizations, the LMG system is always coupled to an environment which cause the damping and thermalisation of the system and can be modelled by a master equation [30,31] with collective decay [32] where D[ Ĵ+ ]ρ = Ĵ− Ĵ+ ρ + ρ Ĵ− Ĵ+ − 2 Ĵ+ ρ Ĵ− is the Lindblad-Dissipator and Ĵ± = Ĵx ± i Ĵy .
To compensate the dissipative effects we assume a time dependent coupling γ x of Pyragas form [21] .Therefore we condition γ x to depend on a difference of Ĵz averages at two different times with time delay τ as following In the following, we investigate the ESQPT signal at the quantum level using the effective Hamiltonian approach and at a semiclassical level using the solution of mean-field equations in thermodynamic limit.In the second case we show the unique effect of the feedback loop especially in context of the ESQPT signal.

III. DISSIPATIVE ESQPT AT QUANTUM LEVEL
The fate of the ESQPT signal in presence of dissipation is still not studied well.In closed systems the ESQPT is hidden in the energy spectrum or is visible in observable averages as a function of energy [10] as a non-analyticity, which can be obtained by quantum mechanical or semiclassical calculations.Its origin can be understood at a semi-classical level as critical points from the energy surface.
In dissipative systems described by a Lindblad master equation, an effective non-hermitian Hamiltonian can be defined in a standard way.Rewriting the master equation Eq. (2) as we obtain the effective non-hermitian Hamiltonian with complex spectrum, which is shown in the Fig. 1 (upper) for different values of κ in the symmetry broken phase.For κ = 0 the spectrum has an imaginary part which scales with κ.The imaginary part can be interpreted as a decay rate at a certain energy level [33,34].To complete the spectral information we show in the lower part of Fig. 1 the corresponding density of states.
For example, open triangles show the density of states along the blue dashed line in the upper figure in the known κ = 0-case [9], where a logarithmic divergence at E = −0.5N is due to the ESQPT.
How does the dissipation affect the ESQPT?The ES-QPT survives and, somewhat surprisingly, it becomes visible not only in the density of states of the nonhermitian Hamiltonian H ef f , but can be also directly seen from its complex eigenvalues, see Fig. 1

(upper).
The ESQPT is hidden in this representation in a -form at the energy of E = −0.5N .Another feature (due to the assumed Lindblad operator (2)) is the vanishing of the imaginary part of H ef f at the north and south poles of the corresponding Bloch sphere, which leads to zero imaginary part at the corresponding energies of ∓0.5N .Thus the decay rate is zero there and the dissipative effects disappear at this points.Note, this effect is also present at the level of the mean-field Eqs.(7), where dis-sipative terms do not contribute at the poles for arbitrary κ values due to conservation of the spin length.
We used two different methods to calculate the density of states in the κ = 0.05h case.First, we used only the real part of the eigenvalues and counted their number in a certain energetic window (blue crosses in Fig. 1).Second, we used with the non-hermitian Hamiltonian (open squares).We emphasize, that the results agrees and have still a logarithmic divergence at the energy E = −0.5Nwhich can be well fitted by a log-function (red dotted curve).Especially for E > −0.5N there is no strong deviation from the κ = 0 case (open triangles).Only for E < −0.5N , there is a deviation from the non-dissipative result, which can be better seen for the curve with a bigger dissipation rate κ = 0.5 (green dot-dashed).
We emphasize, that only the -structure and not the zero imaginary part in the complex spectrum indicates the ESQPT.For γ y = 0 both effects are at the energy E = −0.5, though they can be easily separated for an anisotropic LMG Hamiltonian (γ y = 0), where the ES-QPT can be shifted to other energies and a jump in a density of states can occur on top [9].Setting γ y = 2.5h, we observe a shift of the -signal to the corresponding energy of ESQPT in this case (open circles in the upper part of Fig. 1 ).The corresponding density of states (open circles in the lower part of the figure) has a peak at the ESQPT energy E ≈ −0.54N and a jump at E = −0.5N[9].Note, that the finite size effects smooths the peak and the jump in this case.

IV. DISSIPATIVE ESQPT SIGNAL AT MEAN-FIELD LEVEL A. Mean-field-equations
To derive a mean-field equation for the dissipative LMG system, we use Ô = T r( Ô ρ), assume the factorisation assumptions Ô1 Ô2 ≈ Ô1 Ô2 which is known to hold in the thermodynamic limit and to forecast the same observable averages and the phase transition as the quantum mechanical calculations [24,35].We then obtain a following set of closed semiclassical equations of motion [30] Jx (t) = hJ y (t) − κJ x (t)J z (t), ( 7) with rescaled averages J i = 1 N Ĵi .Without time delay (λ = 0) the QPT, which is one important property in this system, corresponds at this semi-classical level to a pitchfork bifurcation [31]: Eq.( 7) has two stationary states J 0 i , corresponding to a normal phase with (J 0 x , J 0 y , J 0 z ) = (0, 0, 1/2) and a symmetry broken phase whose stability swaps at a critical coupling Note, that even the dissipative model still fulfils the conservation law J 2 x +J 2 y +J 2 z = 1/4, thus the dynamical evolution is restricted to a sphere.Furthermore, even with time delay (λ > 0) the fixed points remain the same, as the Pyragas term vanishes in the steady state.
In the following, we investigate the ESQPT signal in presence of damping with and without control.Later we show that the time delayed coupling γ(t) may affect the linear stability of fixed points and the dynamical evolution of the system in a completely unexpected way, particular in its acting against the dissipation.

B. Dissipative damping of the ESQPT
Using the semiclassical equation of motion Eq. ( 7) we now study the action of dissipation on the ESQPT at this level.Therefore we first look at a dynamical evolution of the LMG system.The spin averages of the system are restricted to the Bloch sphere, which is shown in Fig. 2 (inset).The color represents the energy for a given system configuration for κ = 0 and the white lines represent the paths with the same energy in the symmetry broken phase.Without dissipation the system follows one of this paths keeping the energy fixed, i.e. staying in the eigenstate.For each eigenstate one can compute an averaged value for (J 2 x , J z ).The black (dotted) curve (Fig. 2) shows this values for multiple eigenstates.This is a novel representation of two observable averages, which was recently suggested to visualize the energy independent ES-QPT signal [24], which is visible again as a peak.The continuation of the peak would end by (0, 1/2).In the upper half of a Bloch sphere we see a separatrix, a (white) path which goes through a north pole of the sphere.Due to the symmetry of the path, its averaged values for J 2 x and J 2 y should be zero, whereas the J z average is 1/2.Thus, paths close to the separatrix are responsible for the ESQPT peak in (J 2 x , J z ) diagram.Without dissipation (κ = 0) the period length for fixed energy (Fig. 2, upper right) diverges at the separatrix energy [2].For κ = 0 the energy is not conserved any more and the system tends oscillating around the J z -axis to a steady state Eq.(8).Fig. 2 (inset) shows two examples of the system state evolution (J x , J y , J z )(t) for two different κ values, which where obtained by solving Eq. ( 7) numerically.The ESQPT signal is now hidden in the dynamical evolution of the system.But, as we can see in Fig. 2, especially for big κ values there are only less paths, which are close to the separatrix, thus the ESQPT signal will be damped especially for big dissipation rates κ.The impact of damping to the ESQPT signal in the J z , J 2 x plane can be obtained by calculating the mean values of J 2 x (t), J z (t)-evolution for a period or by finding some optimal effective period ∆t, in a way that the mean values calculated using the definition matches to the closed case as good as possible.Note, as a period length in the first case we take the time for one full rotation around the J z -axis, which changes with time.
In the first case, the results for different damping rates κ are shown in Fig. 2 by colored circles/diamonds and unfilled squares.We see, that the peak is now smoothed, but still visible especially for small damping rates.The blue crosses show the second case with effective period which lead to a much better ESQPT signal.

V. PYRAGAS CONTROL OF THE ESQPT SIGNAL
Now we set λ = 0 in γ x (t) Eq. ( 3) and investigate the impact of time delayed control on the system.

A. Linear stability analysis with time delay
A usual approach to analyse the effects of time delayed feedback is first to check the stability of fixed points in the presence of control [36].Therefore we linearise Eq. ( 7) around the fixed points J i (t) = J 0 i + δJ i , obtaining the following system of linearised equations with δv ≡ (δJ x , δJ y ) T , with with eliminated J z component by use of the spin-length conservation law.
The roots of the corresponding characteristic equation det(Λ1 − B − A exp(−Λτ )) = 0 (13) determines the stability of a fixed point, which is stable if all real parts of all solutions Λ are negative [37].
In Fig. 3 we plot the biggest real part of eigenvalues in the τ − λ− domain.We see that there is a window for λ -values, there the stability of fixed points oscillates from stable to unstable and from unstable to stable while increasing the time delay.Outside this window, the fixed points remains either stable (λ 1h) or lose their stability forever (λ 2.5h).Note, that the boundaries between stable and unstable zones (brown line) can be calculated from an analytical expression, see Appendix A. Next, we analyse the system properties in the unstable regime, use them to obtain a sharp ESQPT signal and show chaotic behaviour for larger time delays τ .

B. Feedback compensates dissipation
Our feedback scheme Eq. (3) modifies the system dynamics in an interesting and unexpected way.Increasing the time delay τ , we cross the boundary and make the fixed point unstable.For smaller τ values the trajectory ends in a new stable state in form of a limit cycle, thus a Hopf bifurcation occurs.Moreover the trajectory of the limit cycle has only small deviations from paths with fixed energy, which the LMG system would take without dissipation.The size of trajectories can be changed by τ , thus tuning the time delay value corresponds to a change of energy in a closed LMG system.Fig. 4 demonstrates the feedback action, showing the trajectory evolution for different values of τ .Here we start close to the fixed point.The red (thick) curve shows the stationary state.Note, that the change of initial condition to the region outside the separatrix leads to the same stationary state for the considered τ values.Thus, using the stationary limit cycle states for fixed τ values offers a possibility to obtain an ESQPT-Signal again, by calculating the and J2 x (with similar definition) averages for fixed values of τ (t 2 > t 1 and with t 1 big enough to become a stationary solution).Fig. 5 shows the results of this calculation and compares them with the closed and dissipa-tive cases without feedback.The (orange) rhombi shows the Jz , J2 x time averages in stationary limit cycle phase of the dissipative LMG system, the black curve shows the ESQPT -signal of a closed system without control.Each rhombus has its own time delay, the inset in Fig. 5 shows the τ -dependence of the averaged values.Note, that up to τ ≈ 0.2/h the fixed point is stable and the averaged values are the values of the fixed point.We see a very good overlap between the control caused and the original signals, thus our feedback schema compensates the dissipative smoothing (unfilled green squares) of the ESQPT-signal very well.

C. Chaotic behaviour
For τ 1 the stationary dynamics can become much more complex than just a creating limit cycles.Oscillations with more than one maximum and minimum appear, which is known as a way to chaos by period doubling [38,39].In Fig. 6 (upper) we plot the maxima and minima of J x (t) oscillations for t 1 as a function of time delay τ and mark different areas by capital letters A-E.In the inset we show a zoom for the area A. Increasing τ from zero, the stable fixed points (arrow) loses their stability (orange dotted line) at the first boundary condition from the Fig. 3 and a Hopf-bifurcation appears (part A).In part B, the fixed points become stable again.But there exist still a stable limit cycle solution with rather big J x -amplitude.Thus, in area B two stable solutions are possible, which is not in contradiction to the stability of the fixed point, as the initial condition is not chosen in a way to fulfil the linearised assumption.For 4 τ 5 this limit cycle disappears.Further increase of time delay leads to a creation of a period doubling structure (area C), which is separated by windows where the solution converges to a limit cycle.In the region D the fixed point becomes stable again and the double period structure is gone, whereas in the Region E it appears again.Note, that the dotted line at J max x = 0 represents the unstable trivial fixed point.
Such chaotic behaviour can also be used to obtain the ESQPT signal, see blue crosses in Fig. 5. Fixing the time delay in the chaotic area of the phase diagram (Fig. 6), one can use the integration method with effective period Eq. ( 10) to obtain the shown curve from the dynamical system evolution on the Bloch sphere.We see, that also this result matches pretty well to the original ESQPT signal.

VI. DISCUSSION
In this paper we have demonstrated the effect of dissipation on the ESQPT signal for the LMG model and showed how to compensate it using time delayed Pyragas feedback modulating the interaction parameter between the atoms.Our results show, that the ESQPT is encoded x in the stationary state for fixed τ with τ ∈ [0.1, 0.6]h, which matches to the signal of the closed system (black curve) very well.Unfilled squares represent the ESQPT-signal under dissipation.The blue crosses shows the ESQPT-signal obtain from a dynamical evolution of the system for one big τ -value using Eq. ( 14).(Inset) Jz, J2 x as a function of τ in a steady state.Parameters: in the spectral properties of the effective Hamiltonian as well as it is visible in the measured averaged values of the spin components.In the last case, smoothing effects appear which can be undone by our feedback scheme.
We think, that an experimental measurement of an ESQPT signal using the time delayed method is easier than in an ideal, closed system.Using time delay one has only to measure the system values for different time delays τ , instead of preparing the system in eigenstates or coherent superpositions to obtain the same information [24].
We also checked the interaction strength γ x depending on other operator differences instead of J 2 z , or a modulated magnetic field h instead of γ x .However, in both cases the general dynamical properties remains the same.The LMG system has still limit cycles close to the separatrix for some fixed time delayed values and has a param- eter range with chaotic behaviour.Both effects might be interesting from an experimental point of view.On the one hand, it is easier to control the magnetic field h, on the other hand, choosing another feedback loop can shift the appearing effects to other τ values, which could be easer to realize.
We also tried to find ESQPT signatures in the correlation function [40] C t (z) ≡ Ĵ+ (t + z) Ĵ− (t) , using the quantum regression theorem [41], but have not succeed, as the ESQPT signal (which should be visible as a peak at zero frequency in the Fourier space) interfere with effects like macroscopic occupation and degenerated spectrum in the symmetry broken phase.However one possible way out could be to drive the LMG system additionally with an external laser and calculating the resonance fluorescence spectrum, but this would require a re-derivation of the master equation.Note, that such spectra have been already calculated for an LMG system, but only in the linearised version [31].
Similar to the Dicke model with feedback [23], the Pyragas controlled LMG model has stable limit cycle phases.In contrast to the Dicke model, it shows chaotic dynamics for bigger time delay.This is surprisingly as one would expect such behaviour especially from the Dicke system, which is chaotic in nature [42].However,the Dicke model has an additional bosonic mode which is not bounded by a conservation law like the spin components.This could be a reason, why no chaotic behaviour appears there.
We think, that the feedback-induced limit cycles are hidden property of the LMG modell.On the one hand limit cycles are a natural property of the closed system.On the other hand, the time dependent limit cycles describe the ESQPT signal pretty well.Though the chaotic behaviour is an artefact of the time delay feedback, as it was neither a part of a closed LMG system.
The feedback scheme is applied at a semiclassical level, thus we have neglected the influence of fluctuations in the thermodynamic limit, which could be important [43].Nevertheless, we think, that the fluctuations have not dramatic contributions to the described effects for N 1 as they should scale with 1/ √ N and the used feedback scheme does not modify γ in a strong way.Furthermore, the semiclassical LMG model predicts in many cases the same results [4,10,24].But the full quantum version of the considered feedback type still remains an open issue.However, a recently published article [20] shows a way to go beyond mean field for a coherent type of feedback where the author describes feedback action via mapping to a bigger system.This would be one possible way to study the role of oscillations in quantum systems with one special feedback type.
Bringing all sin and cos terms in both equations to one side, squaring them and adding together, we eliminate the τ -dependence and obtain the following equation with 2 −2h 2 (J 0 z ) 2 + 4h (J 0 x ) 2 2(J 0 z ) 3 λ − γJ 0 z + γ(J 0 z ) 3 + κ 2 (J 0 x ) 4 − 4κ(J 0 x ) 3 J 0 y γ − 2(J 0 z ) 2 λ + 2(J 0 x ) 2 (J 0 y ) 2 κ 2 + 2γ γ − 4(J 0 z ) which can be solved for s and has then in general 4 different solutions The Eq. (A5) fixes now the eigenvalue Λ = is for a given fixed point and feedback strength λ.Eq. (A1) gives for every fixed s-value the corresponding time-delay τ , solving for example the imaginary part for τ we obtain where z ∈ Z.The choice of z is necessary to get boundary conditions at higher time delays τ .Note, that the s dependence is also hidden in G i .At this step we have still to much solutions.A lot of them are non-physical (if τ < 0 or τ ∈ C) or do not fulfil the real part equation Eq. (A1) and have to be sorted out.The remaining solutions are plotted as a brown line in Fig. 3.

FIG. 1 .
FIG. 1. (upper) The complex spectrum for different values of γ.The -structure indicates the position of the ESQPT.The inset shows the zoom around the ESQPT energy.(lower) The density of states of the effective LMG Hamiltonian shows a logarithmic divergence at E = −0.5N .For comparison, the γy = 0 case is additionally shown in both figures (open circles).Parameters: (γ = 1.5h),N = 1000 − 10000,λ = 0.

FIG. 3 .
FIG. 3. Stability diagram in the τ − λ-plane.The color represents the stability robustness of fixed points in the symmetry broken phase, which are not stable for positive values (red coloured area).The brown line defines a boundary condition and has an analytical expression.Parameters: γ = 1.5h, κ = 0.05h.

FIG. 4 .
FIG.4.Dynamical evolution of the controlled LMG system on the Bloch sphere for the following τ -values (from left to right): τ h = 0.2; 0.25; 0.3; 0.31; 0.35; 0.5; 1; 2.5; The red thick path represents the stationary state.Increase of τ forces the solution to cross the separatrix.In this way the ESQPT signal (Fig.5) is restored from stationary solutions with different τ values.Note, the axes and scaling are as in the inset of Fig.2but are not visible due to breakdown.Parameters: γ = 1.5h, λ = 1, κ = 0.05h.

FIG. 6 .
FIG. 6. Bifurcation and period doubling for different time delays τ .For each τ -values the values of local maxima and minima of Jx(t)-value are plotted in stationary case.The dashed orange lines shows the unstable fixed points.The arrows points to the non-trivial fixed points.The phase diagram is divided the regions A-E with different properties.The inset shows shows a zoom of area A. Parameters: γ = 1.5h, λ = 1h.