Quantum filter reduction for measurement-feedback control via unsupervised manifold learning

We derive simple models for the dynamics of a single atom coupled to a cavity field mode in the absorptive bistable parameter regime by projecting the time evolution of the state of the system onto a suitably chosen nonlinear low-dimensional manifold, which is found by use of local tangent space alignment. The output field from the cavity is detected with a homodyne detector allowing observation of quantum jumps of the system between states with different average numbers of photons in the cavity. We find that the models, which are significantly faster to integrate numerically than the full stochastic master equation, largely reproduce the dynamics of the system, and we demonstrate that they are sufficiently accurate to facilitate feedback control of the state of the system based on the predictions of the models alone.


Introduction
A two-level atom coupled to a driven and observed cavity mode exhibits a variety of interesting dynamics [1], including scenarios where trajectories of the atom-cavity state tend to localize transiently but jump between multiple regions of phase space on longer timescales.Here and in related work, we loosely refer to such behaviour as 'bistability', and we use the term 'stable region' or 'attractor' to refer to the local regions of phase space, where the system spends most of its time.Bistable systems have potential applications as memory units and switches, and this motivates studies of ways to understand and control their dynamics.The phase bistable regime, where the system has two stable regions with different values of the phase of the cavity field, has been investigated in several papers [2,3,4,5], and it has been demonstrated that quantum jumps between the two stable regions can be observed in the photo current from a homodyne detector monitoring the field leaking out of the cavity [4].As shown in [6], there is also an absorptive bistable regime, for which the stable regions have different values of the amplitude of the cavity field mode.Quantum jump behaviour is observed in this case as well, but it is more complicated to obtain simple approximate descriptions of the dynamics due to lack of symmetry between the two stable regions [7], and we thus consider this regime in the following.Examples of experimental investigations of bistability in cavity quantum electrodynamic systems are provided in [8,9,10,11].
The time evolution of the state of a continuously monitored quantum system is governed by a stochastic differential equation, but it is typically a very slow process to integrate this equation numerically due to the large dimensionality of the Hilbert space.This is, in particular, a problem, if we would like to control the system dynamics through feedback, since, in that case, it is necessary to track the state of the system in real time.In many cases, it turns out that the dynamics does not explore all degrees of freedom in the full Hilbert space, and this opens the way to develop simple low-dimensional models, which can, at least approximately, predict the time evolution of the state of the system.One way to obtain such models is to project the system dynamics onto an affine linear subspace of low dimension, and this technique has turned out to be very successful in the case of phase bistability [5,7], while the results for absorptive bistability are less satisfactory [7].It is, however, quite possible that improved results can be obtained by considering the more flexible case of projection onto a nonlinear manifold.In the present paper, we demonstrate that the latter approach provides simple models of absorptive bistable dynamics, which are sufficiently accurate to allow feedback control of the state of the system.
The paper is structured as follows.In section 2, we introduce the system and integrate the stochastic master equation to provide examples of absorptive bistable dynamics and quantum jumps.In section 3 we derive reduced models of the behaviour of the system by first identifying a low-dimensional manifold, which captures most of the dynamics, and then projecting the full system dynamics onto that manifold.The ability of the reduced model to reproduce the results of the full model is investigated in section 4. Finally, in section 5, we demonstrate that a feedback scheme, which builds only on predictions of a reduced model, can be used to hold the system at one of the stable regions.Section 6 concludes the paper.

