Multimode cold-damping optomechanics with delayed feedback

We investigate the role of time delay in cold-damping optomechanics with multiple mechanical resonances. For instantaneous electronic response, it was recently shown in \textit{Phys. Rev. Lett. \textbf{123}, 203605 (2019)}, that a single feedback loop is sufficient to simultaneously remove thermal noise from many mechanical modes. While the intrinsic delayed response of the electronics can induce single mode and mutual heating between adjacent modes, we propose to counteract such detrimental effects by introducing an additional time delay to the feedback loop. For lossy cavities and broadband feedback, we derive analytical results for the final occupancies of the mechanical modes within the formalism of quantum Langevin equations. For modes that are frequency degenerate collective effects dominate, mimicking behavior similar to Dicke super- and subradiance. These analytical results, corroborated with numerical simulations of both transient and steady state dynamics, allow to find suitable conditions and strategies for efficient single or multimode feedback optomechanics.


I. INTRODUCTION
A widespread technique for the removal of thermal noise from a given mechanical degree of freedom involves electronic feedback loops. The procedure is based on the continuous monitoring of a system's observable, followed by the application of an adequate cooling action via the feedback device. For example, in optomechanics [1-11], a cavity field quadrature is detected and the result is applied either optically (as a radiation pressure force) or electrically to the thermally activated mechanical resonator. While generally one aims for the isolation and cooling of a specific vibrational mode, it has been recently shown that efficient simultaneous cooling of a few independent modes is also possible either in the case of using sideband cooling [12], via machine learning [13] or cold damping [14].
We provide here a more in-depth analytical treatment of simultaneous cold-damping of many mechanical resonances [15][16][17] and address a crucial aspect extremely relevant in experiments, i.e. the inherent time delay τ inh that characterizes any electronic feedback loop. It is generally agreed that the delayed action of the feedback loop can lead to unwanted heating eventually leading to an instability [18][19][20][21]. Extending the analytical approach that we have previously introduced in Ref. [14] to include a delay time τ , we show that analytical solutions are still possible to some degree in the fast-feedback-lossy-cavity (FFLC) regime. In this regime, lossy cavities allow for quick read-out and broadband feedback allows for quick cooling. Most importantly, we suggest that in order to counteract unwanted heating effects, the feedback loop delay time could be further delayed by introducing an additional delay time τ add . By fitting the total delay τ = τ inh + τ add to the characteristics of the system, cooling efficiency close to the level of the τ = 0 can in some cases be achieved.
Simultaneous cooling of N mechanical resonances can provide either a wider bandwidth for sensing applications or a stronger optomechanical coupling to a collective mode. We show that a single feedback loop can very efficiently couple to a bright collective mode which, in the near degeneracy case where all mechanical modes lie within a very narrow frequency window, can be up to N times faster damped to a N times lower occupancy than a single mode. This is reminiscent of the superradiance effect as in an increase in the collective radiative rate for a system of N quantum emitters coupled to a single bosonic mode, as in the Dicke model in quantum optics [22,23]. The corresponding effect of subradiance, i.e. strong suppression of radiative rate, is mimicked by the decoupling of the other N − 1 collective dark states from the feedback loop. For efficient cooling of many resonances in a wide frequency window, one then has to instead engineer a linear dispersion relation such arXiv:2006.08430v1 [quant-ph] 15 Jun 2020 that all adjacent modes are separated by more than the damping rate introduced by the feedback loop. In terms of collective modes, the spread of many mechanical resonances over a large frequency interval insures strong bright-dark couplings which in turn leads to sympathetic cooling of all dark modes.
The analytical treatment followed here is based on solving a system of coupled quantum Langevin equations for N mechanical modes coupled to a single cavity optical mode and subjected to a feedback force. In Sec. II we introduce the model that includes the feedback force with time delay and detail the procedure that allows for the linearization of the radiation pressure interaction in the high-amplitude field limit. In Sec. III we derive the simplified equations of motion for N coupled mechanical degrees of freedom based on the FFLC approximation. The dynamics is generally non-Markovian as the task is to solve a system of coupled integro-differential equations; however, for relatively small time delay, we introduce the weak (wFFLC) and the strong (sFFLC) Markovian approximations which simplify the task by turning the dynamics into a set of coupled linear differential equations. In Sec. IV we provide analytical results for the steady state of the system by using either a time domain analysis, particularly useful under the Markovian approximation but also in the Fourier domain with a generality extending into the non-Markovian regime as well. We benchmark important results in Sec. V for a single resonance, highlighting the role of time delay and showing that an additional delay in the feedback signal with a conveniently chosen τ add can improve cooling efficiency. In Sec. VI we extend these results to the two modes case, elucidating the interplay between collective damping and the time delay effects. The different levels of approximations are then tested against exact numerical simulations based on solving time dynamics of a set of stochastic differential equations. Finally, in Sec. VII we present analytical and numerical results for many resonances in particular following a linear dispersion relation and in the case of degenerate modes and illustrate strategies for efficient cooling with adjustable feedback time delay.

