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

Off-resonant optical imaging is the most popular method for continuous monitoring of a Bose-Einstein condensate (BEC). However, the disturbance caused by scattered photons places a serious limitation on the lifetime of such continuously-monitored condensates. In this paper, we demonstrate that a new choice of feedback control can overcome the heating effects of the measurement backaction. In particular, we show that the measurement backaction caused by off-resonant optical imaging is a multimode quantum-field effect, as the entire heating process is not seen in single-particle or mean-field models of the system. Correctly simulating such continuously-monitored systems is only possible using the number-phase Wigner (NPW) particle filter, which is a hybrid between the leading techniques for simulating non-equilibrium dynamics in condensates and particle filters for simulating high-dimensional non-Gaussian filters in the field of engineering. The new control scheme will enable long-term continuous measurement and feedback on one of the leading platforms for precision measurement and the simulation of quantum fields, allowing for the possibility of single-shot experiments, adaptive measurements and robust state-preparation and manipulation.


Introduction
A Bose-Einstein condensate (BEC) in a dilute atomic gas can be well-isolated from its environment, and is highly controllable using a combination of optical, rf and magnetic fields [1]. This gives it the potential to address a broad range of research questions in fundamental and applied science, including studies in quantum nonequilibrium thermodynamics [2,3], entanglement of massive particles [4] and the quantum simulation of both cosmological phenomena, such as Hawking radiation emitted from a black hole's event horizon [5], and phase transitions in condensed matter, including superconductivity and quantum magnetism [6]. Furthermore, Bose-condensed sources are likely to be key components in a range of future technologies, such as improved precision inertial sensors based on atom interferometry, where the sensitivities of current devices are limited by the properties of thermal sources [7,8,9]. This wealth of research opportunity has therefore led to great practical interest in the ability to control the spatial state of the BEC's quantum field. Research has predominantly focussed on open-loop control of the condensate by direct control of optical and magnetic potentials. In contrast, continuous measurement feedback control of BECs is rarely employed due to their perceived fragility under continuous measurement. However, this intuition does not account for the advantages active feedback control can provide. Active feedback can control properties of the BEC using the information gained from the continuous measurement, which includes properties that would have been adversely affected by the backaction. Furthermore, active feedback provides robust and reliable behaviour that cannot be matched by open-loop control schemes [10]. In particular, it has been shown that the linewidth of an outcoupled atom laser can be improved by continuously monitoring the BEC with off-resonant light and applying active feedback [11,12,13,14,15,16]. However, these prior analyses used single-particle or mean-field models of the system, and thus did not fully incorporate the multi-mode quantumfield effects of the measurement backaction. In this paper, we present the first analysis of a feedback-controlled BEC that does include the crucial effects of measurementinduced spontaneous-emission noise, by using the number-phase Wigner (NPW) particle filter to perform our simulations. Furthermore, we demonstrate that the effects of the spontaneous-emission noise can be cancelled with the application of an active feedback control that is tailored to the measurement being applied. This result demonstrates that there is no fundamental limit to the lifetime of a condensate due to off-resonant imaging, and also shows that improving the linewidth of outcoupled atom lasers from condensates is feasible with active feedback.
Continuous off-resonant measurement of BEC has been achieved in experiment [17,18,19], and potentially has numerous and varied applications. Recent proposals that rely upon the continuous off-resonant imaging of ultracold atomic systems include the generation of entanglement in hybrid optomechanical systems [20], the generation of coherence between two initially uncorrelated condensates [21], an investigation into the emergence of classical dynamics from a quantum system [22], and the simulation of condensed matter systems [23,24]. Unfortunately, in previous experiments the decoherence due to the spontaneous scattering of the photons placed hard limits on the lifetime of the condensate [18,25]. Such disruption could be mitigated by placing a BEC in a cavity [26,27,28], but it is was previously unknown if this disruption could be cancelled entirely or if continuous monitoring places a fundamental limit on the lifetime of the BEC. Early work by Haine et al. [29] showed that density disruptions in a condensate could be reduced using a feedback control based on simple adjustments to the BEC's confining potential, albeit in the unrealistic situation where measurement backaction was neglected. The effect of measurement backaction was incorporated in Wilson et al. [13], which studied the effectiveness of simple damping feedback with a single-mode model of a BEC undergoing continuous centre-of-mass position measurement. Most relevant to this paper, previous analyses by Szigeti et al. [14,15], who considered the semiclassical limit of a multi-mode quantum filtering model for an off-resonantly imaged, feedback-cooled BEC, suggested that active feedback could remove the disruption caused by off-resonant imaging entirely. However, these models were analysed under the Hartree-Fock approximation, which assumes the condensate has a fixed number of atoms all in the same single-particle state, which in turn neglects part of the spontaneous-emission noise, making this analysis incomplete. Using the NPW particle filter to perform a multi-mode quantum-field analysis, we avoided such approximations and discovered that an additional control mechanism is required.
The NPW particle filter allows for the simulation of a feedback controlled, continuously-monitored BEC that includes the quantum statistics ignored by the Hartree-Fock method. It has been developed over a sequence of papers by Hush et al., and is currently the most precise numerical method for simulating the quantum statistics of monitored BECs. The NPW particle filter method has two main components. The first is the Number-Phase-Wigner representation [30,31,32] which falls into a class of stochastic methods based on phase-space representations, which includes both truncated Wigner (TW) [33,34,35,36,37], positive-P (P + ) [38,39,35] and other variants [40,41]. These methods have proven to be highly efficient at investigating the dynamics of BECs that cannot be probed by mean-field methods, including the generation of entanglement [4], squeezing [42] and thermal fluctuations [34]. The second aspect of the NPW particle filter method is an application of an engineering technique called particle filters [43,44,45,46]. Particle filters are the most efficient control theory methods available to continually estimate systems that are both high dimensional and nonlinear [45], such as a monitored BEC.
This paper aims to provide a complete investigation into the feedback controls that are required to remove the effects of spontaneous emission due to off-resonant imaging of a BEC. In doing so, it also provides the first practical application of the NPW particle filter. Previous work has verified the accuracy of this simulation technique and compared its performance to other established methods [31]. To make this paper as accessible as possible, we provide an overview of the relevant numerical techniques, including the NPW particle filter, in section 2. Section 3 investigates the heating induced by spontaneous-emission noise due to off-resonant optical imaging, and demonstrates how a novel control scheme can control this heating effect, leading to net cooling of the BEC to steady state.
In section 2.1 we introduce the conditional master equation for a BEC under continuous off-resonant imaging. In anticipation of future applications, we have kept the form of the conditional master equation as generic as possible. Section 2.2 introduces the Hartree-Fock approximation, which is the relevant semiclassical approximation for continuous off-resonant optical imaging of a BEC. Finally, in section 2.3 we give an overview of the NPW particle filter, including the algorithmic details for its implementation.

