Rapid Steady State Convergence for Quantum Systems Using Time-Delayed Feedback Control

We propose a time-delayed feedback control scheme for open quantum systems that can dramatically reduce the time to reach steady state. No measurement is performed in the feedback loop, and we suggest a simple all-optical implementation for a cavity QED system. We demonstrate the potential of the scheme by applying it to a driven and dissipative Dicke model, as recently realized in a quantum gas experiment. The time to reach steady state can then reduced by two orders of magnitude for parameters taken from experiment, making previously inaccessible long time attractors reachable within typical experimental run times. The scheme also offers the possibility of slowing down the dynamics, as well as qualitatively changing the phase diagram of the corresponding physical system.


Introduction
Steady states of open quantum systems, where driving forces and internal dynamics are balanced by dissipation and/or other types of environmental noise, are often of experimental interest. There are indeed a number of platforms currently available, where the experimental control of the individual constituents of interacting quantum systems allows precise preparation of interesting and non-trivial steady states, through measurement and control of well-defined outputs and inputs of the system. These include systems using trapped ions or ultra-cold atoms, opto-mechanical systems, and systems based on cavity or circuit quantum electrodynamics (CQED) (see e.g. Refs. [1][2][3]). This can be an alternative to quantum state preparation by coherent (unitary) evolution, and is potentially of great interest to quantum information processing technologies as a way of preparing computational resources, such as maximally entangled states [4][5][6][7], or even providing for a route to quantum computation [8]. One great advantage of such an approach is that the steady state is robust against variations of the initial state.
In practice, the time-scale on which such a steady state is reached is very important; the generation of the desired state often requires a degree of control that is hard to sustain over time, which can pose a challenge for finite-time experiments [9][10][11]. In the present paper, we propose an all-optical feedback scheme relevant, for example, to CQED systems, that can be used to i) change the stability of long time attractors, so that one can switch between different behaviors, and ii) change the characteristic time-scale for approaching a steady state, thus potentially speeding up the convergence. The scheme we use is based on the time-delayed feedback control method developed by Pyragas [12], and often referred to as the time-delay autosynchronization (TDAS) [13,14]. The control is based on coherent feedback, i.e., no measurement is performed in the feedback loop, which can be advantageous, or even necessary, for stabilizing the high frequency dynamics of optical systems or high speed electrical circuits [13,14]. Another great strength of the approach is that it does not require the steady state to be known a priori. We apply coherent TDAS, to our knowledge, for the first time to a quantum system, with the feedback signal treated quantum mechanically as well [15][16][17]. Delay-times have often been assumed to be negligible in theoretical modelling of coherent quantum feedback, motivated by the fact that a time-delay can introduce undesirable instability to the system, and that it can often be made very small in practice [15,18]. In contrast, with TDAS, the delay is the crucial ingredient for increasing the system's stability. We note that delayed coherent feedback has also been considered in [19,20] for a single-atom single-excitation system (and under the usual rotating wave approximation).
We will demonstrate the potential of our scheme by applying it to a highly topical example, namely an open system version of the Dicke model. The Dicke model is a paradigmatic model in quantum optics, describing the interaction of a collection of two-level atoms with a single cavity mode [21]. This model has been realized and studied in a number of recent experiments [9,22,23] based on a Bose-Einstein condensate (BEC) coupled to a field mode of a high finesse optical cavity. The cavity has a natural dissipative output channel, which has been used to monitor the system in real time. In particular, the system undergoes a quantum phase transition as an effective coupling strength between the BEC and the light field is increased beyond a critical value, which can be observed through the intensity of the output cavity field [9]. Spontaneous symmetry breaking has also been observed through a heterodyne measurement scheme [22], and measuring the correlations of the density fluctuations has been used to observe the diverging time-scale upon approaching the critical point [23]. This type of monitoring of a dissipative channel is non-destructive, and as these experiments have shown, offers a very promising route for the observation of complex many-body quantum dynamics.
We will take advantage of this dissipative channel as well, by using it as the input to a non-invasive feedback loop that can alter the characteristic time-scale for the relaxation of the system, as well as the stability of the long-time attractors. This is particularly relevant for the BEC realization of the Dicke model where, as we will discuss in more detail below, the approach to steady state can be slow compared to typical experimental run times. Adverse effects such as spontaneous emission, or atom loss, will eventually cause deviations from the desired, idealized behavior. This is particularly problematic close to phase boundaries, where we expect critical slowing down. The system is therefore a very interesting test-case for our proposed feedback scheme. We will treat the feedback in a semi-classical approach, where quantum fluctuations are linearized. For the Dicke model, this approximation is valid in the thermodynamic limit of a large number of atoms. The approach we develop should be similarly applicable to a variety of topical quantum-optical systems for which feedback can also have useful and interesting consequences (see e.g. Refs. [24][25][26][27]).
The paper is organized as follows. In Section 2 we introduce our feedback scheme in a general setting, using the standard input-output theory for quantum optical systems. Then, in Section 3, we introduce our primary system of study, the Dicke model as recently experimentally realized, and apply our feedback scheme. We study in detail the effect of the forcing due to the feedback, and find optimal delay-times for rapid convergence to steady state. We compare the performance of the system with and without feedback, and demonstrate improvements in the relaxation time of two orders of magnitude. In Section 4 we consider potential consequences for finite-time experiments. In Section 5 we consider the influence of the feedback force on quantum fluctuations. Finally, in Section 6, we give some concluding remarks.