II. MODEL
We follow the evolution of an optomechanical system at the level of operators subject to both unitary evolution as well as to dissipation (included as optical and thermal quantum fluctuation input noises, i.e. the standard quantum Langevin approach in optomechanics [24]). The system is comprised of an optical cavity mode coupled via the radiation pressure Hamiltonian to N mechanical resonances of a single vibrating end-mirror. The N independent modes of vibrations have effective mass m j and frequency ω j . The quantum Langevin equations of motion [24] for the N + 1 degrees of freedom reaḋ We have introduced dimensionless position and momentum quadratures Q j and P j for each of the N independent membrane oscillation modes with standard commutations [Q j , P j ] = iδ jj . The term ∆ 0 = ω c − ω describes the detuning of the cavity resonance frequency ω c from the laser frequency ω and κ its decay rate. The input laser power is given by = 2Pκ/ ω . The optomechanical coupling is described by the radiation pressure Hamil- OM is the single-photonsingle-phonon coupling rate for the j-th mode. The single cavity mode at frequency ω and loss rate κ is described by the bosonic operator A with [A, A † ] = 1. The zeroaverage noise terms are delta-correlated in the time domain a in (t)a in † (t ) = δ(t − t ).
The parameter γ j describes the damping of the j-th resonator mode and together with the associated zeroaveraged Gaussian stochastic noise term ξ j fulfill the fluctuation-dissipation relation resulting in thermalization with the environment. The noise term can be fully described by the two-time correlation function: where Ω is the frequency cutoff of the reservoir and S th (ω) = ω[coth ( ω/2k B T ) + 1] is the thermal noise spectrum. A standard white noise input with delta correlations both in frequency and time is obtained for sufficiently high temperatures k B T ω j from the correlation function resulting in the approximate form ξ j (t)ξ j (t ) ≈ (2n j + 1)γ j δ(t − t )δ jj , wherē n j = (exp( ω j /k B T ) − 1) −1 ≈ k B T / ω j describes the average occupancy of each vibrational mode.

A. Linearization
Let us rewrite all operators A = A +a, Q j = Q j +q j and P j = P j + p j as a sum of their expectation value and zero-averaged fluctuations. When the cavity field amplitude is large with respect to the fluctuations, one can simplify the equations of motion by neglecting terms such as a † a as being small compared to | A | 2 . Under this approximation, the classical averages satisfy the following equations of motion In steady state (obtained by setting Q j = Ṗ j = Ȧ = 0) the cavity field amplitude can be shown to satisfy the following non-linear equation The non-linearity stems from the intensity dependent cav- | A | 2 /ω j owing to the radiation pressure induced displacement from equilibrium. The steady state value for the displacements are With the omission of the small nonlinear terms a † a and aq j , we obtain the linearized equations of motioṅ where x = (1/ √ 2)(a + a † ) and y = (i/ √ 2)(a † − a) are the fluctuations of the quadratures of the cavity field and x in and y in are similarly defined in terms of field noise operators. The effective optomechanical coupling terms are given by G j = √ 2g (j) OM A and are enhanced by the large cavity field amplitude. We set the condition that the effective cavity detuning OM Q j , containing a collective mechanically-induced frequency shift is kept at zero value.

