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

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 QPT is visible as a non-analyticity of the ground state at some point in the parameter space [1]. In contrast the ESQPT is a quantum criticality of the excited states, hidden as a level clustering [2] in the energy spectrum. Therefore, the ESQPT becomes visible in the system spectrum density as discontinuities or divergence [2, 4], which induce a non-analyticity for some system observables, e.g. for the average over the corresponding eigenstates [5]. Moreover, the eigenstate structure contains the ESQPT signatures, too [6]. For some famous models like Dicke superradiance [7–10] or Lipkin–Meshkov–Glick (LMG) [11], 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 [5, 12] and a ESQPT is connected to a saddle point in the semi-classical energy potential [2]. A critical path with an energy of the saddle point—the so called separatrix—can then be defined in the semiclassical configuration space. The paths which are close to the separatrix constructs then the ESQPT signal: for both models, the ESQPT manifests in the observable average as a peak at a certain energy [4, 13, 14]. Furthermore, a semiclassical treatment allows the analytical calculation of the density of states [4, 15]. ESQPT was for example studied in connection with quenches [5, 6] or chaos [5] and was already observed in molecular systems [16] or microwave Dirac billiards [17].

Feedback control is a promising tool to change the system dynamics in a desired way. Influence of the laser statistics [18], neurosystems [19] or even control at a quantum level [20–24] 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 [25] has an important niche. For instance, it was used to speed up the convergence to a steady state in the dissipative quantum system [26].

In our previous works, we applied time delayed Pyragas control to a Dicke model [7] to create new non-equilibrium phases [27] and suggested a new method to extract the ESQPT-signal from the time evolution in the closed LMG system [11, 28]. Furthermore, closed loop control was already applied to the LMG model to induce new phases or to modify the divergence in density of states [29–31]. Moreover, Oberthaler et al implemented the LMG model using the existing ingenious experimental setup based on ${}^{87}\mathrm{Rb}$ Bose–Einstein condensates [32], 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 [33], though they are always present in experimental realizations.

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 organized as follows. In section 2 we introduce the dissipative LMG model and the feedback. In section 3 we study the dissipative effects on a ESQPT at a quantum level evaluating the complex spectrum of an effective Hamiltonian. In section 4 we show the smoothing effects of the ESQPT signal due to dissipation at a mean-field level and show in section 5 how to compensate them using time delayed Pyragas control. In the last 6 section we discuss the results.

The LMG model [11] describes an interaction between N spins and represents a special case of a Heisenberg model. In general, the LMG Hamiltonian reads

$$\begin{eqnarray}&&\hat{H}=-h{\hat{J}}_{z}-\displaystyle \frac{{\gamma }_{x}}{N}{\hat{J}}_{x}^{2}-\displaystyle \frac{{\gamma }_{y}}{N}{\hat{J}}_{y}^{2},\end{eqnarray} \tag{ 1 }$$