Coherent time-delay auto-synchronization
The Pyragas' time-delay auto-synchronization (TDAS) method is a continuous feedback control method first developed to stabilize unstable periodic orbits and equilibrium states embedded in a chaotic attractor [12][13][14]28]. We will also use the method for manipulating the characteristic time-scale on which a stable fixed point is approached. Briefly, the idea behind TDAS for stabilizing a dynamical system, whose classical state is given by x(t), is to apply one or more continuous feedback forces of the form The feedback force vanishes in steady state, or for a periodic orbit if the delay, τ , is a multiple of the period. This is referred to as noninvasive feedback. The delay, and feedback strength, k, are parameters that should be varied in experiment to achieve driving towards a particular long-time attractor.
We consider a CQED system, as illustrated in Fig. 1. The internal dynamics of a cavity, consisting of a field mode, possibly interacting with other quantum and classical degrees of freedom, is described by a HamiltonianĤ. The cavity is assumed to have two mirrors b and c, corresponding to two distinct pairs of input-output ports. The equation of motion for the cavity mode is, according to the standard input-output Figure 1. A schematic illustration of our feedback control scheme. The internal dynamics of the cavity is described by a Hamiltonian,Ĥ. The cavity consists of two mirrors, b (left) and c (right), with decay rates κ b and κc respectively. The light-blue tilted bars denote beam splitters BS1 and BS2, with transformation properties defined by the unitary matrices S 1 and S 2 , respectively. The circle denotes a delay-time of τ and a φ phase shift to one of the feedback arms. We assume that Faraday isolators (not shown) separate the cavity input and output fields.
theory for optical quantum systems [29], Hereâ is the cavity mode annihilation operator, κ b and κ c are the decay rates for the two mirrors, andb in (t),ĉ in (t) the annihilation operators of the input fields incident on the respective mirrors. Here, and in the following, we will, when convenient, suppress the time argument for any system operator evaluated at time t, but keep the time argument for input and output fields for clarity. We will assume thatĉ in (t) corresponds to a vacuum field, although it would also be of interest to consider a drive here. The input fields obey the commutation relations whereξ in denotes any of the input mode operatorsb in orĉ in . In addition, any vacuum input fieldξ in (t) corresponds to a Gaussian white noise operator with zero mean, and its only non-vanishing correlation function is The corresponding output fields are given aŝ andĉ The output from mirror c is now split into two feedback arms, as illustrated in Fig.  1. One of the feedback arms has a time-delay τ and a phase shift φ, while we set the time-delay and phase shift of the other arm to zero for simplicity. The input and output fields for the two beam-splitters BS1 and BS2 shown in Fig. 1 are related through f 1 (t) where S 1 and S 2 are unitary matrices,ν 1 (t) is a vacuum input field to beam splitter BS1, andν 2 (t) is the (unused) other output field from beam splitter BS2. We want the time-delayed and time-undelayed fields to be incident on mirror b in opposite phase, so as to give the desired destructive interference in steady state. For this purpose, we choose the beam splitter transformations to be given by and where r, s ≥ 0 (real) and r 2 + s 2 = 2. In passing, we remark that these are not the most general choice of beam splitter transformations but sufficient for our purposes. We now find forb in (t):b where we have definedb The mode operatorb in (t) satisfies Eq.
(2) and Eq. (3) and is thus the "vacuum part" of the field incident on mirror b. A control force is generated from the difference between the current cavity field, a(t), and the field at some point in the past,â(t − τ ). This forcing is then fed back into the system. By making use of Eq. (7) in Eq. (1) we obtain where k ≡ rs √ κ b κ c satisfies 0 ≤ k ≤ √ κ b κ c . We remark that the vacuum input fields b in (t) andĉ in (t) are not independent, as can be seen from Eq. (8).
In the following section we will apply this scheme to an open version of the Dicke model, and find optimal delay-times for parameters motivated by recent experiment.