Absorptive bistability and quantum jumps
Our model system is a two-level atom with ground state |g and excited state |e coupled to a cavity field mode with coupling strength g as illustrated in figure 1.The cavity mode, which decays at a rate 2κ, is driven by a coherent laser beam, and light reflected from the cavity is observed with a homodyne detector.The excited state of the atom decays at a rate γ by spontaneous emission, but since the emitted photons travel in random directions, it is so far not experimentally feasible to detect all of them with high efficiency.We thus assume no detection of spontaneously emitted photons in the following and use a density operator ρ to represent the state of the atom and the cavity field mode.
The time evolution of ρ in a frame rotating with the frequency of the drive laser is determined by the stochastic master equation (see, for instance, [12] for a derivation) with Hamiltonian Here, â is the cavity field annihilation operator, σ = |g e| is the atomic lowering operator, ∆ c is the detuning between the cavity resonance frequency and the frequency of the drive laser, ∆ a is the detuning between the atomic transition frequency and the frequency of the drive laser, and E = √ 2κβ, where |β| 2 is the average number of photons in the probe beam arriving at the cavity input mirror per unit time and β is assumed to be real.The phase φ, which is varied experimentally by varying the relative phase of the probe beam and the local oscillator, determines which quadrature of the cavity field is detected.The x-quadrature is measured for φ = 0 and the p-quadrature is measured for φ = π/2.Finally, the Wiener increment dW is a Gaussian stochastic variable with mean 0 and variance dt, which is related to dy, the observed homodyne photo current (in units of photons per time) integrated from t to t + dt divided by the square root of the average number of photons per unit time in the local oscillator beam, through dy = dW + √ 2κ Tr(âρe −iφ + ρâ † e iφ )dt.
In an experiment, the photo current is measured as a function of time, and we can eliminate dW between (1) and (3).In numerical simulations, on the other hand, we use a random number generator to obtain realizations of dW and integrate (1) directly.
Useful insight into the dynamics predicted by (1) can be obtained through various semiclassical approximations [1,13].We shall not pursue such models further here, but simply choose a set of parameters for which the dynamics of the system has been shown to be absorptive bistable [7]: ∆ c /γ ⊥ = 0, ∆ a /γ ⊥ = 0, κ/γ ⊥ = 0.1, g 0 /γ ⊥ = √ 2 and E/γ ⊥ = 0.56, where γ ⊥ = γ/2 is the transverse atomic decay rate.We assume throughout that the initial state of the system is the state with zero photons in the cavity and the atom in the ground state, and we use the second order derivative free predictor-corrector method of [14] to integrate (1).The expectation value of the number of photons in the cavity is typically below 22, and we thus truncate the basis of the Hilbert space of the cavity field at 59 photons, leading to a density matrix of dimension 120 × 120.
In figure 2 we show results for the time evolution of â + â † /2 and the expectation value of the Pauli operators σ x = σ + σ † and σ z = [σ † , σ]/2 for a given realization of the measurement noise dW .The system is seen to jump between two stable regions with different expectation values of the operators even though â + â † /2 fluctuates more in the upper region than it does in the lower region.The symmetry of (1) dictates that σ y = i σ − σ † and −i â − â † /2 are both zero for the case of homodyne detection of the x-quadrature, while they fluctuate randomly around zero for the case of homodyne detection of the p-quadrature.Note that since the time average of −i â−â † /2 is zero for both of the stable regions, our ability to distinguish the regions through a measurement of the p-quadrature relies on the fact that the time evolution of the state is different in the two regions.

Model reduction
We now turn to the problem of deriving a simplified low-dimensional differential equation, which, ideally, contains the same dynamics as the full stochastic master equation (1).To do so, we first need to identify a suitable manifold onto which we can project the dynamics.Several manifold learning strategies have already been investigated in the literature [15,16], and here we use the method of local tangent space alignment (LTSA) [17].An advantage of this method is that it optimizes the choice of low-dimensional space in local areas, which means that it is well suited to describe systems with more than one stable region.LTSA defines the manifold in terms of single  points, but in order to perform the projection, we need a differentiable function, which relates the coordinates of all points in the low-dimensional space to the coordinates of the same points in the full space.This problem is solved by fitting a function to the points obtained from LTSA.