Conditional master equation
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 [57,58,50], which is better known as a filter in the engineering community [59,60,54]. The conditional master equation of a BEC undergoing continuous off-resonant imaging is governed by the Stratonovich equation [13,61,14,15,16,32]: whereρ 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 whereψ(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 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ρ 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 [62]. The effect of the measurement is described by the remaining terms in the conditional master equation (1). The superoperators are the decoherence, innovations, and Stratonovich correction superoperators, respectively, for an arbitrary operatorĉ where · = Tr[·ρ]. The decoherence superoperator describes the effect of the measurement backaction on the system, while the innovations superoperator incorporates the information obtained from the measurement by continuously updating the density matrix. The specifics of the backaction and measurement signal depend upon the set of measurement operatorŝ 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. [61], Szigeti et al. [14,15], and others [12,63,64,65,66]. 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 [53]. Explicitly, if the measurement signal changes by ∆Y i within a time interval ∆t, then 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 [67,49,54,50]. 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 stochastic differential equations (SDEs) (for some examples, see [68]). 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,ρ (p) , here indexed by p) and then average these paths (or trajectories) to get a result. We use the notation A[·] to mean this procedure. For example, in order to find the average density over time, we simulate equation (1) P times and then compute 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ρ (p) can be thought of as a new 'experiment'. We note that if there was no feedback then A[ρ] =ˆ in the limit as P → ∞, whereˆ 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ψ(x) = M m=1 u m (x)â m , where {u m (x)} m are a set of orthonormal functions, and each modeâ m has a set of J energy eigenstates labeled |j m . The density matrix is thenρ 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.

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 whereρ 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 and 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: whereφ(x, t) is related to the normalized wavefunction by and we have defined 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 meanfield 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.

