Optimal control for feedback cooling in cavityless levitated optomechanics

We consider feedback cooling in a cavityless levitated optomechanics setup, and we investigate the possibility to improve the feedback implementation. We apply optimal control theory to derive the optimal feedback signal both for quadratic (parametric) and linear (electric) feedback. We numerically compare optimal feedback against the typical feedback implementation used for experiments. In order to do so, we implement a tracking scheme that takes into account the modulation of the laser intensity. We show that such a tracking implementation allows us to increase the feedback strength, leading to faster cooling rates and lower center-of-mass temperatures.


Introduction
The ability to precisely control and cool the motion of mechanical resonators in order to generate quantum states is of great interest for testing fundamental physics, such as investigating the quantum-to-classical transition [1,2]. A wide variety of resonator systems have shown promise for achieving such goals, including membranes [3,4], micro-and nano-resonators [5,6,7,8] and cantilevers [9,10]. Although ground state cooling as been experimentally realized in optomechanical systems [3,4,8], there is an appetite to reach such states in levitated systems. Levitated nanoparticles are extremely well isolated from their environment, opening up the possibility for very long decoherence times and ground state cooling in room temperature conditions. Indeed, optically levitated silica particles have had their center-of-mass motion cooled to millikelvin [11,12,13,14] and sub-millikelvin [15,16] temperatures, whereas nanodiamonds [17,18] have been used for spin coupling experiments [19,20]. Other levitation mechanisms, such as Paul traps [21], hybrid electro-optical traps [22], and magnetic traps [23,24,25] have also been proposed as candidates for preparing macroscopic quantum states [26,27,28] and testing spontaneous collapse models [29,30]. In order for any of these resonator systems to approach the quantum regime, their motion must first be cooled to close to the ground state, which can be achieved with cryogenically cooling the environment or with active feedback schemes.
In this paper we consider an optically levitated silica nanoparticle, trapped by by the gradient force generated by tightly focusing a 1550nm laser with a high numerical aperture (N.A.) paraboloidal mirror, as shown in figure 1. The optical trap is contained within a vacuum chamber to isolate the particle from its environment as much as possible. Typically, parametric (quadratic) feedback cooling, by modulating the intensity of the trapping laser at twice the particle's oscillation frequency [11,12], is implemented to cool the particle's motion to ∼mK temperatures. Currently, feedback signals are implemented by tracking the phase of the oscillator by locking to the frequency of motion, using either lock-in amplifiers or, more recently, with a Kalman filter [31]. The Kalman filter, a sophisticated filtering technique used in engineering applications [32,33,34,35], can be implemented in real-time to accurately track the particle's motion, before applying the modulating feedback signal [13,36]. Such schemes are very effective for tracking the particle motion for small laser modulation, but above a certain (low) threshold loses track of the particle. This is a limitation as higher modulation results in faster cooling rates and a lower final temperature.
Recently, cooling the motion of charged nanoparticles by applying an electric field which is at the same frequency of the particle's motion has been demonstrated [16,37] and implemented with optimal control protocols [38] for optical traps, as well as proposed for electrical traps [26]. A charged needle, placed in the vacuum chamber close to the laser focus, has been used for force sensing applications [39] and investigations of Fano resonances [40] in levitated optomechanics. To first approximation, the electric field generated by the needle couples linearly to the particle position, making it suitable to  Figure 1.
The experimental setup that we are simulating. The position of the particle is detected by interference between the scattered and divergent light at the photodetector. This detected signal is then passed into an oscilloscope for recording and a field programmable gate array (FPGA) to perform the tracking and generate the feedback signals which can then be sent to the acousto-optic modulator (AOM) to modulate the light and perform quadratic feedback or to the needle to perform linear feedback.
implement linear feedback cooling. By applying a force to oppose the particle motion, the amplitude of motion can be reduced. It is worth noting that for this cooling technique the coupling strength cannot be indefinitely high, as too strong an applied force would drive the particle to hotter temperatures, and could even result in the particle being ejected from the trap.
In this article we consider whether it is possible to implement a feedback protocol which takes into account all the contributions to the particle dynamics, including decoherence and photon recoil, and compare to current feedback schemes discussed previously. We utilize optimal control theory to investigate both quadratic (parametric) and linear (electric) feedback (section 3). Optimal control theory has been applied to other experimental systems [41], including for manipulation of Bose-Eintsein condensates to prepare complex quantum states [42], designing excitation pulses in NMR [43] and tailoring robustness in solid-state spin magnetometry [44]. Additionally, it has been proposed for mixed state squeezing in cavity optomechanics [45], feedback cooling and squeezing of levitated nanopshperes in cavities [46] and recently for feedback cooling in low frequency magnetic traps [27]. To compare optimal cooling with typical feedback cooling, we numerically emulate the system, by solving its equations of motion, and we track it by numerically solving a second set of equations (section 4). This tracking technique is more rubust, for our applications, than the one performed by a basic timedependent Kalman filter, and allows us to increase the laser modulation to achieve better cooling of the trapped nanoparticle. Sections 4.1 and 4.2 are respectively dedicated to quadratic and linear feedback results, while section 4.3 concerns the common features of the two cooling schemes. In the next section we start by introducing the theoretical framework of the setup considered, and of its numerical simulation.
We remark that, although the analysis presented concerns a gradient force optical trap, the improved tracking developed is flexible, and could be also implemented in the other optomechanical setups previously described.