Application to the open Dicke model dynamics
We will first introduce the generalized Dicke model that will be the main object of our study. It describes the interaction of N two-level atoms with a single mode of the electro-magnetic field. The Hamiltonian describing the internal dynamics of the system isĤ where parameters ω 0 and ω are atomic and cavity frequencies, respectively, g is the linear interaction strength, and U is a non-linear coupling constant. The operatorŝ a andâ † are, again, the annihilation and creation operators for the cavity mode, and {Ĵ x ,Ĵ y ,Ĵ z } are collective atomic operators satisfying the conventional angular momentum commutation relations. The Hamiltonian in Eq. (10) is identical to the conventional Dicke model [21] Hamiltonian if U is set to zero. We will not enter into the details on the underlying physics of the BEC experiments realizing Eq. (10), but refer the reader to Ref. [9]. Briefly, the two-level atoms of the Dicke model are realized through pairs of discrete momentum states of a BEC. The linear coupling to the cavity field, (2g/ √ N )Ĵ x (â † +â), effectively describes Rayleigh scattering of photons between the cavity mode and an auxiliary pump laser. Importantly, the effective coupling strength, g ∼ √ P , can be tuned via the laser intensity P . The parameter ω is determined by the detuning of the cavity mode frequency from the pump laser frequency, and is therefore controllable. The collective atomic frequency ω 0 is fixed by the optical wave-vector and atomic mass (i.e., it is set by the recoil energy). It is therefore not readily tunable like the other parameters. The non-linear coupling constant U is given by a dispersive light-shift. We will consider parameters where this non-linear coupling does not play a major role, but refer to Refs. [10,11,30] for discussions on the very interesting dynamics that can result from this term.
We will assume that the cavity consists of two mirrors, both with decay rates κ b = κ c ≡ κ/2 (for simplicity), and implement the feedback scheme introduced in the previous section and illustrated in Fig. 1. We will, furthermore, assume the "thermodynamic" limit N → ∞, where we can employ a semi-classical approach to find expectation values of system operators, as described below. Later in Section 5, we examine quantum fluctuations by linearizing the operator equations of motion around these expectation values.
By making use of Eq. (9) and Eq. (10), we can now write down the Heisenberg equations of motion for the cavity mode and the spin operators: where now 0 ≤ k ≤ κ/2, andâ in (t) = 1/ √ 2 b in (t) +ĉ in (t) is used to denote the sum of the vacuum input fields through the two mirrors and where, for clarity, the time-dependence has been made explicit. Sinceb in (t), which corresponds to a vacuum part of the field incident on mirror b, is not independent of the input fieldĉ in (t) on mirror c,â in (t) does not obey the usual relations Eq. (2) and Eq. (3), but instead we have that which can be verified by making use of Eq. (8).
The non-linear operator equations, Eq. (11), can not be solved directly, and with delayed feedback (τ = 0), no equivalent master equation can be derived. By neglecting fluctuations and factorizing operator products, we can, however, derive a closed set of equations of motion for the five real-valued variables: which take the form Here, again, the time-dependence has been made explicit for clarity. These equations of motion conserve the total length of the spin j 2 x + j 2 y + j 2 z , and we will restrict our study to states on the Bloch sphere and thus always assume the constraint Below we will also make use the rescaled complex variable α ≡ â / √ N = x 1 + ix 2 , and, for notational convenience, we also write x ≡ (x 1 , x 2 , j x , j y , j z ).
In the absence of feedback, i.e. k = 0, this model has been explored theoretically in great detail, both in the thermodynamic limit N → ∞ in Refs. [10,11], and for finite N and including all quantum effects in Ref. [30]. In Refs. [10,11], steady states were found analytically by setting the left-hand sides in Eq. (14) equal to zero, and a further stability analysis was performed by linearizing around the fixed points. A surprisingly rich phase diagram was uncovered, with the appearance of several new phases due to the presence of the non-linear coupling U and the cavity decay parameter κ when compared to the conventional Dicke model [31,32]. One of the key findings in Ref. [11] was that the emergent time-scales of the collective dynamics vary significantly throughout the phase diagram, and that the run-times of current experiments may not be sufficient to reach the long time attractor in all cases.
Here we will now focus our attention on the following set of parameters as taken from the recent experiment in Ref. [23]: With this choice of parameters, a phase transition from the normal phase, x ⇓ ≡ (0, 0, 0, 0, −1/2), to a super-radiant phase with {x 1 , x 2 , j x } = 0, happens at a critical coupling strength [11] Below this critical coupling there are two fixed points: the normal phase and the inverted phase where only the former is stable in the absence of feedback (k = 0). Above the critical coupling, both of these phases are unstable in the absence of feedback, and two new stable fixed points, x SR ± , come into existence, given by [11] The two fixed points x SR ± differ only in the choice of sign for j x and α, related to a duality emerging from the invariance of Eq. (10) under the parity transformation a → −â,Ĵ x → −Ĵ x . We can treat the two solutions simultaneously, and we will refer to both as "the super-radiant phase", and denote them both by x SR , when the difference between them is of no importance. The super-radiant phase transition (in the absence of feedback) is qualitatively illustrated in Fig. 2.
The treatment in Refs. [10,11] showed that the approach to steady state can be exceedingly slow. This is partly related to the relatively small value of ω 0 , and can in fact be interpreted as critical slowing down, due to closeness to phase boundaries in an extended parameter space [10,11]. By careful adiabatic elimination of the cavity mode, under the condition {ω, κ} ≫ ω 0 , one can find an effective rate describing the incoherent dynamics of the atoms. In the normal phase, for example, it is found to be ≃ 4κg 2 (ω − U/2)ω 0 / (ω − U/2) 2 + κ 2 2 [11] (see also equation Eq. (21b) below). At g = g c this is roughly 0.3·2π Hz, which indicates a remarkably slow decay taking on the order of seconds. The parameter ω 0 , describing the collective atomic frequency, is fixed by the optical wave vector and the atomic mass, and can therefore not easily be made larger in practice. In Fig. 3, we show two typical examples of the atomic inversion,  figure), as functions of g, in the absence of feedback (k = 0). Below critical coupling, g < gc, there are two fixed points x ⇑ (green) and x ⇓ (blue), where only the latter is stable. Above critical coupling, g > gc, they are both unstable, and two new stable fixed points, x SR ± (red), come into existence, satisfying Eq. (20). Note that the positions of fixed points stay the same under TDAS feedback control, only their stability might change.
j z (t), as a function of time, below critical coupling, with g/g c = 0.74 and above, with g/g c = 1.1. The other parameters are as in Eq. (16). In both cases the initial state is taken to be x = (0, 0, 1/ √ 12, 1/ √ 12, 1/ √ 12). In Fig. 4, we similarly show the time-evolution of the normalized photon number |α(t)| 2 for the same initial state and parameters. These figures clearly show the exceedingly slow approach towards the stable steady state, in agreement with the predictions from Ref. [11].
In order to obtain the characteristic time-scales governing the approach to a fixed point,x(t) ≡x, when feedback is applied, we consider a small perturbation, y(t) = x(t) −x, and linearize the equations of motion around the solution. By using an Ansatz y(t) = exp(λt)y 0 , we can then derive a characteristic equation for λ. Each solution λ corresponds to a characteristic (inverse) time-scale for the dynamics close to the steady state. In the presence of feedback, k, τ > 0, there are in fact an infinite number of solutions λ k . However, a crucial result in the analysis of delay differential equations is that there are only a finite number in any real half-plane Re λ > σ, σ ∈ R [33]. Thus it becomes feasible to find the slowest λ k that ultimately governs the time-scale for approaching or leaving a fixed point. The details of a stability analysis for Eq. (14) are given in Appendix A.
A steady state solutionx is stable only if all λ k have negative real parts, i.e., Re λ k < 0. The time-scale for the approach to a stable steady state solution is thus governed by the eigenvalue with real part closest to zero, which we will denote by λ 1 .
The key to our control scheme is that the eigenvalues can be manipulated through the variation of k and τ . In particular, both the magnitude and sign of Re λ 1 can be changed. This opens up the possibility of changing the emergent time-scales, as well as qualitatively changing the phase diagram of the system, by changing the stability of a steady state.
We can find λ 1 numerically by solving a transcendental characteristic equation, given in Eq. (A.5). In Fig. 5, we plot Re λ 1 as a function of τ for k = κ/2, g/g c = 0.74 and 1.1, while the other parameters are kept as in Eq. (16). For g/g c = 0.74, we linearize around the normal phase, x ⇓ Eq. (18), and the inverted phase x ⇑ Eq. (19),   18), is the only stable fixed point in the absence of feedback. The blue line shows the time-evolution without feedback. The red line is with k = κ/2 and τ = 50 µs. The green line shows the time-evolution for k = κ/2 and τ = 100 µs, for which the normal phase is unstable, and the system exhibits persistent oscillations in the long time limit. Bottom panels: g/gc = 1.1; the super-radiant phase, x SR as in Eq. (20), is the only stable fixed point in the absence of feedback. The blue line shows the time-evolution without feedback, and the red line is for k = κ/2 and τ = 50 µs. For both panels, all other parameters are given in Eq. (16).
which are always valid fixed points. For g/g c = 1.1 we also linearize around the superradiant fixed point, x SR Eq. (20). This analysis can be used to find optimal values for τ close to the steady state. We observe that a minimum value for Re λ 1 is reached for a delay of around τ ≃ 50 µs, for both the normal phase when g/g c = 0.74 and the super-radiant phase when g/g c = 1.1. The value of Re λ 1 is, without feedback (τ = 0), roughly −0.14 · 2π Hz and −0.35 · 2π Hz, respectively, for the two cases. In comparison, at the first minima, with τ = 52 µs and τ = 50 µs, the values are −26 · 2π and −58 · 2π Hz, i.e., two orders of magnitude larger.
In Fig.  3 and Fig.  4 the time-evolution of the collective inversion, j z (t), and the normalized photon number, |α(t)| 2 , from an initial state x =  The green line shows the time-evolution for k = κ/2 and τ = 100 µs, for which the normal phase is unstable, and the system exhibits persistent oscillations in the long time limit. Bottom panels: g/gc = 1.1; the super-radiant phase, x SR Eq. (20), is the only stable fixed point in the absence of feedback. The blue line shows the time-evolution without feedback, and the red line is for k = κ/2 and τ = 50 µs. For both panels, all other parameters are given in Eq. (16).
(0, 0, 1/ √ 12, 1/ √ 12, 1/ √ 12) T , is shown for various values of τ , with k set to κ/2, and either g/g c = 0.74 or g/g c = 1.1. The other parameters are again as given by Eq. (16). With τ = 50 µs the steady state is reached after ∼ 20 ms, a dramatic improvement when compared to the dynamics without feedback, for which the relaxation takes several seconds. These figures also show the possibility of qualitatively changing the phase diagram by choosing a value of τ that de-stabilizes a fixed point. Note that these results were found by numerically integrating Eq. (14), using an integrator designed for delay differential equations [34].
We note that significant improvement can also be achieved for smaller delay-times than those used in Fig. 3 and Fig. 4. We find an approximate expression for λ 1 , valid Bottom panel: g/gc = 1.1. The green line is for the inverted phase, and the red line is for the super-radiant phase, x SR Eq. (20). The normal phase is not shown for this value of g, but has a value of Re λ 1 ≃ 3.5 · 2π kHz throughout (and is thus unstable). The dashed lines show the corresponding linear approximations given in (21a,21b). For both panels, k is set to κ/2, and all other parameters are given in Eq. (16).
So far, we have not accounted for any loss in the feedback loop, and considered only the ideal case k = κ/2 in our numerical results. However, good results are also achieved for smaller k, as already indicated by the approximate linear dependence in Eq. (21b). In Fig. 6, we show Re λ 1 as a function of k and τ , for linearization around two steady states: the normal phase, for g/g c = 0.74, and the super-radiant phase, for g/g c = 1.1. The other parameters are as before Eq. (16). This shows that great improvements in the convergence are also possible with significant loss in the feedback loop. For the choice k = 0.1·κ/2 (90% loss) we, e.g., still expect an order of magnitude improvement in the relaxation time for the parameters considered.