where ${\hat{J}}_{i}=\frac{1}{2}{\displaystyle \sum }_{j=1}^{N}{\hat{\sigma }}_{i}^{j},\quad i\in \{x,y,z\}$ are collective angular momentum operators, h is an effective parameter for external magnetic field in z-direction and ${\gamma }_{x}$ or ${\gamma }_{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 ${\gamma }_{y}=0$ case, if we do not explicitly point to ${\gamma }_{y}\ne 0$ case.

Especially in the experimental realizations, the LMG system is always coupled to an environment which cause the damping and thermalization of the system and can be modelled by a master equation [33, 34] with collective decay [35]

$$\begin{eqnarray}&&\dot{\hat{\rho }}=-{\rm{i}}[\hat{H},\hat{\rho }]-\displaystyle \frac{\kappa }{2N}{\mathcal{D}}[{\hat{J}}^{+}]\hat{\rho },\end{eqnarray} \tag{ 2 }$$

where ${\mathcal{D}}[{\hat{J}}^{+}]\hat{\rho }={\hat{J}}^{-}{\hat{J}}^{+}\hat{\rho }+\hat{\rho }{\hat{J}}^{-}{\hat{J}}^{+}-2{\hat{J}}^{+}\rho {\hat{J}}^{-}$ is the Lindblad-dissipator and ${\hat{J}}^{\pm }={\hat{J}}_{x}\pm {\rm{i}}{\hat{J}}_{y}.$

To compensate the dissipative effects we assume a time dependent coupling ${\gamma }_{x}$ of Pyragas form [25] . Therefore we condition ${\gamma }_{x}$ to depend on a difference of ${\hat{J}}_{z}$ averages at two different times with time delay τ as following

$$\begin{eqnarray}&&{\gamma }_{x}={\gamma }_{x}(t)=\gamma +\displaystyle \frac{\lambda }{{N}^{2}}({\left\langle{\hat{J}}_{z}(t-\tau )\right\rangle}^{2}-{\left\langle{\hat{J}}_{z}(t)\right\rangle}^{2}).\end{eqnarray} \tag{ 3 }$$

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.

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 [13] as a non-analyticity, which can be obtained by quantum mechanical or semi-classical 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 (2) as

$$\begin{eqnarray}&&\dot{\hat{\rho }}=-{\rm{i}}[{\hat{H}}_{\mathrm{eff}},\hat{\rho }]+\displaystyle \frac{\kappa }{N}{\hat{J}}^{+}\hat{\rho }{\hat{J}}^{-},\end{eqnarray} \tag{ 4 }$$

we obtain the effective non-hermitian Hamiltonian

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{eff}}=\hat{H}-{\rm{i}}\displaystyle \frac{\kappa }{2N}{\hat{J}}^{-}{\hat{J}}^{+}\end{eqnarray} \tag{ 5 }$$

with complex spectrum, which is shown in the figure 1 (left) for different values of κ in the symmetry broken phase. For $\kappa \ne 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 [36, 37]. To complete the spectral information we show in the right part of figure 1 the corresponding density of states. For example, open triangles show the density of states along the blue dashed line in the left figure in the known $\kappa =0$-case [15], where a logarithmic divergence at $E=-0.5N$ is due to the ESQPT. Note, that in some references like [28] the normalization of the system observables is performed by $j\equiv N/2$ instead of here used N normalization, the ESQPT is then located at $E=-1j.$

Zoom In Zoom Out Reset image size

How does the dissipation affect the ESQPT? The ESQPT survives and, somewhat surprisingly, it becomes visible not only in the density of states of the non-hermitian Hamiltonian H_eff, but can be also directly seen from its complex eigenvalues, see figure 1 (left). 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_eff at the north and south poles of the corresponding Bloch sphere, which leads to zero imaginary part at the corresponding energies of $\mp 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 equations (7), where dissipative 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 $\kappa =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 figure 1). Second, we used

$$\begin{eqnarray}&&\nu (E)=-1/\pi {\mathfrak{I}}\mathrm{Tr}\left(\displaystyle \frac{1}{E-{H}_{\mathrm{eff}}}\right)\end{eqnarray} \tag{ 6 }$$

with the non-hermitian Hamiltonian (open squares). We emphasize, that the results agrees and have still a logarithmic divergence at the energy $E=-0.5N$ which can be well fitted by a log-function (red dotted curve). Especially for $E\gt -0.5N$ there is no strong deviation from the $\kappa =0$ case (open triangles). Only for $E\lt -0.5N,$ there is a deviation from the non-dissipative result, which can be better seen for the curve with a bigger dissipation rate $\kappa =0.5$ (diamonds with horizontal line).