Identification of the low-dimensional manifold
In brief, the input to the LTSA algorithm is a set of N vectors x (i) , i = 1, 2, . . ., N, sampled with noise from an unknown d-dimensional (nonlinear) manifold embedded in an m-dimensional space, where m > d, and the objective is to identify the underlying ddimensional manifold.For a linear manifold this is done by computing the d-dimensional affine subspace, which minimizes the sum of the square of the errors between the original vectors and the vectors projected onto the affine subspace.To tackle the more general case, the LTSA procedure aligns local linear structures into a global nonlinear manifold, where the local structures are the affine subspaces obtained by applying the above procedure to subsets of the N points.The ith subset is chosen as the k nearest neighbours of the ith point (including the point itself), where k is a number satisfying d ≪ k ≪ N. The output of the algorithm is the coordinates τ (i) of the points in the d-dimensional space and an approximate map from the d-dimensional space to the full m-dimensional space, from which it is possible to compute corrected coordinates x(i) of the points in the m-dimensional space.The map is, however, only valid in small regions around each τ (i) , and it is not differentiable at all points.
In our case, we start from a set of density matrices sampled from the time evolution of the state of the system.Concretely, we choose the density matrix at times t = 501 γ −1 ⊥ , t = 502 γ −1 ⊥ , . .., t = 2500 γ −1 ⊥ for the trajectory used to compute the results shown in figure 2, and we choose k = 60 as in [7].These density matrices are transformed into real column vectors x (i) , each with 14400 elements, by concatenating the real part of the columns of the upper right triangular part including the diagonal followed by concatenation of the imaginary part of the columns of the upper right triangular part excluding the diagonal, i.e., x x where n = 1, 2, . . ., N in (4) and n = 1, 2, . . ., m − 1 and m = 2, 3, . . ., N in ( 5) and ( 6).This construction method ensures that every vector x in the m-dimensional space corresponds to a Hermitian matrix.Furthermore, the LTSA algorithm ensures that the new points in the m-dimensional space are correctly normalized.There is, however, no guarantee that the constructed points correspond to positive semi-definite matrices, and here we rely on the ability of the time evolution equation to keep the state of the system within the physically acceptable region.Depending on the purpose of the reduced model, it is not necessarily optimal to minimize the projection error with respect to the dot product (x (i) ) T x (j) , and one could, for instance, consider to multiply the density matrix elements with different weight factors [7].We note in particular that (x (i) ) T x (j) is equal to Tr(ρ (i) ρ (j) ) if a factor of √ 2 is included on the right hand side of ( 5) and ( 6), but we have avoided to do so in the following, because we obtain better results without the factor √ 2. Having obtained a set of points τ (i) in the low-dimensional space and the corresponding coordinates x(i) in the full space, we next construct a map from the low-dimensional space to the full space via fitting.We need to compute one fit for each of the m = 14400 coordinates in the full space, and to make this procedure practical, we use the same fitting model for all the coordinates and choose this model to be linear in the fitting parameters, i.e., we assume a map of form x = cf (τ ), where x is a vector in the m-dimensional space, c is an m × r matrix of fitting parameters, and f is an r × 1 vector, whose elements f j are arbitrary functions of the coordinates τ = (τ 1 , τ 2 , . . ., τ d ) T in the d-dimensional space.To ensure that x is correctly normalized for all τ , we choose f 1 (τ ) = 1 and minimize i ( , where v is an m × 1 vector, whose ith entry is one if i = n(n + 1)/2 for some n ∈ {1, 2, . . ., N} and zero otherwise (i.e., v T x = Tr(ρ), where ρ is the density matrix corresponding to the vector x).The result is where c = [(z T z) −1 z T y] T , with z ij = f j (τ (i) ) and y ij = x(i) j , is the standard linear least squares result without constraints.

Projection of the dynamics onto the identified manifold
The result of the last subsection is a relation of form where c j is the matrix obtained by applying the inverse of (4-6) to the jth column of c.
To project the stochastic master equation onto the manifold defined by ( 8), we follow the derivation in [7].We would like to interpret dρ as a vector, but dρ only transforms as a vector if (dW ) where g = g(τ ) is the metric tensor with elements Combining this relation with we obtain an expression for the time evolution of τ i To simplify the notation, we define vectors v 1 and v 2 with elements and matrices M g , M H and M 1 with elements where M = −iH/ − κâ † â − κâ 2 e −2iφ − γσ + σ − .Finally, we insert A[ρ(τ )] and B[ρ(τ )] into (12) and convert back to itô form to obtain where Integrating the low-dimensional equation ( 19), we can now approximately predict the time evolution of ρ through (8).

Performance of the reduced models
Since we would like to use the reduced models to predict the state of the system in a feedback scheme, we should check the performance of the reduced models by generating a realistic photo current using the full stochastic master equation and then use that photo current to integrate the reduced models.Examples of this procedure, using polynomials as fitting models, are provided in figure 3. We have plotted â+â † /2, because this is the quantity we need to predict in the feedback scheme proposed in the next section.For the case of homodyne detection of the p-quadrature we have used a weighted linear least squares method to compute c in (7) to reduce the effect of outliers.The precise initial state of the reduced models is not important, because the observed value of the photo current quickly drags the reduced model to the correct state, and we have thus chosen τ = (0, 0, . . ., 0) T (the average of τ (i) ) for simplicity.In case of instability, we have reset the reduced model to τ = (0, 0, . . ., 0) T whenever the program returns a non-determined value of τ .
The figure shows that simple low-dimensional models are able to provide accurate predictions for the value of â + â † /2 for the case of homodyne detection of the xquadrature.For homodyne detection of the p-quadrature, the two-dimensional model obtained by using a second order polynomial as fitting model is observed to largely reproduce the time evolution of â + â † /2, but the details differ.In particular, the reduced model does not reproduce the highest values of â + â † /2.To improve the agreement between the full and the reduced model, one could try to increase the order of the fitting polynomial or use a higher-dimensional model.For a fourth order polynomial in two dimensions, the reduced model is able to reach the highest values of â + â † /2 and the value of â + â † /2 in the lower state is also more accurate, but the model predicts false jumps to the lower state, which is undesirable in a feedback scheme.False jumps are not observed in the figure for the four-dimensional model with a second order 2 ) T , for (b) and (e) it is a fourth order polynomial in two dimensions, and for (c) and (f) it is a second order polynomial in four dimensions.The left figures are for homodyne detection of the x-quadrature, and the right figures are for homodyne detection of the p-quadrature.polynomial as fitting model, and we thus use this latter model to predict the state of the system in the next section.
The reason why we obtain very accurate results for homodyne detection of the xquadrature is that we use the actual photo current dy obtained from the full stochastic master equation to drive the reduced model.This means that the back action of the measurement on the system is large whenever the predicted value of the measured quadrature differs significantly from the value obtained from the full stochastic master equation.It is, in fact, a much harder test of the performance of the reduced models to check whether they are able to predict the value of quantities that are not related in a simple way to the observed quadrature such as the value of â + â † /2 for homodyne detection of the p-quadrature.The above results thus indicate that the low-dimensional models actually capture most of the full dynamics of the system.As expected, we also find that the reduced models for homodyne detection of the p-quadrature provide more accurate results for −i â − â † /2 than for â + â † /2.