Dynamical models
The optically levitated nanoparticle undergoes continuous monitoring of its motion by the trapping laser [39,47]. Specifically, we consider the experimental situation when the translational and rotational degrees of freedom are decoupled, and a single translation degree of freedom can be identified in the detected signal [48]. We will label the position of this one-dimensional motion with z. We call J the homodyne current that is physically accessible with the experimental setup, i.e. the quantity recorded by the measurement apparatus.
We consider two types of dynamical modelling of the experiment: (i) a tracking model and (ii) an emulation model [49]. In a nutshell, the models of type (i) are used to track the state of the system, and to provide its best estimate, either in real-time or in post-selection, given the input homodyne current J measured by the experimental setup [50,51]. The only role of the models of type (ii) is to generate a trajectory of the system and to output the homodyne current J, i.e. to emulate the system when the experimentally measured J is not available. This proves useful if one wants to investigate numerically the efficacy of a technique before implementing it in the experimental setup.
We now discuss in detail the tracking (section 2.1) and emulation (section 2.2) model, while the feedback mechanism will be discussed in section 3.

System tracking
The dynamics of the continuously monitored trapped particle is described by the following master equation [52]: the first being the harmonic trap provided by the laser, and the second being the feedback Hamiltonian, where u(t) and v(t) are feedback signals that depend on the particle measured position and momentum. The second and third terms of equation (1) describe respectively decoherence and monitoring provided by laser photons, with detection efficiency η and monitoring strength k. The fourth term accounts for decoherence due to residual gas particles, with Γ = 4mk b T γc 2 andL =ẑ +i γc Γp [53,54,55]. We will refer to this model as the tracking model and to the corresponding state by ρ T , i.e. the tracked state. We assume that the particle is described by an initial Gaussian state. Since the dynamics of equation (1) is quadratic, the system state remains Gaussian during its evolution, i.e. it is fully described by mean values of position and momentum operators ( ẑ T , p T ) and their variances ( . By exploiting equation (1) one finds that these evolve according to the following equations: These equations account for the fact that the monitoring strength k is proportional to the laser power, that is modulated by the feedback signal u(t): 0 ω L is the coupling strength that depends on the laser power P [13].

System emulation
In the previous section we have discussed how to track the motion of an optically levitated system. Specifically, we have assumed that we are given an experimentally measured homodyne current J which is then used to continuously update our knowledge of the system. However, to emulate the system and to generate an output homodyne current J we have to use a modified dynamical model. In particular, we will rewrite the update term due to the detected photons as a stochastic back-action term and we will include additional stochastic terms to account for the undetected photons as well as gas collisions [47,49,56]. This is fully analogous to a classical emulation of the system: loosely speaking, each scattering event due to photons (even if undetected) or to gas particles makes the particle recoil, and this is modelled by noise terms. We will refer to such a model as the emulation model and denote the corresponding state of the system by ρ E , i.e. the emulated state.
To emulate the system we consider the following dynamical equation: where the first line (from left to right) includes the unitary evolution, the diffusion terms due to gas scattering, the diffusion term due to photon scattering, and the stochastic back-action term due to the detected photons, where dW is a Wiener process with zero mean and correlation E[dW dW ] = dt. The third line accounts for the nanoparticle recoil due to undetected photons (∝ dV ) and gas particles (∝ dZ), where dV and dZ are additional independent Wiener processes with zero mean and correlations set to This latter term ∝ dZ is responsible for thermalizing the motion of the nanoparticle with the gas particles. The associated homodyne current is given by: where We now, similarly as for the tracking model, limit the discussion to Gaussian states and introduce the vector y = (y 1 , y 2 , y 3 , y 4 , In this case equation (8) can be reduced to the following coupled set of stochastic differential equations:ẏ

Optimal feedback
Our aim is to determine the optimal controls u * (t) and v * (t) (here and in the following the asterisk denotes the optimal realization of a function) that provide the best cooling of the trapped particle, i.e. that minimize mean energy Note that, although Ĥ 0 does not depend explicitly on the control functions u(t) and v(t), y does. Since such a dependence is linear, LQG optimization cannot be applied [57,58], and one needs to tackle the problem differently. We exploit Pontryagin's Minimum Principle (PMP), an important tool of optimal control theory, that allows to find the optimal control that minimizes a given cost function [59]. The one solved by the PMP is a minimization problem with constraint (given by the equations of motion (3)- (7)). It is convenient to introduce the 'co-states vector' λ ≡ (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ) and to define a 'co-state Hamiltonian' as follows: One can check that the evolution equation for the co-states iṡ while equations (3)- (7) can be conveniently rewritten as follows: Pontryagin's principle precisely states that the optimal control u * , v * are those such that Since the equations of motion for the components of y are linear both in u and v, one can check the optimal signals satisfying the condition (19) are where sgn is the sign function that is 1 (−1) when its argument in positive (negative). We remark that the sign function form of the control has a simple intuitive explanation: in the case of quadratic feedback one wants to stiffen (weaken) the trap maximally when the particle moves away (towards) the trap center, and, similarly in the case of linear (electric) feedback one would like to stop the particle in its motion by applying the maximum "breaking" force. The only restriction is thus on the the feedback strength, i.e. on the modulation depth: if the modulation depth is too strong one risks loosing track of the particle,or worse, loosing the particle from the trap. In order to obtain the explicit expressions for the two control functions, one needs to solve the two coupled sets of equations (17)- (18). This is in general a hard task because, while the state equations (18) have initial boundary conditions and propagate forward in time, the costate equations (17) have final boundary conditions (at the measurement time ∆t) and propagate backward in time [59]. As we will discuss in the next section, under certain conditions it is numerically convenient to adopt a different strategy instead of solving the co-state equations.

Numerical analysis of feedback schemes
To emulate the system we first discretize equations (9)- (14), i.e. we consider a time step ∆t E , and the Wiener increments ∆W , ∆V , ∆Z. We set the initial state to be a Gaussian thermal state, i.e.
and y 1 = y 2 = y 5 = 0, where T is the temperature of the gas particles and k B is Boltzmann's constant. Specifically, we set T = 300K, i.e. we assume that the gas of Schematic diagram of the numerical simulation. The emulation part (purple) consists of the system and detector: the system is a Gaussian state described by the vector y and evolves according to equations (10)- (14) with three inputs, i.e. the control (u and v), the environmental noise (dV and dZ), and the back-action noise (dW ), while the detector has two inputs, i.e. the state y and the imprecision noise (dW ), and produces the output current J given by equation (9). The estimation and control part (blue) consists of the tracking and actuator: the estimator of the state consists of the vector x and evolves according to equations (3)-(7) with two inputs, i.e. the control (u and v) and the current J, while the actuator controls the functions u and v in response to the best estimate of the system given by x. In an experimental realization the emulation part and simulated current J are replaced by the experiment and the experimental current, respectively; the estimation and control part remains unchanged.
particles is at room temperature. We set ω = 2π × 70kHz, and m = 9.42 × 10 −19 kg, which are typical values of trapped dielectric silica particles of radius ∼ 50nm. We also set η = 0.003 and α = 4.04 × 10 25 which are typical values for our optical trap [11]. We then propagate the initial state for a time t = t prep ∼ 5ms, with the control functions set to u = v = 0: this preparation procedure ensures that the state of the system at time t = t prep is more realistic, i.e. the noise and dynamics of the emulation model will drive the system to a new state y(t prep ). It also allows the tracking to converge and begin tracking the system well. Specifically, we solve the stochastic differential equations in equations (9)-(14) using the fourth-order stochastic Runge-Kutta method.
To track and control the system we need to solve in parallel also the equations (3)-(7), as well as choose the control functions u and v. However the experimentally available time-step for the tracking and control is limited by the apparatus, e.g. sampling rates, reaction times and time lags. It is thus reasonable to consider larger time-steps, ∆t = M ∆t E and ∆t C = ∆t N , for the measurement and tracking/control respectively, with N, M ∈ N . To simulate the current experimental capabilities presented in [13] we set N = 5, M = 2000 and ∆t E = 0.5ns. We have verified numerically that such a value of ∆t E provides with enough temporal resolution to simulate sufficiently well the evolution of the system.
We set the initial state of the tracking at time t = 0 to be a Gaussian thermal state, i.e. x(0) = y(0), where the non-zero values of y(0) are given in equations (21). We switch on the feedback control at time t = t prep . In order to avoid the difficulties of propagating the co-states equations backward in time, we adopt the following strategy. At each step ∆t we select the optimal control by selecting the optimal trajectory: we propagate the estimated state x forward in time for ∆t using the time-step ∆t C = ∆t N for each possible trajectory of the controls u or v. Since according to equation (20) the value u can have only values ±1 this amounts to 2 N trajectories; the same applies also for v. If both u and v would be controlled simultaneously in such a way we would have a total of 4 N trajectories. We select the trajectory that minimizes the cost function given in equation (15), i.e. the one that minimizes the estimated energy. It turns out that, at least for low values N , the parallelization of the optimal control problem is computationally feasible. In particular, the scenario investigated here is particularly relevant for experiments involving FPGAs; for example, setting N = 5 gives a total of 2 5 = 32 trajectories for a single control function, which is readily solved in parallel using even moderately priced FPGAs.
The schematic diagram in figure 2 gives an overview of the emulation-tracking implementation; for a more detailed introduction see e.g. [58,49]. The feedback details for quadratic and linear cases will be respectively discussed in the following subsections 4.1 and 4.2. Section 4.3 is devoted to the discussion of common features of the two feedback schemes.

Quadratic feedback
Parametric (i.e. quadratic) feedback is widely used in levitated optomechanics. The relevant equations describing this type of feedback can be obtained simply setting δ = 0 in section 2. This type of feedback is typically performed by modulating the laser at twice the phase of the particle, setting u = ω E x 1 x 2 [13], where E = x 2 2 2m + mω 2 2 x 2 1 . For a fixed value of β we can then directly compare the optimal control u * with the simple double phase u. Specifically, figures 3 and 4 show that the cooling obtained with the double phase modulation is effectively equivalent the optimal feedback; this is explained by the fact that the feedback time trace for the two cooling approaches is almost the same (see inset in figure 3). The difference between the sine profile and the square-wave function does not substantially affect the magnitude of the cooling.
However, there is one important difference between the basic tracking technique (exploited in double phase cooling) and the new tracking (exploited in optimal cooling). The first is performed via a Kalman filter that simulates a phase-locked loop (PLL) by exploiting equations (3)-(7) with β = 0 (unmodulated tracking). The latter instead makes use of equations (3)-(7) with the same β as in equations (10)-(14) (modulated tracking), and requires a fully FPGA-based implementation. One of the limitations of the typical (unmodulated) implementation of parametric feedback is that one can reach only strength of about β = 0.01. This is due to the fact that for higher values of β the tracking loses track of the particle because of the larger frequency variation of the system and therefore cooling is not as effective, and for higher modulation depths even heats the system instead of cooling, see Figure 5. One of the merits of our improved   tracking scheme is that it allows us to track the system trajectories for higher modulation strength. Figures 4 and 5 clearly show that β can be increased, allowing faster cooling rates and lower particle temperatures to be obtained.

Linear feedback
Linear feedback can be implemented in the experimental setup by inserting into the vacuum chamber a needle to which an electric voltage is applied [40]. The electric force generated by the needle affects the particle motion as described in section 2, i.e. we set β = 0 and modulate the control function v. Figure 6 compares the optimal linear feedback with cold damping, i.e. a force proportional to the particle velocity, v = y 2 m showing that the latter is comparable to the performance of the optimal feedback. The explanation is the same as the quadratic case: the difference between the square optimal signal and the sinusoidal velocity does not significantly affect the cooling efficiency. An important issue one needs to account for while using linear feedback is that the force kicking the particle might lead to reheating and particle loss if they are too strong. For this reason the simulation of figure 6 makes use of a δ in the "optimal range" identified in figure 7.

Discussion
From our simulation and analysis, it is found that low phonon number states can in principle be reached with both quadratic (figure 3) and linear (figure 6) feedback protocols. We find that the cooling signal obtained via optimal control theory does not outperform typical feedback cooling. This can be explained by the fact that our knowledge of the system is given only by the (measured) position of the particle, and this contains all the information about its dynamics (including decoherence and recoil effects). Since typical feedback controls are based on the measured position and velocity of the particle (i.e. the difference of two subsequent positions), they already contain our best knowledge on the system, and the control shape does not play a crucial role. We remark here that the equations for the mean values (equations (3)-(4) and (10)- (11)) are essentially decoupled from the equations for the variances (equations (5)- (7) and (12)-(14)), the only connection being the arguments of the control functions (see equation (20)). Accordingly, for all practical purposes it is enough to use only equations (3)-(4) and (10)- (11).
We find that for quadratic feedback, increasing the modulation depths β decreases the minimum phonon number (or temperature) achievable, and for linear feedback there is an 'optimal range' of coupling strength δ for when cooling is most effective. Above this optimal region the particle will be cooled less effectively, and eventually heated, as δ increases. The achievable phonon number for increased modulation depth (coupling strength) can be seen in figure 7.
Increasing the cooling strengths β and δ also increases the initial rate that the nanoparticle is cooled, as expected. The initial cooling rate as a function of cooling strength for both cases can be seen in figure 8. As the quadratic modulation depth β increases the initial cooling rate increases linearly, whereas increasing the linear cooling strength δ results in an non-linear increase in the initial cooling rate that approaches an asymptotic value. It was found numerically that this is because the equations (10)- (14) with β = 0, δ = 0 (quadratic feedback) require a smaller sampling rate than the equations (10)- (14) with β = 0, δ = 0 (linear feedback) in order to track the system well. By increasing the sampling rate sufficiently, it was found that the cooling rate also increases linearly for linear feedback, as is shown in the inset of figure 8. Experimentally, where sampling rates cannot be set arbitrarily high, this has practical implications in cooling rates that can be achieved via linear feedback.
In previous experimental works, where parametric feedback has been utilized, the tracking of the particle's motion does not take into account the laser intensity modulation due to the feedback, which results in a maximum modulation depth of ∼ 1.5% [11], after which the effectiveness of cooling decreases, eventually causing heating. This is due to the fact that this tracking does not take into account the varying laser intensity, which effectively changes the oscillation frequency of the particle for a fraction of an oscillation period. For small modulation depths this isn't an issue as the effective frequency is still within the tracking bandwidth of the tracking mechanism, but for larger modulations results in the feedback being out of phase. The tracking scheme presented here overcomes this limitation by factoring in the laser modulation in the tracking equations, allowing access to extremely high modulation depths and cooling rates with quadratic feedback.
After cooling, it was found that the contribution to the mean energy of the particle is dominated by the expectation values of the position and velocity, (x 1 , x 2 in equation (15)) whereas the variances' (x 3 , x 4 ) contribution is found to be negligible. The energy contained in the variances quite rapidly achieve the Heisenberg limit x 3 x 4 = /2, and the energy fluctuations are due to random collisions with gas particles and photons. It is interesting to investigate which of these has more of an effect of the particle dynamics. Note that since the variances are constant , equations (5)-(7) are essentially a "time-dependent Kalman filter" [60,61]. Figure 9 shows that for high pressures the phonon number is mostly affected by the gas particles collisions. This effect can be made negligible by reducing the gas pressure in the vacuum chamber. However, when the pressure is sufficiently low one reaches the photon recoil limit, the regime where fluctuations are mostly given by the photon scattering. This kind of effect is always present in the experiment, and it ultimately limits the achievable phonon value. One might try to decrease photon scattering by lowering the laser power, but this leads to a less stable trapping and to a weaker detected signal.
We eventually remark that linear and quadratic feedback can be combined and used at the same time although in our investigations with modulation depth around 1.5% this has not significantly altered the final temperature.

Conclusions
We have considered a cavityless levitation experimental setup, and we numerically investigated both parametric (quadratic) and electric (linear) optimal feedback cooling. The comparison of optimal feedbacks against the typical implementations (respectively double phase and cold damping) show that, although the feedback profiles are .
Shows the dependency of the average phonon number reached once the energy has converged with the pressure at which the simulation is performed for optimal and double phase quadratic feedback at three different feedback modulation depths. At about ∼ 10 −8 mbar we reach the photon-recoil regime [15]. different, this does not substantially affect the magnitude of the cooling. However, the implementation of optimal feedback forced us to develop a more sophisticated tracking scheme. This allowed us to go beyond one of the limitations of the typical implementations of parametric feedback, namely the low modulation strength limit. One of the merits of the more sophisticated tracking scheme is that it allows us to increase the modulation strength, obtaining a faster cooling rate and reaching lower temperatures. Furthermore, it was found that for linear feedback there exists an 'optimal range' for the coupling strength, that provides most effective cooling. Moreover, although combined (quadratic+linear) feedback does not seem to significantly alter the achieved temperature, it might still be experimentally helpful to make the cooling more stable. Further improvement might be obtained by applying functional non-Markovian techniques [62] to optimal control theory, in order to account for experimental time lags in the derivation of the optimal feedback.