We emphasize, that only the ⋎-structure and not the zero imaginary part in the complex spectrum indicates the ESQPT. For ${\gamma }_{y}=0$ both effects are at the energy $E=-0.5,$ though they can be easily separated for an anisotropic LMG Hamiltonian $({\gamma }_{y}\ne 0),$ where the ESQPT can be shifted to other energies and a jump in a density of states can occur on top [15]. Setting ${\gamma }_{y}=2.5\;h,$ we observe a shift of the ⋎-signal to the corresponding energy of ESQPT in this case (open circles in the left part of figure 1 ). The corresponding density of states (open circles in the right part of the figure) has a peak at the ESQPT energy $E\approx -0.54N$ and a jump at $E=-0.5N$ [15]. Note, that the finite size effects smooths the peak and the jump in this case.

4.1. Mean-field-equations

To derive a mean-field equation for the dissipative LMG system, we use $\left\langle\hat{O}\right\rangle=\mathrm{Tr}(\hat{O}\hat{\rho }),$ assume the factorization assumptions $\left\langle{\hat{O}}_{1}{\hat{O}}_{2}\right\rangle\approx \left\langle{\hat{O}}_{1}\right\rangle\left\langle{\hat{O}}_{2}\right\rangle$ 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 [28, 38]. We then obtain a following set of closed semiclassical equations of motion [33]

$$\begin{eqnarray}{\dot{J}}_{x}(t) & = & {{hJ}}_{y}(t)-\kappa {J}_{x}(t){J}_{z}(t),\\ {\dot{J}}_{y}(t) & = & -h{J}_{x}(t)+2{\gamma }_{x}(t){J}_{x}(t){J}_{z}(t)-\kappa {J}_{y}(t){J}_{z}(t),\\ {\dot{J}}_{z}(t) & = & -2{\gamma }_{x}(t){J}_{x}(t){J}_{y}(t)+\kappa \left({J}_{x}{(t)}^{2}+{J}_{y}{(t)}^{2}\right),\end{eqnarray} \tag{ 7 }$$

with rescaled averages ${J}_{i}=\frac{1}{N}\left\langle{\hat{J}}_{i}\right\rangle.$ Without time delay $(\lambda =0)$ the QPT, which is one important property in this system, corresponds at this semi-classical level to a pitchfork bifurcation [34]: equation (7) has two stationary states ${J}_{i}^{0},$ corresponding to a normal phase with $({J}_{x}^{0},{J}_{y}^{0},{J}_{z}^{0})=(0,0,1/2)$ and a symmetry broken phase

$$\begin{eqnarray}{J}_{x}^{0} & = & \pm \sqrt{\displaystyle \frac{-4\;{h}^{2}(\gamma -\sqrt{{\gamma }^{2}-{\kappa }^{2}})+{\kappa }^{2}(\gamma +\sqrt{{\gamma }^{2}-{\kappa }^{2}})}{8\gamma {\kappa }^{2}}},\\ {J}_{y}^{0} & = & \displaystyle \frac{\gamma -\sqrt{{\gamma }^{2}-{\kappa }^{2}}}{\kappa }\cdot {J}_{x}^{0},\quad {J}_{z}^{0}=h\cdot \displaystyle \frac{\gamma -\sqrt{{\gamma }^{2}-{\kappa }^{2}}}{{\kappa }^{2}},\end{eqnarray} \tag{ 8 }$$

whose stability swaps at a critical coupling

$$\begin{eqnarray}&&{\gamma }_{{\rm{c}}}=h+\displaystyle \frac{{\kappa }^{2}}{4h}.\end{eqnarray} \tag{ 9 }$$

Note, that even the dissipative model still fulfils the conservation law ${J}_{x}^{2}+{J}_{y}^{2}+{J}_{z}^{2}=1/4,$ thus the dynamical evolution is restricted to a sphere. Furthermore, even with time delay ($\lambda \gt 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 $\gamma (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.

4.2. Dissipative damping of the ESQPT

Using the semiclassical equation of motion, equation (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 figure 2 (inset). The colour represents the energy for a given system configuration for $\kappa =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}_{x}^{2},{J}_{z}).$ The black (dotted) curve (figure 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 ESQPT signal [28], 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_x² and J_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}_{x}^{2},{J}_{z})$ diagram. Without dissipation $(\kappa =0)$ the period length for fixed energy diverges at the separatrix energy [2].