NPW particle filter
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 truncated Wigner (TW) phase-space method [33,34,35,36,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 [69]. 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,70], 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 [71,72,73]. Concepts similar to weights have also been used in phasespace methods [74,75], quantum measurement theory [76] and quantum control [77]. Of the many applications, the techniques that are closest in spirit to our approach are so-called particle filters [43,44,45,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 [78,79,80], GPS navigation [81] and financial modelling [82,83]. Particle-filter-like simulations appear to have first been applied to quantum systems in the context of simulating monitored BECs [84]. However, very soon after this publication they were also applied by Jacobs [85] 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 [84]. 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 ‡: For notational compactness, we suppress the time index for the functions α (k) and weights w (k) .
is a 'fictitious' Stratonovich noise increment, which is generated independently for each of the K paths, and 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: For comparison with equation (10), the expectation of the local number squared is given by 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 positionspace 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: and then reassigning weight values as follows: where we use → to mean 'write to' in an algorithmic sense. Immediately after reassigning weights,ω = 1. The final and most important numerical consideration is resampling. As the maximum weight ω max = max k ω (k) becomes exponentially larger than the smallest weight ω min = min k ω (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 ω 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 ωmax , and store the weight's value: (ii) Find the index of the smallest weight, k ω min , and store the weight's value: (iii) If the minimum weight is relatively smaller than our tolerance, w min /w max < 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: (v) Rebalance the weights: w max /2 → w (kω min ) and w max /2 → w (kω 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.

Cancelling spontaneous-emission noise
We consider the continuous monitoring of an harmonically trapped BEC under two common off-resonant imaging systems: a cavity-mediated measurement and phasecontrast imaging, both shown schematically in figure 1 and figure 3. Our aim is to determine whether active feedback control can be used to reduce the excitations in the condensate density caused by the measurement backaction (or, indeed, excitations due to any other process). In particular, we examine the effectiveness of a linear control that alters the trapping position [29,14,15,8] using both the semiclassical Hartree-Fock approximation and the NPW particle filter. We show that in both imaging apparatuses there exist experimentally accessible regimes where the measurement induces additional heating in the condensate that cannot be removed by simple feedback to the trapping position. This additional heating is not modelled by the Hartree-Fock simulations, and thus is a quantum multi-mode effect which we call spontaneous-emission noise (for reasons that become clear shortly). Fortunately, as we also show in section 3.2, it is possible to design a control that cancels the spontaneous-emission noise, allowing the continuously-monitored BEC to be cooled to a steady state.

Cavity-mediated measurement
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 [86,87,88,16], and the measurement has been demonstrated experimentally for optomechanical resonators [89] and ultracold atoms [90,91]. 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 [86]. 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 BEC Standing Wave Potential Control Coherent Light Source Figure 1. Schematic diagram of a BEC under feedback control using a continuous cavity-mediated measurement. The condensate is coupled to the mode of a cavity (depicted by the red cosine-squared surface). Then a homodyne measurement of light output from the cavity gives a measurement signal proportional to a density moment of the atomic cloud. The green coils represent controls that modify the trapping field with the aim of driving the BEC to a low energy steady state, the trapping field itself is not shown.
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 [86,49]. 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 [92]. Specifically, the Stratonovich quantum filtering equation for the conditional state of the atomic field is: whereψ(x) is the one-dimensional field operator satisfying [ψ(x),ψ † (x )] = δ(x − x ), and the single-particle Hamiltonian contains the kinetic energy, potential energy due to the harmonic trap, and the linear feedback of strength u s > 0.P andN are the many-body momentum and number operators, respectively, defined aŝ 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 where c ξ (x) = cos 2 (ξ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 HO = /(mω), and the wavelength λ of the optical cavity [86]. 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'φ(x, t): where All the atoms have the same wavefunctionφ(x) under the Hartree-Fock approximation. Similarly, the NPW particle filter for equation (22) is given by the general result of equations (15): where C Observables were calculated from equations (31) using the weighted average defined in equation (17) and the correspondences (18).
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 where x 0 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. 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 spontaneousemission events can be described. In contrast, the NPW particle filter allows noncollective states of the atoms, and can therefore simulate all spontaneous-emission  Figure 2. (a) Comparison of the Hartree-Fock method and the NPW particle filter for the simulation of a feedback-controlled BEC undergoing a continuous cavity-mediated measurement. The NPW particle filter was integrated using equations (31), averaged over P = 10 paths, and is plotted in red with circles. Hartree-Fock simulations were performed by integrating equation (28), averaging over P = 100 paths, and are plotted in purple with stars. A small value of ξ = 0.1 is plotted in (a,i), while a higher strength of ξ = 0.5 is plotted in (a,ii). Both plots have γ = 5, N = 100 and u s = 1. The Hartree-Fock method only agrees with the NPW particle filter when ξ = 0.1 [see inset of (a,i)]; it fails to predict the additional disruption to the BEC, caused by spontaneous-emission noise, when ξ = 0.5. (b) A comparison of a controlled BEC under cavity-mediated measurement without and with the quantum-noise control, averaged over P = 10 paths. (b,i) shows a comparison of the per-particle energy of the system, while the density fluctuations in the system without and with control are shown in (b,ii) and (b,iii), respectively. In (b,i): integration of the NPW particle filter (31) without quantum-noise control, u c = 0, is plotted in red with circles; integration with the addition of the quantum-noise control, u c = 5, is plotted in orange with diamonds. (b,ii) and (b,iii) are contour plots of the BEC density over an individual path, with (light) blue hues to (dark) red hues indicating low and high density, respectively. The noise in the contour plots is a physical consequence of the measurement process, rather than numerical uncertainty (which is less than 0.1%). All integrations were performed with u s = 1, ξ = 0.5, γ = 5 and N = 100.
events. It is for this reason that we call this additional backaction effect spontaneousemission 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.

Quantum-noise control
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 modedependent. 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 where u c is the strength of the quantum-noise control for the cavity-mediated measurement,ψ (x) ≡ ∂ xψ (x), and c ξ (x) ≡ ∂ x c ξ (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 spontaneousemission 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 groundstate energy. However, the conditional master equation (22) that describes the cavitymediated 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 interatomic 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 [94,95] and for the generation of non-classical atom laser states BEC Homodyne Detector Array Coherent Light Source Potential Control Array Figure 3. Schematic diagrams of a BEC under feedback control using continuous phase-contrast imaging. The condensate is illuminated by off-resonant light. After interacting with the atoms, the light is detected by an array of CCD cameras (which we can model as an array of homodyne measurements). The green wires represent controls that modify the trapping field with the aim of driving the BEC to a low energy steady state, the trapping field itself is not shown. [42]. Fortunately, feedback control is still possible in this regime using a phase-contrast imaging setup.

Phase-contrast imaging
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,96]. In phasecontrast 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: is the Hamiltonian for an interacting BEC under linear feedback [see equation (2)], and the measurement operators arê These operators describe a measurement of the local particle density operator, ψ † (x)ψ(x), convolved with a kernel function µ ν (x), at every position x. Here we use * to indicate a convolution, which is defined as The kernel function µ ν (x) depends on the dimensionless resolution length scale ν = x ⊥ /k 0 x 2 HO , where x ⊥ is the size of the condensate in the tight trapping directions, x HO = /(mω) for loose trapping frequency ω, and k 0 = 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μ ν (k) = exp (−νk 4 ). Since there are multiple measurement channels, indexed by the position x, there are also multiple Stratonovich noises, η(x, t). Note that, like equation (22), equation (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/ω.
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]: The NPW particle filter converts equation (34) to [84,31,32]: ν (x) = (µ ν * |α (k) | 2 )(x). 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 quantumnoise 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.
As with the cavity-mediated measurement, we can introduce controls for higherorder 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 F M (x, u s , u µ ) describing this new feedback is: 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 [97,98,99,100,101]. 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.

Discussion of Experimental Feasibility
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 figure 1 and figure 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 [86,14], although this may change with the development of  Figure 4. (a) Comparison of the Hartree-Fock method and the NPW particle filter for the simulation of a feedback-controlled BEC undergoing continuous phase-contrast imaging. The NPW particle filter was integrated using equations (39), averaged over P = 10 paths, and is plotted in red with circles. Hartree-Fock simulations were performed by integrating equation (38), averaging over P = 100 paths, and are plotted in purple with stars. A coarse-resolution measurement, ν = 10, is plotted in (a,i), while a fine-resolution measurement, ν = 0.1, is plotted in (a,ii). Both plots have γ = 1, N = 100, U/N = 3 and u s = 1. The Hartree-Fock method only agrees with the NPW particle filter when ν = 10 [see inset of (a,i)]; it fails to predict the additional disruption to the BEC, caused by spontaneous-emission noise, when ν = 0.1. (b) A comparison of a controlled BEC under continuous phase-contrast imaging without and with the quantum-noise control averaged over P = 10 paths. (b,i) shows a comparison of the per-particle energy of the system, while the density fluctuations in the system without and with control are shown in (b,ii) and (b,iii), respectively. In (b,i), the integration of the NPW particle filter (39) without quantum-noise control, u µ = 0, is plotted with red circles. A simulation of equations (39) with the addition of the quantum-noise control, u µ = 5, is plotted in orange with diamonds. (b,ii) and (b,iii) are contour plots of the density of the BEC over an individual path, with (light) blue hues to (dark) red hues indicating low to high densities, respectively. The noise in the contour plots is a physical consequence of the measurement process, rather than numerical uncertainty (which is less than 0.1%). All integrations were performed with u s = 1, ν = 0.1, U/N = 3, γ = 1 and N = 100.
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.

Conclusion
This paper has investigated the previously unexplored multi-mode quantum-field dynamics of a feedback-cooled multi-mode BEC undergoing (a) continuous cavitymediated cosine-squared measurement and (b) continuous phase-contrast imaging. Extended quantum-field simulations, performed with the recently developed NPW particle filter, have revealed regimes with additional measurement-induced disruption to the condensate that are not included in single-mode and semiclassical models, and furthermore cannot be counteracted with simple linear and/or quadratic feedback controls. Fortunately, we have developed a more sophisticated quantum-noise control that cancels the effect of this additional measurement backaction, and shown that this allows for successful feedback-cooling of the condensate to a steady state. This has been an important demonstration of the necessity of multi-mode quantum-field simulations, and in particular the NPW particle filter, to the successful design and implementation of measurement-based feedback-control schemes for multi-mode quantum systems. methods between the WSDEs and particle filters is why we termed this final result a NPW particle filter in the main text. The advantage of the NPW particle filter is that it allows for an approximate simulation of the conditional master equation (1) that includes interesting and important quantum correlations in the atomic field, yet still scales more favourably than direct integration of the conditional master equation (see comment on curse of dimensionality in section 2.1).
It is convenient to express equation (1) in terms of the number operator and phase operator (as defined by Susskind and Glowgower [102]) before transforming to the NPW representation. For most of the terms in equation (1) this is straightforward, as we simply use the definition of the local number operatorn(x) =ψ † (x)ψ(x). The most difficult term is the single-particle Hamiltonian h(x, u), as it may contain derivatives and so we cannot simply assume that it commutes with the field operators. Instead, we define the functional h (x, y, u) = h(y, u)δ(x − y), giving where h (x, y, u) does commute withψ(x) andψ † (y). We now change the field operators to number and phase operators using the identitieŝ These correspondences are used to transform the Hamiltonian and measurement operators of equation (1) tô respectively, which is correct up to a constant offset in the single-particle Hamiltonian, h(x, u). The next step is to translate the evolution of the conditional density matrix,ρ, in equation (1) to the evolution of a functional NPW distribution N [n(·), ϕ(·)]. We define the functional NPW representation by applying the single-mode NPW representation defined in [30] to every point in space: where n(·) k(·)=−n(·) denotes a sum over k(x) from −n(x) to n(x) for every point x and x denotes a product over every point x. A functional distribution is required because we are simulating a quantum field. For example, the quasi-probability that there are n 0 particles with phase ϕ 0 at a spatial point x is given by N [n 0 δ(·), ϕ 0 δ(·)]. For additional details on functional quasi-probability distributions, see [103,104,47]. Finding the evolution of N [n(·), ϕ(·)] is achieved using correspondences derived in [30]. We present the functional version of these correspondences below: where ∂ ϕ(x) is a functional partial derivative and dϕ(·) is a functional integration measure.
Applying the correspondences (A.6)-(A.11) to equation (1) using the modified operators (A.3) and (A.4) gives the following evolution for the functional NPW representation: where c.c. denotes the complex conjugate, and f [n(·), ϕ(·)]N [n(·), ϕ(·)] (A.14) is our notation for an expectation using the NPW quasi-probability distribution, where ∞ n(·)=0 denotes a sum over n(x) for every point x. The evolution presented in equation (A.12) is exactly the same as that given by conditional master equation (1), as we have not yet applied any approximations. There is also a one-to-one mapping betweenρ and N [n(·), ϕ(·)], but in practice we rarely need to reconstructρ entirely. Instead, a much more useful connection between the two is through the generation of moments. In [30] it was shown that anti-normally ordered field operators are related to NPW expectations by The solution to the equation of motion for the quasi-probability distribution [equation (A.12)] gives identical physical results to the conditional master equation (1). However, it also suffers from the same curse of dimensionality as the conditional master equation. In order to address this issue we look to classical probability. In classical probability theory, all probability distributions whose evolution is continuous in time and Markovian must obey a Fokker-Plank equation [47]. Fokker-Planck equations suffer from the same curse of dimensionality as quantum systems insofar as the size of the probability distribution grows exponentially with the number of degrees of freedom. Fortuitously, for every Fokker-Planck equation there is an equivalent set of stochastic differential equations (SDEs) that can be simulated efficiently. The moments of the SDEs are simply averaged many times to produce the same expectations as the Fokker-Planck equation. This 'trick' is how phase-space methods such as TW escape the curse of dimensionality. Specifically, when applying TW, higher-order derivatives are truncated such that the quasi-probability distribution then obeys a Fokker-Planck equation, which can be unravelled into a set of SDEs. We apply a similar method here with the NPW particle filter. However, since the system is continuously-monitored we instead aim to produce a Kushner-Stratonovich equation which can then be unravelled into a set of WSDEs, as detailed in [84].
We note that the terms in equation (A.12) generated by the measurement and, perhaps surprisingly, the nonlinear inter-atomic interactions are all in Kushner-Stratonovich form and can be unravelled exactly. However, the terms on the first three lines, which were generated by the single-particle Hamiltonian, are not. Fortunately, in typical BEC experiments, we can safely assume that the number of particles, N , is much larger than the number of occupied modes, M . This allows us to make the following three approximations: (i) By definition of the NPW representation, terms of the form dϕ (·)e 2i(ϕ(y)−ϕ (·))(n(y)+1 /2) are measures of the state of the system's overlap with the vacuum. If N M then there will only be a negligible occupation of the vacuum, and thus we can simply set these terms to zero.
(ii) The square roots on the first line of equation (A.12) can be expanded and truncated to second order in derivatives with respect to ϕ(x). This is essentially the same approximation as required by TW. The approximation is valid provided functions under the square root are sufficiently smooth, which is guaranteed when N M .
(iii) The discrete variable n(·) is taken to the continuous limit, which is an excellent approximation when N M .
These approximations were all developed and verified in [31,32]. After the application of these approximations, equation (A.12) simplifies to dN [n(·), ϕ(·)] = 1 dx dyh (x, y, u) − i∂ n(y) (n(x) + δ(0)/2) (n(y) + δ(0)/2) and dα(·) is a functional integration measure. Since equation (A.17) is in Kushner-Stratonovich form, it can be unravelled into an equivalent set of WSDEs [84]: where k is an integer indexing the K stochastic paths, ζ (k) i (t) is a 'fictitious' Stratonovich noise increment, which is generated independently for each of the K paths, and is the definition of the weighted average over the swarm. Equations (A. 19) and (A.20) are restated in the main text as equations (15) and is what we refer to as the NPW particle filter.
(A. 23) in the limit that the number of paths, K, goes to infinity [84]. If a truly infinite number of paths was required for this equality to hold, we would not have solved the curse of dimensionality as claimed. Fortuitously, it has been demonstrated that this equality approximately holds with vastly less memory than would be required to simulate the NPW representation directly. Typically, choosing K somewhere between 10 3 −10 6 works very well for simulations of monitored BECs [31,32]. We can use equation (A.15) and equation (A.23) to generate any quantum expectation value from the NPW particle filter. For example, the expectations of the field operator and the one-body correlation correlation matrix are respectively.

Appendix B. Initial conditions for simulations
In this appendix we discuss the initial conditions used to simulate equation (1) via both the Hatree-Fock method and the NPW particle filter. To begin, we note that all simulations in this paper assume that the atomic field is in the BEC phase, where all atoms occupy the same quantum state. However, the atoms are not initially in the ground state of the trap, and are instead in a far from equilibrium state with an order parameter described by equation (32). Given that the Hartree-Fock approximation is only truly valid in the zero temperature limit, initializing the field in a BEC state allowed for a fair comparison between semiclassical Hartree-Fock and NPW particle filter simulations. Furthermore, insight into the effectiveness of the feedback control is given by choosing an initial condition with much higher energy than the ground state, as this provides a qualitative picture of the rate at which the feedback removes energy from the system. In principle, this initial condition could be experimentally prepared by lowering interactions in the condensate (through a Feshbach resonance [1]), waiting for the atoms to reach the Gaussian ground state of the trap, and then performing a rapid change in the position of the trap. Of course, this initial condition does not correctly describe finite temperature and strongly-correlated Bose-condensed clouds of atoms. Nevertheless, we anticipate that our conclusions on the effectiveness of the feedback control qualitatively hold more generally, as the stability and robustness of the feedback mechanism ensures that the long-term behaviour of the system is independent of the initial condition.
Although the Hartree-Fock method and the NPW particle filter use the same initial condition [equation (32)], the numerical implementation is very different. The Hartree-Fock method simply sets the value of the order parameter φ(x, t) at t = 0. In contrast, the NPW particle filter requires a sampling of the swarm of fields, α (k) (x), and weights, w (k) , from an appropriate NPW representation of the initial state's density matrix. The remainder of this appendix outlines the details of these implementations.
For a BEC at temperature T ≈ 0 K, which is well-described by the mean-field wavefunction α 0 (x) [55], the density matrix iŝ where |α 0 (·) c means a field of coherent states that has the propertŷ The integral over θ is included because the mean field of a BEC is known only up to an absolute phase. This density matrix can equivalently be written (in first quantized notation) as a mixture of the Hartree-Fock states presented in equation (8): dx dyα 0 (x, t)α * 0 (y, t)|x i i y| , (B.3) where N = dx|α 0 (x)| 2 andα 0 (x) = α 0 (x)/ √ N . Note that the probability distribution for the total number of atoms in the condensate is a Poissonian distribution with mean N . The equivalence of representations (B.1) and (B.3) reaffirms that the Hartree-Fock approximation and coherent-state-based mean-field approximation are closely related.
Since the Hartree-Fock approximation assumes that the total number of atoms in the condensate is fixed, the representation (B.3) cannot strictly be used as the initial state in Hartree-Fock simulations. However, if N 1 then the uncertainty in the total number of atoms is small. For Hartree-Fock simulations, we are therefore justified in setting the total atom number to the mean number N and choosing the initial condition φ(x, t = 0) = α 0 (x) √ N . (B.4) For the NPW particle filter, the initial conditions must be chosen such that weighted averages give the same values as the corresponding quantum expectation values generated from equation (B.1). As with phase-space methods such as TW, the obvious distribution from which to sample from is the NPW quasi-probability representation for the initial density matrix (B.1): N α 0 [n(·), φ(·)] = n(·) k(·)=−n(·) x n 0 (x) n(x) e 2ik(φ(x)−φ 0 (x))−n 0 (x) 2π (n(x) − k(x)))!(n(x) + k(x))! , (B.5) where n 0 (x) = |α 0 (x)| 2 and ϕ 0 (x) = arg[α 0 (x)]. However, the quasi-probability distribution in equation (B.5) can have negative values, and thus cannot be efficiently sampled [105]. Instead, we sample from an approximate distribution that is valid when the number of particles N is much larger than the number of simulated modes M . The details of these approximations and a derivation of the resulting approximate distribution can be found in [31,32]; here we simply present the result: where ψ (1) (x) is the trigamma function. In essence, this probability distribution is the product of a Poissonian distribution for the number variable n(x) multiplied by a Gaussian distribution for the phase variable ϕ(x), for every point in space. When sampling this distribution we keep n(x) quantized, even though it is treated as a continuous variable during evolution. Importantly, this sampling combined with the evolution of the NPW particle filter [equations (A. 19) and (A. 20)] ensures that the total number of atoms per sample remains quantized. Sampling from equation (B.6) only approximately generates quantum moments equivalent to equation (B.5). On the occasions when sampling returns n(x) = 0, the moments generated become much less accurate. This is because the Gaussian distribution for the phase variable is not flat, which means that there will be some phase information in the sample of ϕ(x) = 0 when n(x) = 0. In other words, sampling ascribes some phase information when we occasionally sample vacuum states, which is not physically appropriate. To remedy this situation, we treat times when we sample n(x) = 0 as special cases where ϕ(x) is instead sampled from a flat distribution across phase rather than from a Gaussian distribution. As shown more completely in [32], the application of this technique tends to vastly improve the accuracy of the quantum moments generated. Even with this improvement the sampling technique still only valid N M (meaning the number of particles N is much larger than the number of occupied modes M ); the same approximate limit the evolution is valid in.
To summarize, we present the following algorithm for sampling the swarm that makes up a NPW particle filter when initial state is a BEC with mean field α 0 (x): (i) Run steps (ii) -(vi) for every position x and k ∈ (1, K).
(iv) Sample a random variable from a uniform distribution over the interval [0, 2π), and store it in ϕ (k) (x). Go to (vi).
(v) Sample a random variable from a Gaussian distribution with mean ϕ 0 (x) = arg[α 0 (x)] and variance 1 4 ψ (1) (n(x) + 1), store it in ϕ (k) (x). (vi) Store the weight w (k) = 1 and stochastic field α (k) (x) = (n + 1/2)δ(0)e iϕ(x) . (B.7) Once the swarm has been sampled it can the be integrated using equations (A. 19) and (A.20). We must perform a new and independent sample of the swarm every time we integrate a new virtual measurement record, η i (t), or, in other words, each time we perform a new 'experiment'. Finally, we note that the initial state sampling and the validity of the simulations are intimately related. As previously stated, simulations must have a large particle number N compared with the number of occupied modes M to ensure that both the initial sampling and evolution of the NPW particle filter stay accurate. For all simulations performed in this paper, the number of particles was N = 100 and the number of occupied modes was less than M ≈ N/3, which satisfies the approximation. We were unable to perform simulations with a higher number of atoms while the system was being monitored. Simulating a BEC under measurement involves a process of evolving and resampling sets of weights, which becomes more numerically demanding as the number of atoms increase. Fortuitously, the NPW particle filter is exact for the terms corresponding to the measurement, harmonic trapping potential and the nonlinear inter-atomic interactions. Consequently, the requirement that N M is entirely due to the approximate form of the kinetic energy term. Thus, the validity of the approximations on the kinetic energy term can be verified independent of the measurement. In particular, simulations were performed on a BEC without an interatomic nonlinearity or measurement using the NPW particle filter and compared to an exact solution. They were found to match closely in the N = 100 particle limit considered in this paper.