Controlling spontaneous-emission noise in measurement-based feedback cooling of a Bose–Einstein condensate

We consider the simulation and analysis of a BEC undergoing off-resonant imaging and active feedback control. Continuously monitoring a BEC provides a signal that contains information regarding the density of the BEC. The measurement can be used not only to estimate the observable being measured, but also to give a best estimate (in the least-squares sense) for the full quantum state of the condensate. This can be achieved using a conditional master equation [49, 56, 57], which is better known as a quantum filtering equation in the engineering community [53, 58, 59]. The conditional master equation of a BEC undergoing continuous off-resonant imaging is governed by the Stratonovich equation [13–16, 32, 60]

$$\begin{equation} \partial_t \skew3\hat{\rho} = -\frac{\mathrm{i}}{\hbar}[\skew3\hat{H}(\mathbf{ u}),\skew3\hat{\rho}] + \sum_i \mathcal{D}[\skew3\hat{L}_i] \skew3\hat{\rho} + \sum_i \mathcal{C}[\skew3\hat{L}_i] \skew3\hat{\rho} + \sum_i \mathcal{H}[\skew3\hat{L}_i] \skew3\hat{\rho} \; \eta_i(t), \end{equation} \tag{ 1 }$$

where $\hat {\rho }$ is a density matrix encoding the full quantum field of the condensate conditioned on the measurement result. The first term describes the unitary dynamics of the system, which is determined by the Hamiltonian

$$\begin{equation} \skew3\hat{H}({\bf u}) = \int {\rm d}{{\bf x}} \; \hat {\psi }^\dagger({{\bf x}}) h({{\bf x}},{\bf u}) \hat {\psi }({{\bf x}}) + \frac{U}{2} \int {\rm d}{\bf x} \, \hat {\psi }^\dagger({{\bf x}})\hat {\psi }^\dagger({{\bf x}}) \hat {\psi }({{\bf x}})\hat {\psi }({{\bf x}}), \end{equation} \tag{ 2 }$$