B. Time-delayed feedback
The application of the feedback requires the readout of a cavity field quadrature followed by the appropriate action onto the mechanical resonator. Generally, one can express the applied force as where the convolution term is defined as (g j (t − s)y(s) and depends on the past of the detected quadrature y that is driven by the weighted sum of the oscillator fluctuations q j . The causal kernel g (τ ) contains the feedback gain terms g (j) cd and feedback bandwidth ω fb . The Fourier transform of the feedback kernel is given by which resembles a standard derivative high-pass filter which here additionally contains a delay dependent phase term (similar to the term expressed in [20]) in contrast to the previous case [14,24]. In addition, the parameter τ (neglected previously in Ref. [14]) is the joined feedback delay originating from the measurement signal processing. Within the convolution with θ(t − τ − s) this parameter guarantees that only information up to t − τ can influence the dynamics of the resonator modes. Notice that in the limit ω fb → ∞ the feedback becomes g The quadrature component that is injected into the feedback mechanism y est is the estimated intra-cavity phase quadrature given by which results from a measurement of the output quadrature y out = √ 2κy(t) − y in (t). This follows from the description of a detector with quantum efficiency η which is modeled by an ideal detector preceded by a beam splitter with transmissivity √ η, which mixes the input field with an uncorrelated vacuum field y v (t). Finally, to fully describe the dynamics of an optomechanical system with N mechanical resonances undergoing cold damping with time delay one corrects Eqs. (5) with the following equatioṅ In the next two sections we will describe strategies that can be employed to simplify the equations of motion and deliver Markovian approximations accurately describing the time dynamics of the system for small time delays. The approximations also allow for the derivation of analytical estimates of final occupancies for all modes undergoing cold damping. An extension beyond the Markovian regime will then be obtained by a Fourier analysis of the coupled system of equations in steady state.

III. MULTI-MODE COLD DAMPING: SIMPLIFIED EQUATIONS
The aim of this section is to arrive at a set of 2N coupled equations describing solely the dynamics of the N mechanical resonator modes. To this end we proceed by formally integrating the equations of motion for the optical degree of freedom and replacing them in the equations for the mechanical modes. We first find a general formulation for the integro-differential non-Markovian collective dynamics where feedback-induced damping occurs as a time convolution involving momentum quadratures evaluated at times in the past. Under the FFLC approximation, Markovian collective dissipative dynamics emerges, allowing one to write a linear set of coupled differential equations analytically solvable in steady state.
A. Non-Markovian collective damping Let us start by formally integrating the dynamics of the optical degree of freedom We will aim at computing the estimated quadrature y est (s) which introduces both terms proportional to q j as well as noise terms stemming from the cavity input noise y in (s) and from the vacuum filled port noise y v . This will then give rise to the convoluted force acting on all mechanical momentum quadratures. In a first step we estimate the contribution coming from the intracavity field y as Notice that the feedback force above contains terms proportional to the mechanical displacement quadratures in addition to an extra feedback induced noise. It is however desired to express the effect of feedback as a damping force proportional to a momentum quadrature. To this end we apply integration by parts, noticing that the convolution contains a derivative of the following function h τ (t − s) = e −κ(t−s−τ ) − e −ω fb (t−s−τ ) /(ω fb − κ). In addition, we make use of the relationq j = ω j p j to obtain The feedback force explicitly shows a damping term exhibiting a non-local kernel: this indicates that the action depends on the past behavior of the momentum quadra-tures on timescales defined by the cavity loss rate κ and on the feedback bandwidth ω fb . We can list the set of coupled 2N equations for all resonator modes quadratureṡ The Markovian damping rate γ for each mode is supplemented with a diagonal non-Markovian feedback damping kernel as well as with off-diagonal dissipative couplings to all other modes. The three sources of noise, in addition to the thermal one ξ j , stem from the direct feedback action ξ fb , from the feedback filtered vacuum action in the loss port ξ vac and from the intra-cavity radiation pressure effect ξ rp . They are expressed as The explicit forms of the convolution kernels φ

B. Markovian collective damping
The set of integro-differential equations obtained above is not easily tractable; however, in regimes favorable to cold-damping, a transformation to a much simpler form can be achieved. Let us assume a lossy cavity and relatively fast feedback such that both rates fulfill κ, ω fb ω j , for all j ∈ {1, . . . , n}. In such a case, integration by parts of the non-Markovian kernel can be performed and one can show that, in leading order, the convolution gives rise to a very simple expression We denote this approximation as the strong fast-feedback lossy-cavity (sFFLC) assumption. Under less stringent conditions with κ, ω fb > ω j (which we will refer to as the weak wFFLC condition) a slightly more complicated expression with a larger validity region can be derived and is detailed in Appendix A, E. Using Eq. 16 we derive the non-Markovian equations of motioṅ that depend on past events at time t − τ (higher order delay terms up to p k (t − 2τ ) are derived in Appendix E and indicate a dependence following p k (t − nτ ) for even higher orders). Here, Γ jk = g (j) cd G k ω k /κ is the feedback induced damping rate in the sFFLC approximation (the analogue equations to Eq. (17) in the case of the wFFLC are given in Appendix A). Notice that all noise terms have been gathered into a single term ζ j = ξ j + ξ fb + ξ vac + ξ rp . For τ = 0 the above expression can be plugged back in to give rise to a set of 2N coupled equations where all momentum quadratures are evaluated at the same time t (as detailed in Ref. [14]).
For τ > 0, however, the equations are more complicated to solve as momenta at time t are coupled to momenta in the past at t − τ . The next important simplification that we use implies that during the time τ the motion stays periodic with period ω j . Then one can roughly approx- The approximation is valid as long as damping of oscillations within the interval τ can be neglected. In the sFFLC, the equations of motion become then a simple set of coupled linear quantum Langevin equationṡ Again, for τ = 0 the result reduces to that previously derived in Ref. [14]. Notice that the main effect of non-zero time delay τ > 0 is to modify both the individual mode and mutual damping rates by the cosine factors cos(ω j τ ). This also means that for given time delays the system can exhibit instabilities when the feedback-induced heating rate surpasses the natural decay rate γ. The extra frequency renormalization terms proportional to sin(ω j τ ) are negligible as long as ω j Γ jk , a regime to which we will restrict ourselves in the following.
A particularly interesting regime is that of full frequency degeneracy where ω j = ω. A simplified picture can be used in this case in terms of collective bright Q 1 = k α 1k q k and dark modes Q j = k α jk q k (See Fig. 2). Here, the coefficients from the bright mode α 1j = G j / k G 2 k can help to acquire the coefficients for the N − 1 dark modes via the Gram-Schmidt procedure which satisfies the condition j α lj α kj = δ lk . The bright mode dynamics is described bẏ while the other N − 1 orthogonal dark modes satisfy the following equations of motioṅ The dark state manifold dynamics can be further simplified by injecting the solution for the bright mode resulting inQ where the expression for the compound noise term Ξ j is detailed in Appendix A. The above dynamics shows that in the fully degenerate case the bright mode is damped at a high rate: for equal coupling this rate is directly proportional to N . The dark modes are instead mostly unaffected by the feedback loop (except via the input noise) and decay at the natural decay rates γ. This is reminiscent of the collective monitoring of a collection of quantum emitters by their electromagnetic environment (either free space or cavity) which leads to collective radiative effects known as super-and subradiance [22,23]. The role of the collective bath is played here by the common electronic feedback loop which provides simultaneous dissipative dynamics for all modes.

IV. MULTI-MODE COLD DAMPING: STEADY STATE
To derive the final achievable occupancies for all modes undergoing cold-damping, one can compute the covariance matrix of the system in steady state. We will follow two different paths: i) a time domain analysis suitable to the Markovian case, where the covariance matrix of the system can be computed from the Lyapunov equation and ii) a Fourier domain analysis which only applies to the steady state but presents the advantage of providing exact solutions in the Fourier domain even for the non-Markovian case.