Finite-time experiments and unexplored regions of the phase diagram
We have illustrated the possibility of very slow time-scales associated with the generalized Dicke model when realistic parameters are used, by considering the timeevolution of an initial state, chosen far from equilibrium ( Fig. 3 and Fig. 4). It is, however, important to note that in a typical experiment [9,22,23] the system is prepared in the normal phase for some g < g c , before gradually ramping up g to a value beyond g c . Due to the continuous nature of the phase transition at g c , the change in steady state is gradual enough that the system can react to the small rate of change.
As the experiments have shown, as well as the theoretical modelling in Ref. [11], the measured photon intensity agrees well with the semi-classical steady state value.
We will here follow the approach in Ref. [11] to emulate a finite-time experiment. We will start with the system prepared close to the normal phase, and then ramp up g(t) ∼ √ t to a value beyond the critical point. Specifically, we will choose g(t) = t/t 0 · 1.5g c . To emulate the effect of quantum fluctuations, we prepare initial state of the system close to the normal phase, with j x = j y = 1/ √ N , and N = 10 5 . We now envision an experiment where the goal is to reach the steady state well beyond the critical point, at g/g c = 1.5. For this purpose, we perform a linearization around this steady state, and find a good value for τ , as was illiustrated in Fig. 5. We use, again, k = κ/2, and the other parameters as in Eq. (16). We then find a minimum for Re λ 1 around τ = 16 µs. In Fig. 7, we show the time-evolution of the collective inversion and the photon number as g(t) is ramped up, and compare the situation with and without feedback. Furthermore, we compare two sweeps, one with t 0 = 20 ms, and one with t 0 = 200 ms. We observe that the system, both with and without feedback, responds to the change across the critical point, and closely follows the adiabatic evolution according to the exact steady state values. We, however, note that in the case with feedback, the fluctuations are much smaller. A fully quantum treatment of fluctuations will be given in the next section.
One of the main findings in Ref. [11] was, however, that there are regions of the phase diagram where the situation is far more problematic. This was particularly found to be the case for negative ω (ω < U/2 < 0), where emulation of finite time experiments with sweeps up to 200 ms were not able to approximate the superradiant steady state. This region of parameter space has also not been explored experimentally. Here the normal phase is unstable below threshold, but the time-scale for leaving the normal phase is extremely slow. If the system is prepared close to the normal phase, it is therefore essentially meta-stable on typical experimental runtimes. However, it will be far from steady state when the experiment hits threshold, g = g c , and the system is unable to "respond" to the super-radiant phase transition on an adequate time-scale. In Fig. 8, we exhibit finite-time sweeps with t 0 = 20 ms and 200 ms, but now for ω = −10.0 · 2π MHz. The other parameters are kept as in Eq. (16), and we compare the case with no feedback (k = 0) to the case with (k, τ ) = (κ/2, 16 µs) as before. We see that, without feedback, the phase transition is not visible due to the slow response of the system. In other words, the time-evolution is far from being adiabatic. The results with feedback, on the other hand, are in this case quite interesting. At this value of τ , the system responds rapidly, and does a relatively fast switch from the normal to the inverted phase, before hitting threshold. The evolution subsequently adapts well to the super-radiant phase transition across g c , and for both sweep times the target steady state at g/g c = 1.5 is reached to a very good approximation.