where $\hat {\psi }(\mathbf { x})$ is the field operator that annihilates an atom in the ground state at position x (we assume that the imaging light has been sufficiently detuned from the atomic resonance such that the excited state can be adiabatically eliminated). The field operators obey the commutation relations $[ \hat {\psi }(\mathbf { x}), \hat {\psi }^\dagger (\mathbf { x}')] = \delta (\mathbf { x} - \mathbf { x}')$, where δ(x − x') is a Dirac delta function. The first term in equation (2) is the single-atom energy, governed by the single-particle Hamiltonian h(x,u). Feedback control is included in this Hamiltonian via the vector of control signals, u(t), which in general depends upon the system state $\hat {\rho }$ at time t. The second term in equation (2) is the energy due to the two-body collisions between atoms, assuming a hard-sphere contact scattering potential of strength U. This is an excellent approximation for gases of ultracold alkali atoms [61].

The effect of the measurement is described by the remaining terms in the conditional master equation (1). The superoperators

$$\begin{equation} \mathcal{D}[\skew3\hat{c}]\skew3\hat{\rho} \equiv \skew3\hat{c}\skew3\hat{\rho} \skew3\hat{c}^\dagger -\textstyle\frac{1}{2}\left(\skew3\hat{c}^\dagger \skew3\hat{c}\skew3\hat{\rho}+\skew3\hat{\rho} \skew3\hat{c}^\dagger \skew3\hat{c}\right), \end{equation} \tag{ 3a }$$

$$\begin{equation} \mathcal{H}[\skew3\hat{c}] \skew3\hat{\rho} \equiv \skew3\hat{c} \skew3\hat{\rho} + \skew3\hat{\rho} \skew3\hat{c}^\dagger - \langle \skew3\hat{c}+\skew3\hat{c}^\dagger \rangle \skew3\hat{\rho}, \end{equation} \tag{ 3b }$$

$$\begin{equation} \mathcal{C}[\skew3\hat{c}] \skew3\hat{\rho} \equiv \langle \skew3\hat{c} + \skew3\hat{c}^\dagger \rangle \mathcal{H}[\skew3\hat{c}]\skew3\hat{\rho} -\textstyle\frac{1}{2} \mathcal{H}[\skew3\hat{c}^2]\skew3\hat{\rho}+\langle \skew3\hat{c}^\dagger \skew3\hat{c} \rangle \skew3\hat{\rho}-\skew3\hat{c} \skew3\hat{\rho} \skew3\hat{c}^\dagger \end{equation} \tag{ 3c }$$

$\hat {c}$

$\langle \cdot \rangle = {\textrm {Tr}}[\cdot \hat {\rho }]$

$$\begin{equation} \skew3\hat{L}_i = \int {\rm d}{{\bf x}} \; \hat {\psi }^\dagger({{\bf x}}) l_i({{\bf x}}) \hat {\psi }({{\bf x}}), \end{equation} \tag{ 4 }$$

where l_i(x) are the density-moment generators, which depend upon the precise setup of the imaging optics. Note that since off-resonant imaging is a scattering process, the measurement operators are restricted to the form of equation (4) to ensure that the particle number of the quantum gas is conserved. The stochastic nature of the random wavefunction collapse due to the set of measurement operators is included via the noises η_i(t) (more on these shortly). The validity of the conditional master equation (1) is corroborated by previous theoretical examinations into the off-resonant imaging of BEC by Dalvit et al [60], Szigeti et al [14, 15] and others [12, 62–65].

Physically, the conditional master equation (1) gives an estimate of the condensate as a function of the measurement signal(s). In principle, this conditional master equation could be integrated during an experiment to give a real-time estimate of the current quantum state of the BEC. Furthermore, this information could then be used to apply the required controls to the condensate, which would complete an active feedback loop. In this case, η_i(t) corresponds to the difference between each measurement signal Y_i(t) and the current estimate of the observable being measured [52]. Explicitly, if the measurement signal changes by ΔY_i within a time interval Δt, then

$$\begin{equation} \Delta \eta_i = \Delta Y_i - {\textrm{Tr}}[ (\skew3\hat{L}_i + \skew3\hat{L}_i^\dagger) \skew3\hat{\rho}] \Delta t. \end{equation} \tag{ 5 }$$

However, in this paper we only simulate equation (1). Fortuitously, this can be achieved by using equation (1) and simply treating η_i(t) as a Stratonovich stochastic integral [47, 49, 53, 66]. More specifically, we sample Δη_i(t) from a Gaussian distribution with variance 1/Δt, and then numerically integrate equation (1) using an algorithm appropriate for Stratonovich SDEs (for some examples, see [67]). Each stochastic path of η_i(t) now corresponds to a virtual measurement record corresponding to a single run of an 'experiment'.

In practice, we simulate equation (1) many times (each with a different realization of the stochastic noises, η^(p)_i(t), resulting in different density matrices, $\hat {\rho }^{(p)}$, here indexed by p) and then average these paths (or trajectories) to get a result. We use the notation $\mathbb {A}[ \cdot ]$ to mean this procedure. For example, in order to find the average density over time, we simulate equation (1) P times and then compute

$$\begin{equation} \mathbb{A}\left[\langle \hat {\psi }^\dagger({{\bf x}}) \hat {\psi }({{\bf x}}) \rangle \right] = \frac{1}{P}\sum_{p=1}^P {\textrm{Tr}}\left[ \hat {\psi }^\dagger({{\bf x}}) \hat {\psi }({{\bf x}}) \skew3\hat{\rho}^{(p)}\right]. \end{equation} \tag{ 6 }$$

These ensemble averages describe the typical behaviour of the system. They are generally more useful to examine than individual stochastic paths, as conclusions based on individual paths may be spurious due to the stochastic nature of the paths. However, we emphasize that every integration of equation (1) is physically meaningful, and thus each $\hat {\rho }^{(p)}$ can be thought of as a new 'experiment'. We note that if there was no feedback then $\mathbb {A}[\hat {\rho }]=\hat {\varrho }$ in the limit as P → ∞, where $\hat {\varrho }$ is the density matrix governed by the unconditional version of equation (1). This is not the case when feedback is active, since feedback can produce non-Markovian behaviour that cannot be modelled by an unconditional (or traditional) master equation.

Unfortunately, simulating equation (1) directly and exactly is impractical for even a modest number of particles. Suppose we limit the number of modes we simulate to M, and then truncate the number of excitations in each mode to J. Then the field operator can be written as $\hat {\psi }(\mathbf { x}) = \sum _{m=1}^M u_m(\mathbf { x}) \hat {a}_m$, where $\left \{u_m(\mathbf { x})\right \}_m$ are a set of orthonormal functions, and each mode $\hat {a}_m$ has a set of J energy eigenstates labelled |j_m〉. The density matrix is then

$$\begin{equation} \skew3\hat{\rho} = \sum^{J}_{ {j_1,\ldots, j_M = \, 0} \atop {j_1',\ldots, j_M' =\, 0}} c_{j_1, \ldots, j_M,j_1', \ldots, j_M'} |j_1, \ldots ,j_M \rangle \langle j_1', \ldots ,j_M'|. \end{equation} \tag{ 7 }$$

The memory required to store this density matrix is J^2M times the memory required to store each complex number. This exponential scaling of the size of the quantum field with the number of modes is an explicit example of the curse of dimensionality that generally restricts the direct simulation of large quantum systems. Both the Hartree–Fock method and NPW particle filter are approximate techniques that circumvent this problem. However, as shown shortly, simulations of equation (1) with the NPW particle filter retain more of the important quantum statistics than simulations using the semiclassical Hartree–Fock approximation.

In many BEC experiments, the dynamics of the condensate depend minimally on correlations in the atomic field, and can thus be well-described by a mean-field theory. This is analogous to the case in classical optics, where the quantum fluctuations in the electric field can be ignored. There are a number of different approaches to constructing a semiclassical theory of BEC dynamics, which in most situations give identical equations of motion. Common to all approaches is the assumption that all the atoms of the condensate occupy the same single-particle quantum state. This allows the system to be modelled as a single macroscopic wavefunction (or order parameter) with a size independent of the total number of particles, thereby circumventing the curse of dimensionality associated with the full quantum field.

The semiclassical theory used in this paper is based on the Hartree–Fock approximation. This assumes that (a) the BEC is always in a state of fixed total number N, and (b) that all N atoms occupy the same mode. This restricts the conditional density matrix in equation (1) to the form

$$\begin{equation} \skew3\hat{\rho}_{\rm HF}(t) =\bigotimes_{i=1}^N \left(\int {\rm d}{{\bf x}} \, {\rm d}{\bf y} \, \phi^*({{\bf x}},t)\phi({\bf y},t) |{{\bf x}} {\rangle}_{ii}{\langle}{\bf y} |\right), \end{equation} \tag{ 8 }$$

where $\hat {\rho }_{\mathrm { HF}}$ is written in first quantized form with |x〉_i the position eigenstate of the ith atom and ϕ(x,t) the position-basis wavefunction. This approximation fixes all the quantum statistics of the condensate; for example

$$\begin{equation} \langle \hat {\psi }^\dagger({{\bf x}}) \hat {\psi }({\bf y}) \rangle = N \phi^*({{\bf x}},t)\phi({\bf y},t) \end{equation} \tag{ 9 }$$

and

$$\begin{equation} \langle (\hat {\psi }^\dagger({{\bf x}}))^2 (\hat {\psi }({{\bf x}}))^2 \rangle = N(N-1) |\phi({{\bf x}},t)|^4. \end{equation} \tag{ 10 }$$

Thus, we cannot describe squeezed states, thermal states or any other exotic quantum states with the Hartree–Fock approximation. As described in Szigeti et al [15], equations (8) and (9) can be used to find the equation of motion for the one-body correlation matrix (9), which then gives the following equation of motion for the unnormalized wavefunction:

$$\begin{eqnarray}&&\fl \partial_t \tilde\phi({{\bf x}},t) = \left\{ -\frac{\mathrm{i}}{\hbar}\left(h({{\bf x}},{\bf u}) + (N-1)\frac{U}{n_\phi(t)} |\tilde{\phi}({{\bf x}},t)|^2\right)\right. \nonumber \\ &&\left.+ \sum_i \left(2 l_i({{\bf x}}) L_i^\phi(t)-l_i({{\bf x}})^2 + l_i({{\bf x}}) \; \eta_i(t)\right) \right\} \tilde{\phi}({{\bf x}},t), \end{eqnarray} \tag{ 11 }$$

where $\tilde {\phi }(\mathbf { x},t)$ is related to the normalized wavefunction by

$$\begin{equation} \phi({{\bf x}},t) = \frac{\tilde{\phi}({{\bf x}},t)}{\sqrt{n_\phi(t)}} \end{equation} \tag{ 12 }$$

and we have defined

$$\begin{equation} n_\phi(t) \equiv\int {\rm d}{{\bf x}}\, |\tilde{\phi}({{\bf x}},t)|^2, \end{equation} \tag{ 13 }$$

$$\begin{equation} L_i^\phi(t) \equiv \int {\rm d}{{\bf x}}\, l_i({{\bf x}}) \frac{|\tilde{\phi}({{\bf x}},t)|^2}{n_\phi(t)}. \end{equation} \tag{ 14 }$$

Numerical integration of the unnormalized evolution equation is preferred as it is much faster and possesses increased numerical stability over the normalized equation of motion. Equation (11) must be integrated for each virtual measurement record η_i(t), although this does not need to be done simultaneously. Each integration corresponds to a path or an 'experiment'. We note that this equation reduces exactly to the equation of motion for a single particle when N = 1.

The Hartree–Fock approximation is only valid for condensates well below the phase transition temperature. Even if this condition is satisfied for the initial state, it may cease to be true due to the heating caused by the off-resonant imaging. However, we will see that for measurements that probe the condensate weakly, the Hartree–Fock approximation remains valid.

It might be wondered what advantages the Hartree–Fock method has over the coherent-state-based mean-field approximation that is traditionally used to derive the Gross–Pitaevskii equation (GPE), given that both approximations result in simulations of identical dimensionality. The answer is that applying a coherent-state-based mean-field approximation to equation (1) produces a GPE-like equation that is numerically unstable [15]. This is thought to be due to the effect of the measurement. Off-resonant imaging rapidly gathers information about the total number of atoms in the condensate, projecting the system into a subspace with a well-known total atom number. The Hartree–Fock approximation has a fixed total number of atoms and continues to be valid during this measurement process. In contrast, the mean-field approximation does not have a fixed total number of atoms, and thus gives an inaccurate description of the system. Note that in the case when a BEC is not being monitored, the Hartree–Fock and mean-field approximations are virtually identical. For if we set l_i = 0 and assume that N ≫ 1, then equation (11) reduces to the GPE. Despite these important differences, both the Hartree–Fock approximation and coherent-state-based mean-field approximation are in the same class of semiclassical approximations where the quantum statistics are fixed. We now examine a simulation method that does not fix the quantum statistics, and is therefore able to simulate a much larger class of physical phenomena.

The NPW particle filter is currently the most precise method for simulating the long timescale dynamics of a continuously monitored BEC without fixing the quantum statistics [30]. Instead of integrating a single complex function, as required by the Hartree–Fock approximation, the NPW particle filter requires the simultaneous integration of multiple complex functions, each with its own real-valued weight and additional stochastic noise. The weights determine which one of the complex functions best match the current measurement record, while the noise acts as a correction for the backaction due to the measurement. Observables are then calculated by using a weighted average over these complex functions.

The approximations used in the derivation of the NPW particle filter are detailed in appendix A. Like the TW phase-space method [33–37], it is valid when the number of atoms N is much larger than the number of modes M being simulated. For coherent evolution, the NPW particle filter and TW phase-space method produce essentially identical simulations [31]. TW can also be used to simulate open quantum systems, an early example being [68]. However, the NPW particle filter is more numerically stable and accurate than TW when a system is continuously monitored. This is because the TW representation generates non-semi-positive definite diffusion when applied to equation (1), which makes it numerically unstable [31, 32]. This might come as a surprise, since TW is guaranteed to produce semi-positive definite diffusion when applied to deterministic master equations. However, the assumptions required for this proof fail to hold in the case of stochastic master equations. In contrast, the NPW particle filter has no such issue. We can think of the NPW particle filter as being an additional member of the so-called c-field techniques for simulating the non-equilibrium dynamics of BECs [34, 69], and one which is of most use when the condensate is being monitored.

Phase-space methods typically require multiple integrations of stochastic equations to produce quantum averages. However, the addition of monitoring to the problem requires a set of weighted stochastic differential equations (WSDE) that are simulated in parallel. The concept of a weighting and weighted averages is common in statistics and meta-studies [70–72]. Concepts similar to weights have also been used in phase-space methods [73, 74], quantum measurement theory [75] and quantum control [76]. Of the many applications, the techniques that are closest in spirit to our approach are so-called particle filters [43–46], which have become a leading technology in the field of engineering for classical filtering problems where Gaussian approximations do not apply. Examples include image tracking [77–79], GPS navigation [80] and financial modelling [81, 82]. Particle-filter-like simulations appear to have first been applied to quantum systems in the context of simulating monitored BECs [83]. However, very soon after this publication they were also applied by Jacobs [84] in the simulation of conditional master equations with imperfect detection efficiency. More generally, the emergence of quantum control as an important topic in both theory and experiment is seeing a growth in the application of advanced techniques from engineering, such as particle filters, to quantum mechanical problems.

The NPW particle filter corresponding to equation (1) is derived by mapping this conditional master equation to the NPW representation [30], applying a series of approximations valid in the limit where N ≫ M (the details can be found in [31]), and finally transforming this representation into the NPW particle filter [83]. An overview of this process is presented in appendix A and in more detail in [32]. The final result is a swarm of K complex functions α^(k)(x) with corresponding weights w^(k), indexed with k∈(1,K), that obey the following differential equations⁶:

$$\begin{equation} \partial_t \alpha^{(k)}({{\bf x}}) = -\mathrm{i} \left(\frac{1}{\hbar}(h({{\bf x}},{\bf u}) + U|\alpha^{(k)}({{\bf x}})|^2) + \sum_i l_i({{\bf x}}) \zeta_i^{(k)}(t) \right) \alpha^{(k)}({{\bf x}}), \end{equation} \tag{ 15a }$$

$$\begin{equation} \partial_t w^{(k)} = 2 \sum_i \left( 2 L_i^{(k)} \mathbb{W}[L_i^{(\cdot)}] - (L_i^{(k)})^2 + 2 L_i^{(k)} \eta_i(t)\right)w^{(k)}, \end{equation} \tag{ 15b }$$

$$\begin{equation} L_i^{(k)} = \int {\rm d}{{\bf x}} \; l_i({{\bf x}}) \left(|\alpha^{(k)}({{\bf x}})|^2 - \frac{1}{2}\delta({{\bf x}})\right) \end{equation} \tag{ 16 }$$

and

$$\begin{equation} \mathbb{W}[f(\alpha^{(\cdot)}({{\bf x}}))] \equiv \frac{\sum_{k=1}^K w^{(k)} f(\alpha^{(k)}({{\bf x}}))}{\sum_{k=1}^K w^{(k)}} \end{equation} \tag{ 17 }$$

is the definition of the weighted average over the swarm. The Stratonovich noises η_i(t) are virtual measurement records (corresponding to separate repetitions of the 'experiment'), and are the same for all k. Importantly, no information is exchanged between different 'experiments', so each swarm can be integrated independently. This is not strictly true for each element of the swarm. However, the different k are only 'weakly' coupled via the weighted averages.

Expectations of observables are calculated by using the weighted average defined in equation (17). As specified in [30], each moment generated is related to a specific ordering of operators. In this paper, we only consider observables that can be generated from the one-body correlation matrix, calculated using

$$\begin{equation} \langle \hat {\psi }^\dagger({{\bf x}}) \hat {\psi }({\bf y}) \rangle = \mathbb{W}[\alpha^*({{\bf x}})\alpha({\bf y})] - \textstyle\frac{1}{2} \delta({{\bf x}}-{\bf y}). \end{equation} \tag{ 18 }$$

For comparison with equation (10), the expectation of the local number squared is given by

$$\begin{equation} \langle (\hat {\psi }^\dagger({{\bf x}}))^2 (\hat {\psi }({{\bf x}}))^2 \rangle = \mathbb{W}[|\alpha({{\bf x}})|^4] - 2\mathbb{W}[|\alpha({{\bf x}})|^2] \delta(0) + \textstyle\frac{1}{2} \delta(0)^2. \end{equation} \tag{ 19 }$$

In practice, the presence of the infinite quantity δ(0) causes no difficulty. For if the field α(x) is approximated as a discretized grid of points (which is required for numerical simulation), then δ(0) ≈ 1/Δx, where Δx is the spacing between points on the position-space grid. This feature is common to all phase-space simulations of quantum fields [37]. Unlike the Hartree–Fock method, expectations computed via equation (19) require weighted averages over the swarm. Hence, higher-order moments generated by the NPW particle filter are not locked to a mean field and can fluctuate dynamically. This freedom is what allows the NPW particle filter, and more generally phase-space methods, to investigate the dynamic quantum statistics of a system.

Simulating equation (1) via the NPW particle filter (15) requires appropriate initial conditions for the swarm of complex functions and the weights. Much like the TW method, these initial conditions are sampled from a quasi-probability distribution, and so in general are not deterministic. Appendix B gives the quasi-probability distribution from which to sample the initial conditions for a BEC, and further outlines an algorithm for this initial sampling. Analogous to the Hartree–Fock method, each sample has a quantized and fixed number of particles.

The evolution of the weights (equation (15b)) typically generates an exponential divergence between the largest and smallest weights. This can lead to poor sampling and thus an inaccurate estimate of the current quantum state of the BEC. There are a few numerical tricks to compensate for this issue. Firstly, it is more efficient to integrate the weights in log space, meaning w_l = ln(w) is integrated instead of directly integrating w. Secondly, it is important to frequently renormalize the weights to prevent them from becoming too small or too large. This is achieved by first calculating the current norm

$$\begin{equation} \bar{w} = \sum_{k=1}^K w^{(k)} \end{equation} \tag{ 20 }$$

and then reassigning weight values as follows:

$$\begin{equation} \frac{w^{(k)}}{\bar{w}} \rightarrow w^{(k)}, \end{equation} \tag{ 21 }$$

where we use → to mean 'write to' in an algorithmic sense. Immediately after reassigning weights, $\bar{w} = 1$. The final and most important numerical consideration is resampling. As the maximum weight $w _{\max} = \max _k w ^{(k)}$ becomes exponentially larger than the smallest weight $w _{\min } = \min _k w ^{(k)}$, the complex field with the smallest weight becomes negligible. This causes poorer sampling, and furthermore, the memory being used to store this complex field is wasted. We can therefore neglect the fields which have very small weights (compared to $w _{\max }$) and replace them with an appropriately rebalanced sample from the maximum weighted complex field. A simple algorithm for this process is presented below. We assume the user sets a tolerance _tol such that if the ratio of a weight to the largest weight is below _tol, then that weight is neglected.

• (i)  
  Find the index of the largest weight, $k_{w _{\max }}$, and store the weight's value

  $$\begin{eqnarray*} &&w^{(k_{w_{\max}})} = \max_k w^{(k)} \rightarrow w_{\max}. \end{eqnarray*}$$
• (ii)  
  Find the index of the smallest weight, $k_{w _{\min }}$, and store the weight's value

  $$\begin{eqnarray*} &&w^{(k_{w_{\min}})} = \min_k w^{(k)} \rightarrow w_{\min}. \end{eqnarray*}$$
• (iii)  
  If the minimum weight is relatively smaller than our tolerance, $w_{\min } / w_{\max } < \epsilon _{\mathrm { tol}}$, continue to (iv), otherwise skip to (vii).
• (iv)  
  Overwrite the amplitude possessing the minimum weight with the amplitude whose corresponding weight is the maximum weight:

  $$\begin{eqnarray*} &&\alpha^{(k_{w_{\max}})}({{\bf x}}) \rightarrow \alpha^{(k_{w_{\min}})}({{\bf x}}). \end{eqnarray*}$$
• (v)  
  Rebalance the weights: $w_{\max }/2 \rightarrow w^{(k_{w _{\min}})}$ and $w_{\max}/2 \rightarrow w^{(k_{w _{\max }})}$.
• (vi)  
  Go back to (i), repeat steps until all negligible weights have been removed.
• (vii)  
  Finish and return to the integration of equations (15).

This algorithm, with a few minor enhancements, was used in the simulations presented in this paper. Note that this algorithm is not the optimal solution to the resampling problem. Considerably more advanced and efficient resampling techniques have been developed, and can be found within the particle filtering literature. A tutorial by Arulampalam et al [45] gives an excellent introduction to these resampling techniques.

The numerical techniques required to simulate equation (1) using either the Hartree–Fock approximation or the NPW particle filter have been presented. We now apply and compare both of these methods to a condensate undergoing off-resonant imaging and active feedback damping. This comparison, and the use of the NPW particle filter to design an effective quantum-noise feedback control, are the key results of this paper.

We first consider active feedback damping of a BEC of two-level atoms undergoing continuous cavity-mediated measurement, as shown in figure 1. Damping via feedback to the trapping position has been studied extensively for neutral atoms and optomechanical resonators [16, 85–87], and the measurement has been demonstrated experimentally for optomechanical resonators [88] and ultracold atoms [89, 90]. Light is passed through an optical cavity containing the trapped BEC. We consider the regime of 'good measurement', where the light is far-detuned from the atomic resonance and the cavity mirrors are lossy (i.e. large cavity linewidth), thereby allowing the atomic excited state and cavity dynamics to be adiabatically eliminated, respectively [85]. Then the phase of the transmitted light contains a signal proportional to the overlap of the atomic density with the spatial envelope of the optical field in the cavity. A homodyne measurement of this phase can be used to continually update the estimate of the atomic state. We model this estimate of the system, conditioned on the measurement record, using a conditional master equation (also called a quantum filtering equation). Energy can be removed from the BEC by continuously adjusting the mean position of the trap by an amount proportional to the negative of the bulk momentum of the BEC (which can be estimated from the quantum filter). The equation of motion for the filter that governs the single-atom version of this system is found in [47, 85]. This single-atom conditional master equation is a good model for a weakly interacting BEC undergoing a continuous cavity-mediated measurement provided a collective mode of the atomic-field can be coupled to a single mode of the cavity [91]. Specifically, the Stratonovich quantum filtering equation for the conditional state of the atomic field is

$$\begin{equation} \partial_t \skew3\hat{\rho} = -\mathrm{i}[\skew3\hat{H}_{\mathrm{F}}(u_{\mathrm{s}}),\skew3\hat{\rho}] + \gamma \, \mathcal{D}[\skew3\hat{C}_\xi] \skew3\hat{\rho} + \gamma \, \mathcal{C}[\skew3\hat{C}_\xi] \skew3\hat{\rho} + \sqrt{\gamma} \mathcal{H}[\skew3\hat{C}_\xi] \skew3\hat{\rho} \; \eta(t). \end{equation} \tag{ 22 }$$

For notational convenience, we have written energy in units of ℏω, position in units of $\sqrt {\hbar /(m\omega )}$, and time in units of 1/ω, where ω is the frequency of the harmonic trapping potential and m is the mass of an individual atom. The unitary dynamics are governed by the atomic-field Hamiltonian

$$\begin{equation} \skew3\hat{H}_{\mathrm{F}} = \int \mathrm{d}x \; \hat {\psi }^\dagger(x) h_{\mathrm{F}}(x,u_{\mathrm{s}}) \hat {\psi }(x), \end{equation} \tag{ 23 }$$

where $\hat {\psi }(x)$ is the one-dimensional field operator satisfying $[\hat {\psi }(x), \hat {\psi }^\dagger (x')] = \delta (x - x')$, and the single-particle Hamiltonian

$$\begin{equation} h_{\mathrm{F}}(x,u_{\mathrm{s}}) = \frac{1}{2}\left(-\partial_x^2 + x^2\right) + u_{\mathrm{s}} \frac{\langle \skew3\hat{P} \rangle}{\langle \skew3\hat{N} \rangle}x \end{equation} \tag{ 24 }$$

contains the kinetic energy, potential energy due to the harmonic trap and the linear feedback of strength u_s > 0. $\hat {P}$ and $\hat {N}$ are the many-body momentum and number operators, respectively, defined as

$$\begin{equation} \skew3\hat{P} = \int \mathrm{d}x\; \hat {\psi }^\dagger(x) (-\mathrm{i} \partial_x) \hat {\psi }(x), \end{equation} \tag{ 25 }$$

$$\begin{equation} \skew3\hat{N} = \int \mathrm{d}x \; \hat {\psi }^\dagger(x) \hat {\psi }(x). \end{equation} \tag{ 26 }$$

The strength of the measurement is given by γ, which is proportional to the effective atom–cavity coupling rate and inversely proportional to the cavity linewidth. Increasing γ gathers more information about the system per unit time, however it also increases the rate of heating due to the measurement backaction. The observable that is measured is given by the measurement operator

$$\begin{equation} \skew3\hat{C}_\xi = \int \mathrm{d}x \; \hat {\psi }^\dagger(x) c_{\xi}(x) \hat {\psi }(x), \end{equation} \tag{ 27 }$$

where c_ξ(x) = cos²(ξx − π/4) is proportional to the intensity of the intracavity optical field. The dimensionless length ξ = 2πx_HO/λ is the ratio of the natural length scale of the trap, $x_{\mathrm { HO}} = \sqrt {\hbar /(m\omega )}$, and the wavelength λ of the optical cavity [85].

In the limit ξ ≪ 1, c_ξ(x) ≈ 1/2 + ξx, and so the cavity-mediated measurement is approximately a measurement of the centre-of-mass position of the condensate. A BEC undergoing a continuous position measurement was used as a testing ground for the NPW simulation method due to the existence of reliable benchmarks [31]. However, in contemporary experiments a very small ξ can be difficult to achieve, and so we consider the more moderate values of ξ = 0.1 and 0.5.

We integrated equation (22) using both the restrictive Hartree–Fock approximation and the more complete NPW particle filter. Since equation (22) is a special case of equation (1), we can use the general result of equation (11) to find the equation of motion for the unnormalized macroscopic 'wavefunction' $\tilde {\phi }(x,t)$:

$$\begin{equation} \fl \partial_t \tilde{\phi}(x,t) = \Big\{-\mathrm{i} h_{\mathrm{F}}(x, u_{\mathrm{s}}) + \gamma \left(2 c_\xi(x) C_\xi^\phi(t) -c_\xi(x)^2 \right) + \sqrt{\gamma} \; c_\xi(x) \; \eta(t) \Big\} \tilde{\phi}(x,t), \end{equation} \tag{ 28 }$$

where

$$\begin{equation} C_\xi^\phi(t) = \int \mathrm{d}x\; c_{\xi}(x) \frac{|\tilde{\phi}(x,t)|^2}{n_\phi(t)}, \end{equation} \tag{ 29 }$$

$$\begin{equation} n_\phi(t)= \int \mathrm{d}x\; |\tilde{\phi}(x,t)|^2. \end{equation} \tag{ 30 }$$

All the atoms have the same wavefunction $\tilde {\phi }(x)$ under the Hartree–Fock approximation. Similarly, the NPW particle filter for equation (22) is given by the general result of equations (15):

$$\begin{equation} \partial_t \alpha^{(k)}(x) = -\mathrm{i}\Big(h_{\mathrm{F}}(x,u_{\mathrm{s}}) + \sqrt{\gamma} c_\xi(x) \; \zeta^{(k)}(t)\Big) \alpha^{(k)}(x), \end{equation} \tag{ 31a }$$

$$\begin{equation} \partial_t w^{(k)} = 2 \left\{\gamma \left(2C_\xi^{(k)} \mathbb{W}[C_\xi^{(\cdot)}] - (C^{(k)}_\xi)^2 \right) + \sqrt{\gamma} C_\xi^{(k)}\; \eta(t) \right\}w^{(k)}, \end{equation} \tag{ 31b }$$

$C_\xi ^{(k)} = \int \mathrm {d}x\, c_\xi (x) \left (|\alpha ^{(k)}(x)|^2 - \delta (0)/2\right )$

A comparison of the Hartree–Fock method and NPW particle filter simulations⁷ of a feedback-controlled BEC undergoing a continuous cavity-mediated measurement are shown in figure 2(a). For both the Hartree–Fock wavefunction and NPW particle filter, the initial state was a BEC with a constant phase and amplitude, given by the displaced Gaussian

$$\begin{equation} \alpha_0(x)= \sqrt{\frac{N}{\sqrt{2 \pi} \sigma}}\exp\left(-\frac{(x-x_0)^2}{4\sigma^2}\right), \end{equation} \tag{ 32 }$$

where x₀ is the Gaussian's displacement from the origin, σ is its variance and N is the average total number of the condensate. More details on the choice of initial condition, a technique for sampling this initial state, and a discussion on how it affects the convergence of the numerics can be found in appendix B. When the cosine measurement is close to a position measurement, ξ = 0.1, we see the NPW particle filter and Hartree–Fock simulations match well, with both predicting damping of the BEC to a steady state. For larger values of ξ, the cosine measurement is no longer 'position-like', in that it has some curved structure on the length scale of the BEC. When ξ = 0.5, we see that the NPW particle filter predicts significantly more disruption to the BEC than the Hartree–Fock solution. This disruption, due to the measurement backaction, causes additional heating that cannot be counteracted by simple linear feedback damping. Consequently, if a feedback-control experiment was designed solely on the basis of the Hartree–Fock simulation, then this experiment would likely fail, for it would not account for this additional backaction effect neglected by the Hartree–Fock simulation.

Zoom In Zoom Out Reset image size

The Hartree–Fock method misses the additional measurement backaction effect on the BEC because it neglects certain spontaneous-emission events. Recall that the derivation of conditional master equation (22) assumes that the light interacting with the two-level atoms is sufficiently far-detuned from the atomic resonance that the atomic excited state can be adiabatically eliminated. Thus, all spontaneous-emission events in the system are now encoded as scattering events between the light and BEC. These spontaneous-emission events change the phase of the outgoing probe beam, which is how the cosine-squared moment of the BEC is measured. As the Hartree–Fock approximation assumes all atoms are in the same single-particle state, it restricts which spontaneous-emission events can be described. In contrast, the NPW particle filter allows non-collective states of the atoms, and can therefore simulate all spontaneous-emission events. It is for this reason that we call this additional backaction effect spontaneous-emission noise. Note also that spontaneous-emission noise arises due to non-classical correlations in the atomic field. Hence, in some sense, it is also a form of quantum noise.

Although the effects of spontaneous-emission noise cannot be removed via simple linear damping, they can be cancelled with a more exotic choice of feedback control. We present this control, and demonstrate its effectiveness, in the next section.

The Hartree–Fock simulation converged to the NPW particle filter result when ξ = 0.1, which shows that the net effect of the measurement backaction is partially mode-dependent. Whether the BEC cools to a steady state depends on the competition between the heating due to the scattering and the damping due to the feedback. The linear feedback used in the simulations of figure 2(a) damped the centre-of-mass motion of the BEC, which for ξ = 0.1 is very close to the spatial mode of the disruption due to measurement backaction. When ξ = 0.5, energy was added to modes other than the position moment mode, so we hypothesized that the solution was to add an extra control that targets the cosine-squared mode of the BEC.

A methodology used to damp specific modes of a BEC is described in [29]. We can apply this methodology to create a control that targets the cosine-squared, c_ξ(x), mode of the BEC. The single-particle Hamiltonian that includes this new quantum-noise control is

$$\begin{equation} h_{\mathrm{FC}}(x,u_{\mathrm{s}},u_{\mathrm{c}}) = h_{\mathrm{F}}(x,u_{\mathrm{s}}) + \frac{u_{\mathrm{c}} c_\xi(x)}{\langle \skew3\hat{N} \rangle} \int \mathrm{d}y \; c_\xi'(y) \mbox{Im}\left[\langle \hat {\psi }^\dagger(y) \hat {\psi }'(y) \rangle\right] , \end{equation} \tag{ 33 }$$

where u_c is the strength of the quantum-noise control for the cavity-mediated measurement, $\hat {\psi }'(x) \equiv \partial _x \hat {\psi }(x)$, and c_ξ'(x) ≡ ∂_xc_ξ(x).

NPW particle filter simulations of a BEC undergoing continuous cavity-mediated measurement without and with this novel quantum-noise control are shown in figure 2(b). In part (b, i) we can see that the disruption due to the spontaneous-emission noise is completely cancelled by the quantum-noise control. Comparing parts (b, ii) and (b, iii), we see that the quantum-noise control vastly improves the stability of the BEC under continuous monitoring. In (b, ii) we see that a continuous cavity-mediated measurement causes the BEC to break apart, even in the presence of feedback. In contrast, (b, iii) shows that the additional quantum-noise control cancels the spontaneous-emission noise, thereby preventing this spread.

We have demonstrated that spontaneous-emission noise can be cancelled in a weakly interacting BEC undergoing continuous cavity-mediated measurement, thus showing that a BEC can be feedback-cooled to a steady state close to the ground-state energy. However, the conditional master equation (22) that describes the cavity-mediated measurement is only valid when inter-atomic interactions in the condensate are negligible. Typically, inter-atomic interactions make the BEC larger than the optical wavelength, bringing it out of the near-position measurement regime. Indeed, the coupling between the cavity mode and collective position mode requires negligible inter-atomic interactions. Although a non-interacting condensate can be created with a dilute atomic sample or via a Feshbach resonance [1], it is technically demanding. Furthermore, there are instances where large inter-atomic interactions are desirable, such as for stable atom laser operation [93, 94] and for the generation of non-classical atom laser states [42]. Fortunately, feedback control is still possible in this regime using a phase-contrast imaging setup.

We consider the feedback control of a BEC with inter-atomic interactions of arbitrary strength U under phase-contrast imaging, which has been realized experimentally and has allowed for the real-time continuous measurement of a BEC [18, 95]. In phase-contrast imaging, an off-resonant laser beam passes through the BEC, gaining a phase shift that is proportional to the atomic column density (see figure 3). This light is measured by an array of CCD cameras, which we model as an array of homodyne detectors. The conditional master equation for a quasi-one-dimensional 'cigar-shaped' BEC under phase-contrast imaging is derived in [14], and in Stratonovich form is

$$\begin{eqnarray}&&\fl \partial_t \skew3\hat{\rho} = -\mathrm{i} [\skew3\hat{H}_{\mathrm{FN}}(u_{\mathrm{s}}),\skew3\hat{\rho}] + \gamma \int {\rm d}x \,\mathcal{D}[\skew3\hat{M}_{\nu}(x)] \skew3\hat{\rho} \nonumber \\ &&+ \gamma \int {\rm d}x \, \mathcal{C}[\skew3\hat{M}_{\nu}(x)] \skew3\hat{\rho} + \sqrt{\gamma} \int {\rm d}x \, \mathcal{H}[\skew3\hat{M}_{\nu}(x)] \skew3\hat{\rho} \; \eta(x,t), \end{eqnarray} \tag{ 34 }$$

where

$$\begin{equation} \skew3\hat{H}_{\mathrm{FN}}=\skew3\hat{H}_{\mathrm{F}}+\frac{U}{2}\int {\rm d}x \;\hat {\psi }^\dagger(x)\hat {\psi }^\dagger(x)\hat {\psi }(x)\hat {\psi }(x) \end{equation} \tag{ 35 }$$

is the Hamiltonian for an interacting BEC under linear feedback (see equation (2)), and the measurement operators are

$$\begin{equation} \skew3\hat{M}_\nu(x)=(\mu_\nu * \hat {\psi }^\dagger \hat {\psi })(x). \end{equation} \tag{ 36 }$$

These operators describe a measurement of the local particle density operator, $\hat {\psi }^\dagger (x) \hat {\psi }(x)$, convolved with a kernel function μ_ν(x), at every position x. Here we use * to indicate a convolution, which is defined as

$$\begin{equation} (f*g)(x) \equiv \int_{-\infty}^{\infty}{\rm d}y f(y) g(x-y). \end{equation} \tag{ 37 }$$

The kernel function μ_ν(x) depends on the dimensionless resolution length scale ν = x_⊥/k₀x²_HO, where x_⊥ is the size of the condensate in the tight trapping directions, $x_{\mathrm { HO}} = \sqrt {\hbar /(m \omega )}$ for loose trapping frequency ω and k₀ = 2π/λ is the wavenumber of the off-resonant light. Typically, ν ≪ 1. In the regime where light is predominantly scattered in the direction of the laser beam, the kernel function in Fourier space is given by $\tilde {\mu }_\nu (k) = \exp (-\nu k^4)$. Since there are multiple measurement channels, indexed by the position x, there are also multiple Stratonovich noises, η(x,t). Note that, like equations (22), (34) has been written in harmonic oscillator units, where energy is in units of ℏω, position is in units of x_HO, and time is in units of 1/ω.

Zoom In Zoom Out Reset image size

In this system, the measurement signal provides a resolution-limited density measurement of the BEC. Unlike the cavity-mediated measurement, in no limit does the signal provide a simple position measurement, although multiple position moments can be estimated from the image, whose resolution is defined by the parameter ν. We investigate two values for the resolution length scale: ν = 10 and 0.1, which provide coarse-resolution and fine-resolution measurements of the BEC density, respectively.

Under the Hartree–Fock approximation, equation (34) reduces to [15]

$$\begin{eqnarray} &&\fl \partial_t \tilde{\phi}(x,t) = \Bigg\{-\mathrm{i} \Big(h_{\mathrm{F}}(x,u_{\mathrm{s}}) + U(N-1) \frac{| \tilde{\phi}(x,t)|^2}{n_\phi(t)}\Big) \nonumber\\ &&+ 2 \frac{\gamma}{n_\phi(t)} (\mu_{\nu} * \mu_{\nu} * |\tilde{\phi}|^2)(x,t) + \sqrt{\gamma} (\mu_{\nu} * \eta)(x,t) \Bigg\} \tilde{\phi}(x,t). \end{eqnarray} \tag{ 38 }$$

The NPW particle filter converts equation (34) to [31, 32, 83]:

$$\begin{equation} \fl \partial_t \alpha^{(k)}(x) = -\mathrm{i} \Big( h_{\mathrm{F}}(x,u_{\mathrm{s}}) + U |\alpha^{(k)}(x)|^2 + \sqrt{\gamma} (\mu_{\nu} * \zeta^f)(x) \Big) \alpha^{(k)}(x), \end{equation} \tag{ 39a }$$

$$\begin{equation} \fl \partial_t w^{(k)} = 2 \int {\rm d}x ( \gamma ( 2 M_\nu^{(k)}(x) \mathbb{W}[M_\nu^{(\cdot)}(x)] - M_\nu^{(k)}(x)^2 ) + \sqrt{\gamma} M_\nu^{(k)}(x,t)\; \eta(x) )w^{(k)}, \end{equation} \tag{ 39b }$$

A comparison of Hartree–Fock method and NPW particle filter simulations of equation (34) is shown in figure 4(a). This analysis is performed without a quantum-noise control. When ν = 10 we see that the Hartree–Fock method agrees with the NPW particle filter, and both predict cooling of the BEC to a steady state. This shows that for coarse-resolution density measurements, higher-order modes are weakly excited such that feedback to the position mode provides sufficient damping, and the energy of the BEC stays low. Unfortunately, reaching this parameter limit can be experimentally difficult [15]. When the BEC is probed at a finer spatial resolution (ν = 0.1), the Hartree–Fock solution diverges from the NPW particle filter solution, which shows that the spontaneous-emission noise excites higher-order modes of the BEC that are not cooled by the linear feedback control.

Zoom In Zoom Out Reset image size

As with the cavity-mediated measurement, we can introduce controls for higher-order modes that cancel the heating caused by spontaneous-emission noise. Since our model of phase-contrast imaging gives a multi-channel measurement, we require multiple channels of feedback, effected via the trapping potential. The single-particle Hamiltonian h_FM(x,u_s,u_μ) describing this new feedback is

$$\begin{equation} h_{\mathrm{FM}} (x,u_{\mathrm{s}},u_\mu) = h_{\mathrm{F}}(x,u_{\mathrm{s}}) + \frac{u_\mu}{\langle \skew3\hat{N} \rangle}\left(\mu_\nu * \mu_\nu' * \mbox{Im}[\langle \hat {\psi }^\dagger \hat {\psi }' \rangle]\right)(x) , \end{equation} \tag{ 40 }$$

where μ_ν'(x) = ∂_xμ_ν(x) and u_μ is the strength of the quantum-noise control for the direct measurement. Note that this feedback requires a high degree of spatial and temporal control of the BEC's trapping potential. Fortunately, there are several proposals for implementing such potentials experimentally [96–100].

The effectiveness of the quantum-noise control defined in equation (40) is shown in figure 4(b). In (b, i) we can see that the quantum-noise control completely cancels the disruption caused by a strong phase-contrast measurement of the BEC. In parts (b, ii) and (b, iii) we compare the density of the BEC without and with the quantum-noise control, respectively. In (b, ii) we see that the spontaneous-emission noise causes the BEC to spread rapidly. In contrast, (b, iii) shows that the quantum-noise control (40) prevents the density excitations that cause the condensate to break apart.

We have shown that feedback control can be used to stabilize a continuously monitored BEC close to the ground-state energy for both measurement apparatuses shown in figures 1 and 3, and for a wide range of parameter regimes. Experimentally, it is simplest to implement only linear feedback control via adjustments to the trapping minimum⁸. It might therefore be tempting to try and feedback-cool a BEC in the 'weak-probing' regime (e.g. ξ = 0.1 for the cavity-mediated measurement and ν = 10 for phase-contrast imaging), since spontaneous-emission noise is negligible and only the bulk feedback mechanism, governed by parameter u_s, is required to achieve net cooling to steady state. However, operating in this parameter regime requires relatively small, tightly trapped condensates [14, 85], although this may change with the development of trapping technologies. If the target condensate does not satisfy these requirements, more complex measurement and feedback is required. For moderately sized condensates, a cavity-mediated measurement with an extra quantum-noise control appears promising. As demonstrated in section 3.2, the shape of the quantum-noise control must match the measurement. Thus, in the case of a continuous cavity-mediated measurement, it might be possible to implement the quantum-noise control via the probe beam itself. For very large, loosely trapped condensates, it is likely phase-contrast imaging will be the only option, in which case full spatial control of the trapping potential will be required.