A. Time domain analysis
The linearized quantum Langevin equations presented in Eq. (18) can be rewritten in compact vector forṁ with a vector of fluctuations v = (q 1 , p 1 , . . . q N , p N ) and the corresponding input noise vector n in = (0, ζ 1 , . . . , 0, ζ N ). Here, the elements of the matrix M are defined by the coefficients for q j and p j in Eq. (18). Under the condition that the system is stable, i.e. if all eigenvalues of M have negative real parts, one can find the covariance matrix This can be greatly simplified as all two-time correlations of noise terms n in (s)n in (s ) are delta-like in the strong fast-feedback-lossy-cavity regime with ω fb , κ ω j . One can show that and zero otherwise. These terms can be gathered into a diffusion matrix with elements D in,2i,2j = (2n i + 1)γ i δ ij + G i G j /κ and zero otherwise. The final step involves solving a Lyapunov equation for the covariance matrix Notice that the diffusion matrix shows no dependence on the time delay.
One can proceed to solve the Lyapunov equation by introducing the following notations for momentum correlations, X ij = p i p j + p j p i , position correlations Z ij = q i q j + q j q i and cross-terms Y ij = q i p j + p j q i . From the diagonal elements Z ii and X ii one can then estimate the final occupancy of each mode. With these notations, analytical results can be obtained by solving the following set of algebraic equations (simplified below under the sFFLC approximation) For the wFFLC analogue equations can be obtained in the case k B T ω j which are presented in Appendix A, where we can ignore the contribution from feedback, measurement and radiation pressure noise terms which do not show delta-like correlations in this regime. The form above is more complex than the τ = 0 case studied in Ref. [14] showing the influence of the delay in the weight factors c i = cos(ω i τ ) and s i = sin(ω i τ ) multiplying with the rates Γ ij . Setting τ = 0, quasi-exact but cumbersome expressions for the final occupancy of each mode can be obtained (as detailed in Appendix A and presented in [14]). The non-zero delay case however is more complicated and analytical expressions are harder to obtain except in the single mode and two adjacent modes case which we detail in the next sections.

B. Fourier domain analysis
The time-domain analysis provided above has a domain of validity restricted by the Markovian assumptions implied in the sFFLC and wFFLC approximations. However, in steady state, one can turn the integro-differential set of coupled equations into a simple set of algebraic equations by transforming to the Fourier domain. This allows for solutions inside the non-Markovian regime. Let us write the equations and proceed by eliminating the Fourier components of the field quadratures. One then obtains a set of N equations which allows for the derivation of each mode's response to the input noise The matrix χ χ χ describes the susceptibility matrix of the system. We have introduced a frequency-dependent effective resonance frequency ω 2 j,eff (Ω) = ω 2 j + Ωδω j (Ω) cos(Ωτ ) + ΩΓ j (Ω) sin(Ωτ ) and the corresponding frequency-dependent effective decay rate γ j,eff (Ω) = γ j +Γ j (Ω) cos(Ωτ ) − δω j (Ω) sin(Ωτ ). Additionally, we obtain the terms It is interesting to note that one can immediately obtain the wFFLC steady state computed in the time domain, under the approximation ω 2 j,eff (Ω) ≈ ω 2 j,eff (ω j ) and To obtain the steady state solution of the resonator mode occupation from the Fourier transform we need to calculate where S S S(Ω) describes the position-fluctuation spectrum as presented in the Appendix C, which for high temperatures k B T ω j can be approximated by S jk (Ω) ≈ γ j (2n j + 1)δ jk . In the following we will refer to the term S qj (Ω) = χ χ χ(Ω)S S S(Ω)χ χ χ † (Ω) jj as the position spectrum of the j-th mode.