Quantum Fluctuations
So far, we have considered the semi-classical amplitudes of the relevant observables. It is of importance to consider how quantum fluctuations are influenced by our feedback control scheme. We consider the thermodynamic N → ∞ limit, and follow the treatment in Refs. [31,32,35] by introducing a Holstein-Primakoff representation for the collective spin, i.e.
whereĴ ± =Ĵ x ± iĴ y , andb (b † ) are bosonic annihilation (creation) operators satisfying [b,b † ] = 1. We next introduce fluctuation operators by expanding the fields around their semi-classical amplitudes,â ≡ â + δâ,b ≡ b + δb, where we use b = ± N/2 + Ĵ z . For the normal phase we have â = b = 0, whereas  in the super-radiant phase we have â / √ N = α SR , b / √ N = ± 1/2 + j SR z , where α SR , j SR z are given in Eq. (20). We observe that a positive choice for b corresponds to a negative choice for â , and vice versa. In the following, we will just consider the positive choice for b , as the calculation with the other choice is essentially identical, with a few changes in sign.
After a lengthy calculation, one finds that the fluctuations satisfy the following equations of motion, in the limit N → ∞ (see Appendix B): and and where These expressions are valid in both the normal and the super-radiant phase, whereᾱ andj z refer to the semi-classical amplitudes in the respective phases. The set of euqations (23a-23b) are linear, and as was done in Ref. [35], we solve these equations by making use of Fourier-transform techniques as explained in more detail in Appendix C. The mean fluctuations in the intra-cavity photon number in steady state can then be computed through where δã(ν) is the Fourier transformed cavity field fluctuation. We solve this integral numerically, while varying g, to investigate the diverging fluctuations upon approaching the critical value g c . In Fig. 9, we compare the case with no feedback (k = 0) to the case with (k, τ ) = (κ/2, 50 µs). We see that the feedback can significantly influence the size of the quantum fluctuations. The change of about two orders of magnitude is consistent with the results from the previous sections. Interestingly, we observe that a series of new instabilities develop in the super-radiant phase for increasing g when feedback is applied. We also note that the scaling of (g c /g) δâ † δâ ss with |1 − g/g c |, as g → g c , is found to be the same in the two panels of Fig. 9, corresponding to a universal "photon flux exponent" of 1.0 [23,36]. Fluctuations in theb-mode, δb † δb ss , show similar divergences at the same values of g.