Zoom In Zoom Out Reset image size

For $\kappa \ne 0$ the energy is not conserved any more and the system tends oscillating around the J_z-axis to a steady state equation (8). Figure 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 equation (7) numerically. The ESQPT signal is now hidden in the dynamical evolution of the system. But, as we can see in figure 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}_{x}^{2}$ plane can be obtained by calculating the mean values of ${J}_{x}^{2}(t),{J}_{z}(t)$-evolution for a period or by finding some optimal effective period ${\rm{\Delta }}t,$ in a way that the mean values calculated using the definition

$$\begin{eqnarray}&&\bar{O}(t)=\displaystyle \frac{1}{{\rm{\Delta }}t}{\displaystyle \int }_{t=a}^{a+{\rm{\Delta }}t}O(t){\rm{d}}t,\quad a\in \{0,{T}_{\mathrm{max}}\},\end{eqnarray} \tag{ 10 }$$

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 figure 2 by coloured 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.

Now we set $\lambda \ne 0$ in ${\gamma }_{x}(t),$ equation (3) and investigate the impact of time delayed control on the system.

5.1. 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 [39]. Therefore we linearize equation (7) around the fixed points ${J}_{i}(t)={J}_{i}^{0}+\delta {J}_{i},$ obtaining the following system of linearized equations with $\delta {\bf{v}}\equiv {(\delta {J}_{x},\delta {J}_{y})}^{T},$

$$\begin{eqnarray}&&{\partial }_{t}\delta {\bf{v}}(t)={\bf{B}}\cdot \delta {\bf{v}}(t)+{\bf{A}}\cdot \delta {\bf{v}}(t-\tau ),\end{eqnarray} \tag{ 11 }$$

with

$$\begin{eqnarray}{\bf{A}} & \equiv & -4\lambda {J}_{z}^{0}\left(\begin{array}{cc}0 & 0\\ {({J}_{x}^{0})}^{2} & {J}_{x}^{0}{J}_{y}^{0}\end{array}\right),\\ {\bf{B}} & \equiv & \left(\begin{array}{cc}\kappa \displaystyle \frac{{({J}_{x}^{0})}^{2}}{{J}_{z}^{0}}-\kappa {J}_{z}^{0} & h+\kappa \displaystyle \frac{{J}_{x}^{0}{J}_{y}^{0}}{{J}_{z}^{0}}\\ 4\lambda {J}_{z}^{0}{({J}_{x}^{0})}^{2}-2\gamma \displaystyle \frac{{({J}_{x}^{0})}^{2}}{{J}_{z}^{0}}+\kappa \displaystyle \frac{{J}_{y}^{0}{J}_{x}^{0}}{{J}_{z}^{0}}-h+2\gamma {J}_{z}^{0} & \kappa \displaystyle \frac{{({J}_{y}^{0})}^{2}}{{J}_{z}^{0}}+4\lambda {J}_{x}^{0}{J}_{z}^{0}{J}_{y}^{0}-2\gamma \displaystyle \frac{{J}_{x}^{0}{J}_{y}^{0}}{{J}_{z}^{0}}-\kappa {J}_{z}^{0}\end{array}\right),\end{eqnarray} \tag{ 12 }$$

with eliminated J_z component by use of the spin-length conservation law.

The roots of the corresponding characteristic equation

$$\begin{eqnarray}&&\mathrm{det}({\rm{\Lambda }}{\bf{1}}-{\bf{B}}-{\bf{A}}\mathrm{exp}(-{\rm{\Lambda }}\tau ))=0\end{eqnarray} \tag{ 13 }$$

determines the stability of a fixed point, which is stable if all real parts of all solutions Λ are negative [40].

In figure 3 we plot the biggest real part of eigenvalues in the $\tau -\lambda -$ 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 $(\lambda \ll 1h)$ or lose their stability forever $(\lambda \gtrsim 2.5h).$ Note, that the boundaries between stable and unstable zones (brown line) can be calculated from an analytical expression, see the appendix. 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 τ.