V. SINGLE MODE COOLING
We provide here an analytical and numerical treatment of the time dynamics and steady state final occupancies for a single mode undergoing cold damping with a variable time delay. In steady state, the solutions to the Lyapunov equations under the weak or strong FFLC approximation provide simple, intuitive results for the final achievable occupancies. We provide a numerical validity check for the steady state solutions and extend the analytical calculations to regions of more general validity by solving the coupled set of non-Markovian equations in the Fourier domain.

A. Time domain analysis of steady state
For a single resonator mode, the solution to the Lyapunov equation leads to the following expressions for the momentum and position variances in steady state: We have made use of the optomechanical cooperativity defined here as C = G 2 /(κγ) [1]. This term brings an extra heating contribution owing to the back-action of the continuous monitoring of the field quadrature. As opposed to the cavity self-cooling scheme [1], where one aims at large cooperativities, here we aim to keep this term small. This is easily achieved by making the cavity lossy, i.e. κ G. Notice that, generally, as pointed out also previously for the single mode feedback cooling without time delay [24], the equipartition theorem does not hold and therefore the damped state is not in thermal equilibrium. However, we define an approximate final occupancy quantity that assumes the following analytical expression n eff (τ ) = 1 2 .
(32) In the limit of zero time delay, this reduces to the expected result n eff = γ(n + 1/2 + C/2)/(γ + Γ) (as in Ref. [24]), where the effective damping rate is γ + Γ. With non-zero time delay the rate γ + Γ cos(ωτ ) is reduced and eventually can become negative when the feedback acts completely out of phase and an instability can occur. If the delay is small such that ωτ 1, the effect is minimal as the damping rate is still almost optimal Γ cos(ωτ ) ≈ Γ. For values of ωτ close to π/2 first inefficient cooling and then an instability will occur as it is presented in Fig. 3b, c. In such a case, a good choice is to further delay the feedback response setting τ opt = n × 2π/ω. This is exemplified in Fig. 3d, e, where a delay of τ = 4π/ω results in a low final occupancy that is close to for τ < Γ −1 . For larger delay additional non-Lorentzian features in the power spectrum of the position quadrature appear (as described analytically in the next subsection); these features brought in by the feedback delay are not captured by the approximation and regions of insufficient cooling emerge, as can be seen from Fig. 3f. This can be explained by the results in Fig. 3e where we can see that for increasing delay τ the position spectrum deviates strongly from a Lorentzian form by additional superposed oscillations that follow ∼ exp(iΩτ ), where since we are in Fourier space τ determines the period of the oscillations. Thereby, the larger τ is the smaller is the period of the oscillations. If the period is close (see Fig. 3e inset) or coincides with the resonance condition we witness a high occupation even for τ being a multiple of 2π/ω.

B. Fourier analysis of damping rates
Applying a Fourier transform to the coupled set of integro-differential equations allows one to provide an exact analysis of the steady state even in the non-Markovian regime (see Appendix for details). In the frequency domain one can then compute the variances of the position and momentum as Here, we refer to the term S q (Ω) as the position spectrum of the resonator. The integration goes over the whole power spectrum of the noise which can be split into four contributions , stemming from the radiation pressure force S rp (Ω) = G 2 κ/(κ 2 + Ω 2 ), the feedback back-action noise S fb (Ω) = |g (0) (Ω)| 2 /(4κη) as well as its interference S fb,rp (Ω) and most importantly and dominantly from the thermal fluctuations S th (Ω) ≈ γ(2n + 1) (forn 1). The effective susceptibility appearing above describes the modified response of the position quadrature to the external noise and takes a quasi-Lorentzian form The poles are shifted from the original position ±ω to the effective frequency ±ω eff ω 2 eff (Ω) = ω 2 + Ωδω(Ω) cos(Ωτ ) + ΩΓ(Ω) sin(Ωτ ), (36) while the frequency dependent effective damping rate is also modified from γ to γ eff (Ω) = γ +Γ(Ω) cos(Ωτ ) − δω(Ω) sin(Ωτ ).
Notice that in the absence of feedback the susceptibility is simply that of a damped harmonic oscillator with damping rate γ and resonance frequency ω. The feedback for zero delay time adds a damping rateΓ(Ω) and a shifted resonance frequency ω 2 + Ωδω(Ω). Both quantities are then strongly dependent on the time delay as evidenced in Fig. 3. As long as the effective mechanical susceptibility is close to a Lorentzian the integration of the spectra above becomes trivial as the only considerable contribution comes from the spectrally flat thermal power spectrum (see Fig. 3c,e). We will use the fact that the integral 1/(2π) ∞ −∞ dΩ|χ cd eff | 2 ≈ 1/[2γ eff (ω)] for reasonably small delay times where we can setΓ(Ω) ≈Γ(ω) = Γ and δω(Ω) ≈ δω(ω) = δω and we can approximate matching the expression obtained in case of the wFFLC approximation.
In case of the wF-FLC we can witness a regime of instability for {τ |γ + Γ cos(ωτ ) − δω sin(ωτ ) ≤ 0} as is shown in Fig. 3b which will converge to the frequency intervals n × π/2 ≤ τ ≤ n × 3π/2 for n ∈ N and for ω fb , κ ω j and thereby converge to the sFFLC.