Conclusions and outlook
We have modelled Pyragas' time-delay auto-synchronization feedback control applied to a quantum system. For a topical cavity QED many-body system, the Dicke model as realized in recent experiments, we have investigated in what manner a feedback force can influence the characteristic time-scales governing the relaxation of the system as well as the characteristic size of quantum fluctuations. With optimal feedback times, the relaxation-time can be reduced by two orders of magnitude, with a corresponding decrease in quantum fluctuations. Even with significant loss in the feedback loop (90% loss), one expects an order of magnitude improvement. The scheme put forward also offers the possibility of changing the stability of long time attractors, thus qualitatively changing the phase diagram of the system. A fixed point might, e.g., be de-stabilized at a critical value of the delay-time τ and thus inducing a novel feedback-driven phase transition.
Although we have focused specifically on the Dicke model in this paper, and particularly on how to reduce the relaxation time, we believe the scheme could be applied to a variety of topical CQED systems. It might, e.g., be interesting to consider systems with optical bistability, where the scheme may be used to change the stability of the fixed points. In other systems, it may also be of interest to consider how the feedback control can slow down the dynamics, instead of making the convergence more rapid, if the goal is to preserve quantum information in an initial state. In general, the scheme presented in the present paper offers a novel and non-invasive way to control the overall time-scale governing the dynamics of open quantum systems, where noninvasive here means that the positions of fixed points in parameters space remain unchanged, only their relative stability is changed. Here we summarize the semi-classical stability analysis, leading to the characteristic time-scale governing the time-evolution close to a fixed point. For notational convenience, we will write x ≡ (x 1 , x 2 , j x , j y , j z ) and define a real vector-valued function f such that Eq. (14) can be written in the form where the matrix B ≡ diag(1, 1, 0, 0, 0). We now consider a small perturbation, y(t) ≡ x(t) −x(t), from a solutionx(t) to Eq. (A.1), and linearize the equations of motion, i.e., Before proceeding, it is convenient to eliminate one of the spin variables from Eq. (A.2) by using the constraint Eq. (15). By differentiating this constraint, and using the Ansatz y(t) = exp(λt)y 0 , we eliminate j z , leading to a linear set of equations for where B ′ ≡ diag(1, 1, 0, 0) and the matrix A ′ reads Here we have introduced the notationω ≡ ω + Uj z andω 0 ≡ ω 0 + U (x 2 1 +x 2 2 ) ≡ ω 0 + U |ᾱ| 2 , whereᾱ ≡x 1 + ix 2 . Each solution to the characteristic equation for λ, corresponds to a characteristic (inverse) time-scale for the dynamics close to the steady state. Here I 4 is the 4 × 4 identity matrix. For all the fixed points we are interested in, we have j y = 0, which simplifies Eq. (A.5). Explicitly, we find for j y = 0 that The transcendental equation (A.7) has an infinite number of roots when k, τ > 0. However, a crucial result in the analysis of delay differential equations is that there are only a finite number of roots in any real half plane, Re λ > σ, σ ∈ R (see e.g. Ref. [33]). Thus it becomes possible to find the smallest root, which we denote by λ 1 , that ultimately governs the time-scale for approaching or leaving a fixed point. In Fig. 5 and Fig. 6, λ 1 was found numerically for varying τ and k and different fixed pointsx. This was done using a MATLAB tool for analyzing linear delay differential equations [37].