Zoom In Zoom Out Reset image size

5.2. Feedback compensates dissipation

Our feedback scheme (equation (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. Figure 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.

Zoom In Zoom Out Reset image size

Thus, using the stationary limit cycle states for fixed τ values offers a possibility to obtain an ESQPT-Signal again, by calculating the

$$\begin{eqnarray}&&{\bar{J}}_{z}=\displaystyle \frac{1}{{t}_{2}-{t}_{1}}{\displaystyle \int }_{{t}_{1}}^{{t}_{2}}{J}_{z}(t){\rm{d}}t\end{eqnarray} \tag{ 14 }$$

and ${\bar{J}}_{x}^{2}$ (with similar definition) averages for fixed values of τ (${t}_{2}\gt {t}_{1}$ and with t₁ big enough to become a stationary solution). Figure 5 shows the results of this calculation and compares them with the closed and dissipative cases without feedback. The (orange) rhombi shows the ${\bar{J}}_{z},{\bar{J}}_{x}^{2}$ 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. 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. The inset in figure 5 shows the τ-dependence of the averaged values. We have checked that this dependency is qualitatively the same as the energy dependency of the corresponding expectation values in eigenstates with the energy E in the closed system without time delay. Note, that up to $\tau \approx 0.2/h$ (inset) the fixed point is stable and the averaged values are the values of the fixed point.

Zoom In Zoom Out Reset image size

5.3. Chaotic behaviour

For $\tau \gg 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 [41, 42]. In figure 6 (left) we plot all local maxima and minima between which J_x(t) oscillates in the stationary state (for $t\gg 1$) for fixed values of time delay τ. Note, that the steady state of the system can depend on the initial condition, too, which has to be changed, to get all possible end states. First, we show an example how to obtain the data-points which are located on the vertical dash–dotted line in the bifurcation diagram. Figure 6 (right) shows possible end states of the system in the $({J}_{x},{J}_{y})$ plane for different initial conditions. The maxima and minima of the J_x oscillation (arrows) are lying on the dash-dotted line in the bifurcation diagram. Next, we discuss the systems dynamic as a function of τ, therefore we mark different areas in the bifurcation diagram by capital letters A–E. In the inset we show a zoom for the area A, which contains three parts A1–A3. Increasing τ from zero, both stable fixed points (arrow, orange line, part A1) loses their stability (orange dotted line) at the first boundary condition from the figure 3 and a Hopf-bifurcation appears (part A). There are two limit cycle possible, for positive and negative J_x values, however both limit cycles grow with τ and merge to one big limit cycle which covers both J_x sides on the bloch sphere, part A3 (see also figure 4 for visualization). 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 linearized assumption. The system will then converge either to the fixed point or to the limit cycle, depending on the initial condition. For $4\lesssim \tau \lesssim 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}_{x}^{\mathrm{max}}=0$ represents the unstable trivial fixed point.

Zoom In Zoom Out Reset image size

Such chaotic behaviour can also be used to obtain the ESQPT signal, see blue crosses in figure 5. Fixing the time delay in the chaotic area of the bifurcation diagram (figure 6, left), one can use the integration method with effective period (equation (10)) to obtain the shown curve. Therefore we solve the corresponding equation of motion (7) and use equation (10) starting from t = 0 till $t\gg 1.$ We see, that also this result matches pretty well with the original ESQPT signal.

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 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 [28].

We also checked the interaction strength ${\gamma }_{x}$ depending on other operator differences instead of ${J}_{z}^{2},$ or a modulated magnetic field h instead of ${\gamma }_{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 parameter 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 [43]${C}_{t}(z)\equiv \left\langle{\hat{J}}^{+}(t+z){\hat{J}}^{-}(t)\right\rangle,$ using the quantum regression theorem [44], 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 linearized version [34].

Similar to the Dicke model with feedback [27], 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 [45]. 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 [46]. Nevertheless, we think, that the fluctuations have not dramatic contributions to the described effects for $N\gg 1$ as they should scale with $1/\sqrt{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, 13, 28]. But the full quantum version of the considered feedback type still remains an open issue. However, a recently published article [24] 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.

We thank Georg Engelhardt and Mathias Hayn for useful discussions. The authors gratefully acknowledge financial support from the DAAD and DFG Grants BR 1528/7-1, 1528/8-2, 1528/9-1, SFB 910, and GRK 1558.

To determine the boundaries in figure 3 we choose only imaginary Λ values in equation (13), thus ${\rm{\Lambda }}\equiv i\cdot s,s\in {\mathbb{R}}.$ Splitting then the equation (13) in an imaginary and real part, we obtain the correspondingly conditions

$$\begin{eqnarray}0 & = & {G}_{0}+{G}_{1}\cdot \mathrm{cos}(s\tau )+{G}_{2}\cdot \mathrm{sin}(s\tau ),\\ 0 & = & {G}_{3}+{G}_{4}\cdot \mathrm{cos}(s\tau )+{G}_{5}\cdot \mathrm{sin}(s\tau ),\\ {G}_{0} & = & -\displaystyle \frac{\kappa {({J}_{x}^{0})}^{2}\;s}{{J}_{z}^{0}}+\displaystyle \frac{2\gamma {J}_{x}^{0}{J}_{y}^{0}\;s}{{J}_{z}^{0}}-4{J}_{x}^{0}{J}_{y}^{0}{J}_{z}^{0}\lambda s-\displaystyle \frac{\kappa {({J}_{y}^{0})}^{2}\;s}{{J}_{z}^{0}}+2\kappa {J}_{z}^{0}s,\\ {G}_{1} & = & 4{J}_{x}^{0}{J}_{y}^{0}{J}_{z}^{0}\lambda s,\quad {G}_{2}=-4h{({J}_{x}^{0})}^{2}{J}_{z}^{0}\lambda -4\kappa {J}_{x}^{0}{J}_{y}^{0}{({J}_{z}^{0})}^{2}\lambda ,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}{G}_{3} & = & {h}^{2}-{s}^{2}+\displaystyle \frac{2\gamma h{({J}_{x}^{0})}^{2}}{{J}_{z}^{0}}-2\gamma {{hJ}}_{z}^{0}-{\kappa }^{2}{({J}_{x}^{0})}^{2}-{\kappa }^{2}{({J}_{y}^{0})}^{2}+{\kappa }^{2}{({J}_{z}^{0})}^{2}-4h{({J}_{x}^{0})}^{2}{J}_{z}^{0}\lambda -4\kappa {J}_{x}^{0}{J}_{y}^{0}{({J}_{z}^{0})}^{2}\lambda ,\\ {G}_{4} & = & 4h{({J}_{x}^{0})}^{2}{J}_{z}^{0}\lambda +4\kappa {J}_{x}^{0}{J}_{y}^{0}{({J}_{z}^{0})}^{2}\lambda ,\quad {G}_{5}=4{J}_{x}^{0}{J}_{y}^{0}{J}_{z}^{0}\lambda s.\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&0={s}^{4}+{F}_{1}{s}^{2}+{F}_{0},\end{eqnarray} \tag{ 3 }$$

with

$$\begin{eqnarray*}{F}_{0} & = & \displaystyle \frac{1}{{({J}_{z}^{0})}^{2}}\left\{\left({h}^{2}{J}_{z}^{0}+2\gamma h\left({({J}_{x}^{0})}^{2}-{({J}_{z}^{0})}^{2}\right)+{\kappa }^{2}{J}_{z}^{0}\left({-({J}_{x}^{0})}^{2}-{({J}_{y}^{0})}^{2}+{({J}_{z}^{0})}^{2}\right)\right)\right.\\ & & \cdot \left({h}^{2}{J}_{z}^{0}+2h\left({({J}_{x}^{0})}^{2}\left(\gamma -4{({J}_{z}^{0})}^{2}\lambda \right)-\gamma {({J}_{z}^{0})}^{2}\right)\right.\\ & & +\left.\left.\kappa {J}_{z}^{0}\left(-\kappa {({J}_{x}^{0})}^{2}-8{J}_{x}^{0}{J}_{y}^{0}{({J}_{z}^{0})}^{2}\lambda +\kappa \left({({J}_{z}^{0})}^{2}-{({J}_{y}^{0})}^{2}\right)\right)\right)\right\},\\ {F}_{1} & = & \displaystyle \frac{1}{{({J}_{z}^{0})}^{2}}\left\{-2{h}^{2}{({J}_{z}^{0})}^{2}+4h\left({({J}_{x}^{0})}^{2}\left(2{({J}_{z}^{0})}^{3}\lambda -\gamma {J}_{z}^{0}\right)+\gamma {({J}_{z}^{0})}^{3}\right)+{\kappa }^{2}{({J}_{x}^{0})}^{4}\right.\\ & & -4\kappa {({J}_{x}^{0})}^{3}{J}_{y}^{0}\left(\gamma -2{({J}_{z}^{0})}^{2}\lambda \right)+2{({J}_{x}^{0})}^{2}\left({({J}_{y}^{0})}^{2}\left({\kappa }^{2}+2\gamma \left(\gamma -4{({J}_{z}^{0})}^{2}\lambda \right)\right)-{\kappa }^{2}{({J}_{z}^{0})}^{2}\right)\\ & & -4\kappa {J}_{x}^{0}{J}_{y}^{0}\left({({J}_{y}^{0})}^{2}\left(\gamma -2{({J}_{z}^{0})}^{2}\lambda \right)+2{{J}_{z}^{0}}^{4}\lambda -2\gamma {({J}_{z}^{0})}^{2}\right)\\ & & +\left.{\kappa }^{2}\left({({J}_{y}^{0})}^{4}-2{({J}_{y}^{0})}^{2}{({J}_{z}^{0})}^{2}+2{({J}_{z}^{0})}^{4}\right)\right\},\end{eqnarray*}$$

which can be solved for s and has then in general 4 different solutions

$$\begin{eqnarray}&&s=\pm \displaystyle \frac{1}{\sqrt{2}}\sqrt{-{F}_{1}\pm \sqrt{{F}_{1}^{2}-4\cdot {F}_{0}}}.\end{eqnarray} \tag{ 4 }$$

The equation (4) fixes now the eigenvalue ${\rm{\Lambda }}={is}$ for a given fixed point and feedback strength λ. equation (1) gives for every fixed s-value the corresponding time-delay τ, solving for example the imaginary part for τ we obtain

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1}{s}\mathrm{arctan}\left(\displaystyle \frac{\displaystyle \frac{\pm \displaystyle \frac{{G}_{1}\sqrt{-{G}_{2}^{2}\left({G}_{0}^{2}-{G}_{1}^{2}-{G}_{2}^{2}\right)}}{{G}_{1}^{2}+{G}_{2}^{2}}+\displaystyle \frac{{G}_{0}{G}_{1}^{2}}{{G}_{1}^{2}+{G}_{2}^{2}}-{G}_{0}}{{G}_{2}}}{\displaystyle \frac{\pm \sqrt{-{G}_{0}^{2}{G}_{2}^{2}+{G}_{1}^{2}{G}_{2}^{2}+{G}_{2}^{4}}-{G}_{0}{G}_{1}}{{G}_{1}^{2}+{G}_{2}^{2}}}\right)+\displaystyle \frac{2\pi }{s}\cdot z,\end{eqnarray} \tag{ 5 }$$

where $z\in {\mathbb{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 $\tau \lt 0$ or $\tau \in {\mathbb{C}}$) or do not fulfil the real part equation (1) and have to be sorted out. The remaining solutions are plotted as a brown line in figure 3.