Stabilization of one attractor through feedback control
A natural feedback scheme to hold the system within one of the stable regions is to increase or decrease the intensity of the drive laser depending on whether the value of â + â † /2 is below or above some suitably chosen target value x 0 .This is achieved by adding a proportional feedback term to the stochastic master equation (1), where s p is a parameter determining the strength of the feedback.Integration of the resulting equation confirms that the feedback term has the desired effect, and it is possible to decrease the standard deviation of e(t) to a value, which is very small compared to typical fluctuations of â + â † /2 without feedback.The question is now whether this is still the case if we use a reduced model (with the feedback term included) to predict the value of e(t).
Example trajectories are shown in figure 4. Comparing these trajectories to those in figure 3, it is apparent that the feedback term affects the time evolution, and that the system stays close to the upper or the lower stable region for the considered values of x 0 .For the upper stable region, â + â † /2 is seen to oscillate with a relatively large amplitude, which reflects the fact that the upper stable region is relatively broad as observed in figure 2. For the lower stable region, â + â † /2 is roughly constant except for sudden spikes.By choosing a value of x 0 , which is slightly below the average value of â + â † /2 predicted by the reduced model, we can ensure that the feedback term almost always acts to decrease â + â † /2.This keeps the full model away from the transition region, and as seen in the figure the spikes tend to point towards negative values of â + â † /2 rather than towards larger positive values of â + â † /2 as observed for higher values of x 0 .
For the case of homodyne detection of the p-quadrature, we note that the predictions of the reduced model for x 0 = 3.5 fluctuate less than the results obtained from the full stochastic master equation, while the fluctuations of the full and the reduced model are approximately the same for the case of homodyne detection of the x-quadrature.This is because we use the reduced model to evaluate the error e(t), which means that the feedback term always acts to reduce e(t) for the low-dimensional model.For the full model, on the other hand, e(t) may have the opposite sign, in which case the full model is pushed away from x 0 .The resulting change in the photo current  3, but with the feedback term in (22) included.The error e(t) is computed from the value of â + â † /2 predicted by the four-dimensional reduced model with a second order polynomial as fitting model (green), and as in the last section we have used the photo current obtained from the full stochastic master equation to integrate the lowdimensional model.In (a), x 0 = 3.8 and s p = 0.75 for the upper curves and x 0 = 0.4 and s p = 1 for the lower curves.In (b), x 0 = 3.5 and s p = 0.75 for the upper curves and x 0 = 0.5 and s p = 1 for the lower curves.
drives the reduced model in the same direction as the full model, but this mechanism is more efficient in the case of homodyne detection of the x-quadrature than for homodyne detection of the p-quadrature.
Discrepancies between the full model and the reduced model lead to a feedback of noise into the system, and even though a large value of s p may reduce the variance of e(t) for the low-dimensional model, we observe that the predictions for â + â † /2 obtained from the full model fluctuate over a range that is broader than the distance between the upper and the lower stable region if s p is chosen too large.The power spectrum of the time evolution of â + â † /2 for a trajectory computed from the full stochastic master equation with homodyne detection of the p-quadrature and the power spectrum of the difference between the predictions of the reduced model and the full model in figure 5(a) show that the relative error of the reduced model in predicting â + â † /2 is smaller at low frequencies.This appears because the reduced model is able to predict quantum jumps of the system but does not capture the details of the dynamics in the stable regions, and it suggests that it might be an advantage to mainly feed back the low frequency behaviour, which can be achieved by adding an integral term to the controller dρ fb = s p e(t) where s i and ζ are constants.Rough estimates of reasonable choices of s p , s i and ζ can be obtained as follows.If the measurement noise is turned off by setting dW = 0 in the full stochastic master equation, the state of the system decays to steady state, which is an incoherent mixture of the states corresponding to the upper and the lower stable regions.(Note that dW in equation ( 19) may be nonzero, since we still require that the the photo current is the same for the reduced and the full model.)The behaviour of the system when a small input drive field term dρ u = u(t)[â − â † , ρ]γ ⊥ dt is added can then be investigated by linearizing (1) where δτ is the deviation of τ from the steady state value, δρ is a vector of the deviations of the density matrix elements from their steady state values with one diagonal element omitted, C is a 14403 × 14403 matrix and Q is a 14403 × 1 vector.Choosing x 0 to be the steady state value of â + â † /2 obtained from the reduced model, we can write the error as where v e is a 14403×1 vector for which all but the first four elements are zero.Integrating the linear equation (24), we have an expression for the time evolution of e(t), which we insert into (23) to obtain the feedback term for a given input u(t).In a closed loop setting, the feedback term is used as input, and the system behaviour is thus characterized by the loop transfer function which is the Laplace transform of the coefficient in square brackets in (23) divided by the Laplace transform of u(t).For s = iω, the norm of K is the loop gain at angular frequency ω, and according to the power spectra in figure 5(a) this should be close to zero for angular frequencies above approximately 2π × 0.02 γ ⊥ = 0.126 γ ⊥ .The norm of v T e (C − iω) −1 Q is determined completely by the system and is plotted in figure 5(b) for the case of homodyne detection of the p-quadrature.Since this factor is substantially different from zero for a range of angular frequencies above 2π × 0.02 γ ⊥ , we choose s p = 0. To ensure that |s i (ζ + iω) −1 | has a large negative derivative for ω ≈ 2π × 0.02 γ ⊥ , we set the angular cross-over frequency to ζ = 2π × 0.02 γ ⊥ .Finally, we choose The resulting norm and phase of the loop transfer function are also plotted in figure 5(b).The phase is seen to be well above −π at the cross-over frequency, and we thus expect the feedback to be stable.One should, however, not read too much into (26) as the above derivation is very crude.A similar analysis for the case of homodyne detection of the x-quadrature leads to the parameters s p = 0, s i = 0.51 γ ⊥ and ζ = 2π × 0.1 γ ⊥ and suggests that the feedback may be unstable for s i 0.7 γ ⊥ .
Keeping s p and ζ fixed and searching in the neighbourhood of the above values of s i , we find that s i ≈ 0.51 γ ⊥ is close to optimal for the case of homodyne detection of the x-quadrature.A too small value of s i (for instance, below 0.25 γ ⊥ ) leads to increased fluctuations in the predictions of the full model, because the feedback is insufficient to keep the system within the stable region, and a too large value of s i (for instance, above 1 γ ⊥ ) tends to course instability.For the case of homodyne detection of the pquadrature, we obtain improved results by decreasing s i to 0.15 γ ⊥ .Trajectories for these parameters are shown in figure 6.For the upper stable region the standard deviation of the difference between the reduced and the full model relative to the standard deviation of the trajectory obtained from the full model for the time interval from t = 100 γ −1 ⊥ to t = 750 γ −1 ⊥ is reduced relative to the value obtained for the trajectories in figure 4. For homodyne detection of the x-quadrature it decreases from 0.71 to 0.46, and for homodyne detection of the p-quadrature it decreases from 1.01 to 0.82.This confirms that the integral term feeds back less noise.On the other hand, the integral term is less efficient in keeping â + â † /2 close to the desired value, and the standard deviation of â + â † /2 obtained from the full model is observed to increase.The conclusion is the same for the lower stable region for homodyne detection of the x-quadrature, but for homodyne detection of the p-quadrature we observe a decrease in both the standard deviation of the difference between the reduced and the full model and in the standard deviation of â + â † /2 obtained from the full model.