Appendix B. Quantum fluctuations in the thermodynamic limit
Here we outline in some more detail the calculations leading to Eqs. (23a), (23b) and (24). After introducing the Holstein-Primakoff representation, i.e., and neglecting constant energy terms, we have that the Dicke Hamiltonian Eq. (10) can be rewritten in the form Next, we expand the fields asâ = a 0 + δâ,b = b 0 + δb, where a 0 ∈ C, b 0 ∈ R, and write:Ĥ Let us first consider the normal phase, where we set a 0 = b 0 = 0, before considering the more involved super-radiant phase.
Appendix B.1. The normal phase In this case, we have that Eq. (B.4) becomeŝ Next, we make use of the approximation and by inserting this into the Hamiltonian above, and neglecting the term proportional to U/N that vanish in the limit N → ∞, we obtain an effective Hamiltonian for the normal phase:Ĥ Comparing this with Eq. (24), we see that they agree upon insertingj z = −1/2,ᾱ = 0.

Appendix B.2. The super-radiant phase
For the super-radiant phase, we have that a 0 , b 0 behave like O( √ N ). We therefore expandξ to order O(1/ √ N ), i.e.
where k ≡ N − b 2 0 . Neglecting constant terms, and terms that vanish as N → ∞, we then arrive at the following effective Hamiltonian: − gb 0 (a 0 + a * 0 ) 4k 2 and that Hence, the terms in H ′ that are linear in δâ ( †) , δb ( †) will vanish in the Heisenberg equations of motion. Thus, after inserting the expressions for a 0 , b 0 given above, we find that we can use Eqs.(23a)-(23b) with the Hamiltonian as given in Eq. (24).

Appendix C. Time-delayed feedback in linear quantum systems
Many systems of interest in quantum optics have linear Heisenberg equations of motion, and they are typically treated by introducing Fourier transformed fields [38]. We extended this treatment to linear systems with feedback, as described by the following general form: where a = (â 1 ,â † 1 , . . . ,â n ,â † n ) is a vector of field modes and their adjoints, A is a matrix coupling the different fields, Γ = diag(κ 1 , κ 1 , . . . , κ n , κ n ) is a diagonal matrix of decay rates, similarly √ 2Γ = diag( √ 2κ 1 , √ 2κ 1 , . . . , √ 2κ n , √ 2κ n ) gives the coupling to the input fieldsâ in = (â 1,in ,â † 1,in , . . . ,â n,in ,â † n,in ), and the matrix K i couples the system to the feedback fieldsâ(t − τ i ).
Next, we introduce the Fourier transforms for any operatorÔ. Since 1/ √ 2π The system fields are thus solved in terms of the input fields by inverting the matrix on the left hand side.