VI. SIMULTANEOUS COOLING OF TWO ADJACENT MODES
Let us assume two adjacent resonator modes undergoing simultaneous cold-damping: analytical expressions for the final occupancies can still be obtained from the Lyapunov equation. Under the sFFLC approximation one can express the occupancies as with j = i. The expressions for the off diagonal covariance terms X ij and Z ij as well as the Λ ij terms are more cumbersome and are therefore relegated to Appendix B. From Eq. (39) we see that next to the term expressing the single mode solution, additional terms describing mode to mode coupling emerge, which describe mutual heating effects. For independent modes (such that |ω 1 −ω 2 | → ∞) the mutual heating vanishes. Also notice that as the two modes are subjected to the same damping channel, correlated damping occurs leading to momentum-momentum correlations in the off-diagonal elements. Numerical results involving time domain and Fourier-transform solutions are displayed in Fig. 4. Such effects are particularly evident in the degenerate case with ω 1 = ω 2 . For G 1 = G 2 and g (1) cd = g (2) cd we define the collective bright Q 1 = (q 1 + q 2 )/ √ 2 and dark mode Q 2 = (q 1 − q 2 )/ √ 2 and express them in the Fourier space The susceptibility of the bright mode is modified by the frequency dependent resonance ω 2 B,eff = ω 2 + 2(Ωδω(Ω) cos(Ωτ ) + ΩΓ(Ω) sin(Ωτ )) and damping rate γ B,eff = γ + 2(Γ(Ω) cos(Ωτ ) − δω(Ω) sin(Ωτ )). This results in a damping rate twice as large as compared to the single mode solution found in Eq. (38b) for the bright mode. Instead, the dark mode is fully decoupled from the feedback. We can also analytically list the final phonon occupancies n B,eff (τ ) ≈ 1 2 γ(n + 1/2) γ + 2(Γ cos(ωτ ) − δω sin(ωτ )) × 1 + ω ω + 2(δω cos(ωτ ) + Γ sin(ωτ )) (41a) showing again the decoupling of the dark mode from the feedback loop and that the achievable temperature for the bright mode is half of that of an individual mode.

VII. MULTIMODE COOLING:
The collective basis of one bright and N − 1 dark modes is particularly useful for the case of many resonances. Starting from the equations of motion in Fourier space Eq. (28) we can derive an expression for the bright mode  which solely depends on the noise terms for input. From this result for the bright mode Q 1 the expressions for the N − 1 dark modes Q k can be easily obtained from the relations (details in Appendix C) A simple solution can be found in the case of N identical resonator modes with ω j = ω and γ j = γ. Here, the term for the bright mode given in Eq. (42) becomes where the effective frequency of the bright mode follows ω 2 B,eff (Ω) = ω 2 + j Ωδω j (Ω) cos(Ωτ ) + ΩΓ j (Ω) sin(Ωτ ) with an effective decay rate of γ B,eff (Ω) = γ + jΓ j (Ω) cos(Ωτ ) − δω j (Ω) sin(Ωτ ). This can even further be simplified in the case of identical coupling to the cavity mode G j = G and identical coupling to the feedback force g  dark modes the solutions where ω 2 B,eff (Ω) = ω 2 + N (Ωδω(Ω) cos(Ωτ ) + ΩΓ(Ω) sin(Ωτ )) and γ B,eff (Ω) = γ + N Γ (Ω) cos(Ωτ ) − δω(Ω) sin(Ωτ ) . For zero time delay the bright mode damping rate γ B,eff (Ω) = γ + NΓ(Ω) is N -times larger than in the case of an individual resonator mode, which can be seen from the position spectra in Fig. 5b. The expressions in Eq. (45) show N -uncoupled resonator modes. Here, each mode can be treated independently and following the procedure introduced in section V B and here exemplified for a delay of τ = 0 we obtain for the occupation number of the bright and dark modes showing that only the bright mode experiences cooling, but this mode can reach a much lower occupation number in comparison to the single mode case which is growing with the number N of identical modes that are addressed (see Fig. 5a). In the case that couplings and frequencies vary, the bright mode which is addressed by the feedback mechanism and cooled directly couples with the dark modes and cools them indirectly, resulting in a sympathetic cooling process of the collective modes. More details of the derivation are given in Appendix C.