Conclusion
In conclusion, we have shown that it is possible to derive simple low-dimensional models for the dynamics of a single atom interacting with a cavity field mode in the absorptive bistable regime, and we have demonstrated that the models can be used to construct a feedback scheme, which is able to hold the system within one of the stable regions.This is an important result, because it is unrealistic to integrate the full stochastic master equation in real time.
The suggested feedback scheme relies on predictions of the expectation value of the x-quadrature of the cavity field, and we have considered both the case of homodyne detection of the x-quadrature of the output field from the cavity and homodyne detection of the p-quadrature.The former case is relatively easy to handle, because the estimated quantity is directly related to the observed quantity.For homodyne detection of the p-quadrature, on the other hand, the expectation value of the x-quadrature has to be inferred from the precise time evolution of the p-quadrature, and the reduced model has to do significantly more work.It is thus promising for the method that we also obtain reasonable results in this case.
There are many degrees of freedom in the modelling procedure and, in general, it may require some trial and error to find reduced models that are able to reproduce the system dynamics with sufficient accuracy.In a feedback scheme, one should concentrate on optimizing the ability of the reduced model to predict the quantity that determines the feedback, since errors in this quantity lead to a feedback of noise into the system, which limits the performance of the feedback.
A further line of research could be to find systematic methods to optimize the models.One could, for instance, consider different ways to construct the vectors used as input to the local tangent space alignment procedure, which corresponds to different criteria for the optimal choice of low-dimensional manifold.It is also possible that improved models could be obtained by choosing other kinds of fitting models than polynomials of low order.The selection of the density operators used to compute the low-dimensional manifold could be adjusted according to the purpose of the model.In the case of application of feedback, one could, for instance, include more points at one attractor than the other in order to obtain a better description of the dynamics in the neighbourhood of the attractor we intend to stabilize.

Figure 1 .
Figure1.Two-level atom in a cavity probed with a coherent laser beam.BS is a beam splitter of low reflectivity, and LO is the local oscillator.

Figure 3 .
Figure 3.Time evolution of â + â † /2 obtained from the full stochastic master equation (blue) and from various reduced models (green) assuming the same photo current.For (a) and (d) the fitting model is a second order polynomial in two dimensions, i.e., f= (1, τ 1 , τ 2 , τ 2 1 , τ 1 τ 2 , τ22 ) T , for (b) and (e) it is a fourth order polynomial in two dimensions, and for (c) and (f) it is a second order polynomial in four dimensions.The left figures are for homodyne detection of the x-quadrature, and the right figures are for homodyne detection of the p-quadrature.

Figure 4 .
Figure 4. Time evolution of â + â † /2 obtained from the full stochastic master equation (blue) for homodyne detection of the x-quadrature (a) and homodyne detection of the p-quadrature (b) for the same noise realization as in figure3, but with the feedback term in (22) included.The error e(t) is computed from the value of â + â † /2 predicted by the four-dimensional reduced model with a second order polynomial as fitting model (green), and as in the last section we have used the photo current obtained from the full stochastic master equation to integrate the lowdimensional model.In (a), x 0 = 3.8 and s p = 0.75 for the upper curves and x 0 = 0.4 and s p = 1 for the lower curves.In (b), x 0 = 3.5 and s p = 0.75 for the upper curves and x 0 = 0.5 and s p = 1 for the lower curves.

Figure 5 .
Figure 5. (a) Power spectrum of a trajectory corresponding to the one in figure3(f), but integrated to tγ ⊥ = 4000.The blue curve is the power spectrum of â + â † /2 obtained from the full stochastic master equation, the green curve is the power spectrum of the difference between the reduced and the full model, and ω is the angular frequency.In both cases we have subtracted the mean value before computing the power spectrum.(b) Plots of |v T e (C − iω)Q| (dashed) and of the norm (solid) and the phase (dotted) of the loop transfer function K(iω).

Figure 6 .
Figure 6.Same as in figure4, but for the feedback term in (23).The parameters are s p = 0, s i = 0.51 γ ⊥ and ζ = 2π × 0.1 γ ⊥ for the case of homodyne detection of the x-quadrature and s p = 0, s i = 0.15 γ ⊥ and ζ = 2π × 0.02 γ ⊥ for the case of homodyne detection of the p-quadrature.