VIII. ADJUSTABLE TIME DELAY
As in the single mode case, a strategy can be devised to improve the efficiency of cold-damping with time delay by counter-intuitively further delaying the action of the feedback loop. To this end we fix the condition cos ω j τ ≈ 1 by providing an additional time τ add such that τ = τ inh + τ add . The condition springs from the analytical results with Markovian dynamics under the wFFLC or sFFLC approximations (and with a regime of validity implying that τ < (max(Γ j )) −1 is fulfilled for all modes). This demands that ω j τ ≈ 2πn j with n j ∈ N. In the case that all frequencies are multiples of a common frequency ω this can be obtained by bringing the total delay to fulfill τ = 2π/ω. In the general case this is quite a difficult matter and for an increasing number of modes can result in additional delay times τ add exceeding the time regime where the wFFLC and sFFLC approximations are valid, which even can be seen for a single mode presented in Fig. 3f for very large delay times. From color-maps as presented in Fig. 6 showing the total energy of the system it is possible to find regions with efficient cooling and avoid instabilities. From Fig. 6b showing the total energy for two or four modes we also see that that the region where one or many modes are unstable increases for increasing mode number. This results from the fact that many more constraints differing between each mode have to be fulfilled simultaneously. Nevertheless, the condition min τ with ω j τ ≈ 2πn j can be a good guide to find approximate solutions. Additionally one can employ various nonlinear optimization schemes to obtain a minimum with suitable parameters. The results for two-and four modes following a linear dispersion with δω being the frequency difference between two neighboring modes are compared indicating an increase of the area of instability for increasing mode number. We have used the parameters ω1 = 1 with ω k+1 = ω1 + kδω, g cd = 0.6, κ = 4, ω fb = 4.5,n1 = 1 × 10 5 withn k =n1/ω k , where in (a) we chose G = 0.3, 0.2 and γ = (4, 6) × 10 −5 , while in (b) we have used G k = 0.2 + (k − 1) × 0.1 and γ k = (4 + 2(k − 1)) × 10 −5 .

IX. CONCLUSIONS
We have analytically and numerically shown that efficient simultaneous cold-damping of many mechanical resonances is achievable as long as frequency degeneracy is avoided. Furthermore, detrimental effects stemming from the feedback's intrinsically delayed response can be mitigated by introducing an additional variable delay τ add that can be adjusted to optimize the cooling efficiency.
For example, for a sequence of frequencies which are multiples of a common frequency ω efficient cooling is obtained again for a total delay of τ = τ inh + τ add = 2π/ω. Another approach to solve this problem could be the implementation of a machine learning scheme either to find preferable settings for the delay or by using a machine learning procedure to provide for the full feedback mechanism [13]. In the latter case the feedback would adjust the phase to accommodate for the delayed signal and minimize the final temperature.
A main aspect of our treatment is the transformation to a collective basis. In particular for a number of quasi-degenerate modes, a bright/dark mode analysis shows that the damping can be N times faster while the occupancy is N times lower for a suitable identified collective bright mode (as illustrated in Fig. 5a)). This is a remarkable result in itself as it shows that collective optomechanics can be employed to provide more efficient cooling of an engineered collective mode. For applications aiming instead at better sensing capabilities for wider frequency intervals, the alternative, as also indicated in Ref. [14] is to engineer mechanical resonators with a dispersion relation close to linear such that mutual heating is inhibited.
In terms of methods used, the time domain treatment has been very successful in the case of zero delay [14] allowing for fully analytical results for all mode occupancies. The complication of non-zero time delay can be dealt with more efficiently in the Fourier space where final occupancies for an arbitrary number of modes and a large variety of cases can be obtained exactly. Moreover, the effort is computationally less costly for numerical simulations compared to a brute force approach that involves solving the full set of stochastic differential equations.

X. ACKNOWLEDGMENTS
We acknowledge financial support from the Max Planck Society and from the German Federal Ministry of Education and Research, co-funded by the European Commission (project RouTe), project number 13N14839 within the research program "Photonik Forschung Deutschland". We acknowledge initial discussions with David Vitali which have lead to identifying the feedback delay time problem as essential in cold-damping optomechanics. For cold damping with many resonator modes in the quantum mechanical treatment we start with the equations of motion given byq where the effective cavity detuning is kept at zero ∆ = 0. We can eliminate the cavity field quadratures by formally integrating their equations of motion to obtain:

Non-Markovian dynamics
The y est (s) term introduces both terms proportional to the q j as well as noise terms stemming from the cavity input noise y in (s) as well as from the vacuum filled port noise y v . We can first work out the terms coming from y as: cd ω fb y in (s) .
To obtain a dependence with respect to p j we apply integration by parts noticing that the convolution above contains a derivative of the following function: and the relationq j = ω j p j and we obtain We can now write in simplified notation the reduced set of equations of motion for the 2N resonator modes quadratureṡ The three sources of noise are owed to the direct feedback action, to the feedback filtered vacuum action in the loss port and to the intra-cavity radiation pressure effect: with the following definitions φ (τ ) for the convolution kernels.
The fast-feedback-lossy-cavity (FFLC) approximations: Markovian dynamics In the limit of a lossy cavity and relatively fast feedback with κ, ω fb ≥ ω j for all j ∈ {1, . . . , n} where we can approximate p j (s) ≈ A j cos(ω j s + φ j ) and q j (s) ≈ A j sin(ω j s + φ j ) we obtain an estimate to first order given by and we end up with a set of coupled linear differential equations for the mechanical mode quadratures: We refer to this as the weak fast-feedback-lossy-cavity (wFFLC) approximation. Here, the rate terms are defined as and the frequency shift terms are defined as In the case that κ, ω fb ω j the rate terms converge to Γ jk = g (j) cd G k ω k /κ and δω jk = 0 and we obtain the set of coupled differential equations for the mechanical mode quadratures: We refer to this approximation for the dynamics formed by this conditions as the strong fast-feedback-lossy-cavity (sFFLC) approximation. In case we have N -resonator modes with the same frequency ω it is favorable to write the equations of motions in the collective basis where Q 1 = k α 1k q k = l G 2 l −1 k G k q k describes the bright mode while the N − 1 dark modes are given by Q j = k α jk q k with k α ik α jk = δ ij and can be obtained via a Gram-Schmidt procedure. Starting from the wFFLC the equation of motion for the bright mode results iṅ where all noise terms have been gathered into a single term ζ j = ξ j + ξ fb + ξ vac + ξ rp . The N − 1 dark modes are given byQ By solving Eq. (A14) in the case the system approaches steady state we obtain where ω B = ω + k δω kk cos(ωτ ) + Γ kk sin(ωτ ) and Γ B = γ + k Γ kk cos(ωτ ) − δω kk sin(ωτ ). This allows us to reduce the expressions for the dark states toQ where the modified noise terms are given by

Solving the Lyapunov equation
The set of differential equations presented in Eq. (A8a) can be cast into the forṁ with v = (δq 1 , δp 1 , . . . δq N , δp N ) and n in = (0, ζ 1 , . . . , 0, ζ N ). From the general solution we obtain the correlation matrix where we have ignored the transient solution which will decay strongly for large times t. Regarding the noise correlation term n in (s)n in (s ) component wise we obtain n in,i (s)n in,j (s ) = 0 if i and j are both even numbers, which is resulting in for i, j ∈ {1, . . . , N }. Here, the delay τ only appears in the cross terms of the x in , y in correlations. For ω fb , κ ω j , Γ j we can approximate δ(t) ≈ (ω fb /2)e −ω fb |t| as well as δ(t) ≈ (κ/2)e −κ|t| resulting in that shows no dependence on delay. For δ-correlated noise we can simplify the correlation matrix to where D in,2i,2j = (2n i + 1)γ i δ ij + G i G j /κ for even index numbers and is zero otherwise. The Lyapunov equation for the N -oscillator system which determines the steady solution of the correlation matrix is given by Evaluating the individual components we obtain the set of equations with X ij = p i p j + p j p i , Y ij = q i p j + p j q i and Z ij = q i q j + q j q i and additionally we define c i = cos(ω j τ ) and s i = sin(ω j τ ). In the case that γ j (g cd G j ω j )/κ and Γ ij = (g cd )Γ jj , we can simplify the expression in Eq. (26f) and we obtain Here, we define to simplify expressions.
For τ = 0 we can get, with respect to the approximations introduced above, exact solutions for the final energies of each mode as has been reported in [14] and is presented below For τ > 0 it is possible to obtain analytic solutions of Eqs. (A30) for a single and two oscillator modes.
The equations derived from the Lyapunov equation in case of the wFFLC where we need to consider the regime k B T ω j since here the non delta-like noise correlation terms can be ignored are stated below with ∆Γ ij (τ ) = Γ ij c j − δω ij s j and ∆ω ij (τ ) = Γ ij s j + δω ij c j .

Appendix B: Cooling of two adjacent modes
In the sFFLC and under the approximation γ Γ jj carried out in the drift matrix, we can find analytic solutions for two modes. First we express the diagonal elements as This results in for the final occupation of each mode. The off diagonal elements X 12 and Z 12 are given by where K 22 = [(ω 2 + Γ 22 s 2 )ω 2 + ω 2 Γ 11 t 2 c 1 ] and K 11 = [(ω 1 + Γ 11 s 1 )ω 1 + ω 1 Γ 22 t 1 c 2 ] and t j = tan(ω j τ ). i where the driving noise term is given by For high temperatures (k B T ω j ) this can be approximated by ζ j (Ω) ≈ ξ j (Ω).