Synchronization along quantum trajectories

We employ a quantum trajectory approach to characterize synchronization and phase-locking between open quantum systems in nonequilibrium steady states. We exemplify our proposal for the paradigmatic case of two quantum Van der Pol oscillators interacting through dissipative coupling. We show the deep impact of synchronization on the statistics of phase-locking indicators and other correlation measures defined for single trajectories, spotting a link between the presence of synchronization and the emergence of large tails in the probability distribution for the entanglement along trajectories. Our results shed new light on fundamental issues regarding quantum synchronization providing new methods for its precise quantification.


I. INTRODUCTION
Synchronization is one of the most universal manifestations of emergent cooperative behavior, observed in a broad range of physical, chemical and biological systems [1,2]. It can arise spontaneously as a progressive adjustment of rhythms between oscillatory units due to their weak interaction and despite their different natural frequencies. Appealing examples with interesting applications comprise synchronization between hearth cardiac pacemaker cells [2], chaotic laser signals [3] or micromechanical oscillators [4][5][6].
In classical systems, spontaneous synchronization is usually characterized through the trajectories in phasespace [2]. In contrast, measuring synchronization in open quantum systems becomes more challenging and different avenues have been explored. For instance, temporal correlations in local observables can be quantified by using the Pearson correlation coefficient [10] or global quantum correlations can be addressed through the synchronization error [12]. Quantitative measures of phase-locking based on the expectation values of different non-local correlators [12,17,19,38] have been proposed, but they are often not indicative of the underlying processes [39]. Phase distributions computed from the Wigner quasiprobability distribution [14] or using phase states [38] have been used to gain extra insights in this context. Finally, information measures of correlations like the mutual information [40] or Renyi-entropies [41] have also been employed. In all these approaches, synchronization is computed through the expectation values of different (local or global) observables on the system density operator, as given by the solution of some suitable master equation.
In this paper we aim to go beyond the average effects of noise, and characterize synchronization along individual quantum trajectories in Hilbert space. The quantum trajectory approach describes the stochastic evolution of the pure state of the system of interest when environmental monitoring is available [42,43]. This formalism allows for a deeper notion of synchronization in the quantum regime, and enables one to explore a hidden link between the emergence of synchronization and the generation of entanglement along single stochastic realizations of the process, which cannot be inferred from the density operators.
The impressive development of experimental techniques in the last decade allowed the generation and recording of quantum trajectories in a number of platforms, including ultrahigh-Q Fabry-Perot cavities [44,45], superconducting qubits [46][47][48][49][50] and optomechanical systems [51,52]. Quantum trajectories have been used to detect phase transitions in the steady-state dynamics of dissipative quantum systems [53]. Recently, Ref. [39] provided a first clue on the potential of quantum trajectories in synchronization by using them to detect the presence of different phase-locking regimes. Here we aim arXiv:1910.03325v3 [quant-ph] 10 May 2020 to exploit at maximum the extra information that environmental measurements may offer us to give a deeper characterization of synchronization and phase-locking in the quantum regime.
We consider one of the most paradigmatic setups for the study of quantum synchronization, namely, a couple of (self-sustained) Van der Pol (VdP) oscillators weakly interacting through a dissipative coupling [32,54]. The two VdP oscillators reach limit-cycles in the long time run, where phase locking may appear depending on the trade-off between the oscillators detuning and their coupling strength. We use the statistics of phase-locked trajectories as well as other natural indicators to study synchronization, therefore extending the concept to the single trajectory case. Synchronization may strongly manifest in the shape of the distribution of phase-differences and other synchronization indicators, whose variances drop in its presence. Even if our findings are mainly illustrated using a simple system of two quantum Van der Pol oscillators, we expect our method to provide similar results in other setups.

II. MODEL AND QUANTUM TRAJECTORIES
The VdP oscillator is a nonlinear dynamical system with two different dissipative contributions: a nonlinear damping term and a pumping term powering selfoscillations. This model has been largely studied in the context of synchronization and Hopf bifurcations of classical systems [1,2]. In the quantum case, the model of two quantum VdP oscillators interacting through dissipative coupling can be described with the help of the following Lindblad master equation ( = 1) [32,54] where ρ is the density operator of the two oscillators, i a i is the system Hamiltonian with frequency detuning ∆ω ≡ ω 2 − ω 1 , and we denoted the dissipators as D[L]ρ = LρL † − 1 2 {L † L, ρ} for any Lindblad operator L. The positive rates V , γ (i) ↓ and γ (i) ↑ , stand respectively for the coupling strength between oscillators, and the rates at which nonlinear damping and pumping processes occur. The angle θ will determine the phase difference between oscillators at which synchronization occurs.
The classical equations of motion for the oscillators amplitude are recovered for the annihilation operator expectations α i = a i ρ = Tr[a i ρ] (first order moments) in the infinite photon limit γ The region of parameters (∆ω, V ) for which phase-locking emerges for two VdP oscillators in the classical limit displays the usual Arnold tongue V-shape centered around ∆ω = 0 [32]. For symmetric local damping rates γ (1) ↑,↓ = γ (2) ↑,↓ , it is simply given by V = 2|∆ω| [54].
On the contrary, the quantum limit is achieved when γ [14,32]. In this case the steady state solution π of Eq.(1), obtained by solving L(π) = 0, has been interpreted as a limit-cycle [32]. The presence of offdiagonal elements in π (but not in the local states after partial tracing) is a hint of phase correlations and therefore of the presence of synchronization between the VdP oscillators, as can be indeed checked from the qualitative behavior of the approximated Wigner function [32,54]. In the following we propose a quantum trajectory approach to gain a deeper look into this issue.
The quantum trajectory formalism describes the stochastic evolution of the pure state of the system |ψ(t) , conditioned on measurements obtained from the continuous monitoring of the environment [42,43]. It has been largely used in atomic physics and quantum optics, for which the formalism was originally developed [43]. Within this approach, we can unravel the dynamical evolution given by Eq. (1) by including the backaction of the continuous measurement process of the different environmental contributions (more details about the derivation are given in Ref. [55]). We identify five Lindblad operators in Eq. (1): ↑ a † 2 , and the collective operator L 5 = √ V (a 1 − e iθ a 2 ) (notice that here we introduced the rates inside the definition of the Lindblad operators). The evolution can then be described by the following diffusive stochastic Schrödinger equation: where H eff = H −i k L † k L k /2 is a non-Hermitian (effective Hamiltonian) operator and we introduced the generalized quadrature operators X k = L k + L † k . Here we denoted A ψ(t) ≡ ψ(t)| A |ψ(t) as the quantummechanical expectation values over trajectories at time t. The random variables dW k (t) are Wiener stochastic increments associated with the continuous measurement of the operators X k . They follow Gaussian statistics with zero average over trajectories dW k = 0 and obey dW 2 k = dt. The associated currents from continuous measurements read: where ξ k (t) ≡ dW k (t)/dt correspond to a white noise contribution [42].
It is worth pointing out here that among different ways of unraveling the master equation dynamics (1), we choose the diffusive approach with continuous measurements of X k because it best provides information about the oscillators phases. Other approaches like the ones achieved by direct observation of the quantum jumps correspond to the projection of the system state in the product of local Fock basis, thereby leading to a randomization of the oscillators phases. Still, the persistence of signatures of synchronization in quantum jumps would be interesting to explore.

III. MEASURING SYNCHRONIZATION
In order to characterize synchronization between the two VdP oscillators along a single trajectory |ψ(t) generated by Eq. (2), we introduce two different quantities which will help us to characterize phase-locking and synchronization of observables. The first one is the complexvalued correlator where we recall that the expectation values are taken using the stochastic wave function, |ψ(t) . The angle of the correlator C ψ ≡ |C ψ |e i∆φ ψ characterizes the phase difference between the two oscillators. In the classical limit, when quantum fluctuations can be neglected and the annihilation operators are replaced by the ampli- In general, the best quality of phase-locking |C ψ | → 1 is reached when the two oscillators are completely corre- , indicating that ∆φ ψ is a well defined phase. The minimum value |C ψ | = 0 is instead reached when the operators are completely uncorrelated (| a † 1 a 2 | = 0) and therefore ∆φ ψ contains no information about the oscillators phases.
The statistics of phase-locking along single trajectories calculated from Eq. (4) can be compared with the phase information retrieved from the steady state solution of the master equation π. From now on, we restrict ourselves to the limit γ ↓ /γ ↑ → ∞ where the master equation can be analytically solved and we can compute the correlator C in Eq. (4) (see Appendix A). Assuming for simplicity equal rates in both oscillators γ with the average phase-difference in the steady state ∆φ π defined through tan(θ − ∆φ π ) = ∆ω/(3γ ↑ + V ), independent of non-linear damping. Averaging the indicator C ψ (t) in Eq. (4) for any t over many trajectories, we recover C π in Eq. (5). A second, complementary, measure of synchronization considers the dynamics of local observables and the corresponding Pearson correlator [26]. Focusing on the position quadratures oscillators, this reads: where δ x i ≡ x i ψ(t) − x i ψ(t) and the bar stands for the time-average over the time-window ∆t around t, that is, x i ψ(t) ≡ t+∆t/2 t−∆t/2 ds x i ψ(s) /∆t. The Pearson indicator takes values between 1 and −1 corresponding respectively to perfect temporal synchronization and antisynchronization in the dynamics of x 1 ψ(t) and x 2 ψ(t) . For completely uncorrelated signals it becomes 0. It is worth noticing that the Pearson indicator (6) does not capture synchronization of the positions of the two VdP oscillators in the average steady state dynamics as given by the density operator π, since x i π = Tr[x i π] = 0 for i = 1, 2. One could consider higher moments [10,11], but here we will see how quantum trajectories offer deeper insight in the dynamical evolution of positions, even if these vanish on average in the steady state π.
We performed numerical simulations of the two VdP oscillators system [Eq. (2)] using quantum-trajectory Monte Carlo methods [60]. In order to investigate the steady state dynamics of the system we compute Eq. (2) for pure initial states |π n sampled from the steady state distribution, π = n π n |π n π n |, according to the probabilities π n , where π n |π m = δ nm . When averaging over measurement currents, Eq. (2) reduces to the Lindblad master equation (1), where the steady state π is recovered.
In Fig. 1 we show an example of the time evolution of the modulus and phase of the correlator C ψ (t) in Eq. (4) as well as the Pearson indicator r x1,x2 (t|∆t) over a sin-gle trajectory |ψ(t) as a function of time. We focus on the transition regime to phase locking. The corresponding average values obtained from Eq. (5) are respectively the top and bottom dashed lines. The Pearson indicator changes during the evolution and drops down whenever the phase difference departs from θ = 0, consistently with the local observables on the oscillators trajectories, x i ψ(t) (i = 1, 2). Indeed the relative phase is not locked to a fixed value, displaying instead a slow time dependence, ∆φ ψ (t), which can highly depart from its average value ∆φ π = −0.008 (bottom dashed line). Still, the modulus |C(t)| shows a significant correlation in the trajectories of the VdP oscillators during this time interval, even if the average value is moderate (upper dashed line), indicating a high accuracy of the phase-difference between the oscillators. This means that trajectories which are actually not phase-locked to θ may instead contribute with a high value to |C π |, giving the (wrong) impression that the system is synchronized, and therefore spotting the necessity for looking at synchronization indicators beyond average values.
The present approach also enables us to explore the relation between the emergence of synchronization and the entanglement shared between the two VdP oscillators during single trajectories, as first considered in Refs. [56,57]. The quantum state of the two oscillators remains pure during the whole trajectory [Eq. (2)] due to the incorporation of the environmental measured currents J k (t) in Eq. (3). Therefore the entanglement entropy is a unique measure of entanglement [58,59], namely with the reduced state of oscillator 1 during a stochastic trajectory ρ ψ (t) = Tr 2 [|ψ(t) ψ(t)|], and where we denote by Tr i the partial trace with respect to degrees of freedom of oscillator i. The entanglement entropy S ψ (t) for the VdP oscillators in the quantum regime takes values between S ψ = 0 (no entanglement) and S ψ = log 2 ∼ 0.69 (maximally entangled state). The average of S ψ (t) among trajectories defines the average entanglement associated with the environmental monitoring scheme [56,57]. It is lower bounded by (but it does not necessarily correspond to) the entanglement of formation in the (mixed) steady state π, i.e. the minimum average entanglement of any possible decomposition of π into pure states. Nevertheless, in the following we will focus more on higher order moments of the S ψ distribution rather than on its average.

IV. STATISTICS OF SYNCHRONIZATION
Evaluating the different measurements of phaselocking and synchronization of observables reported above along many trajectories, we are able to numerically reconstruct the full probability distributions of the measures |C ψ (t)|, ∆φ ψ (t), r x1,x2 (t|∆t) and S ψ (t) at any given instant of time t. We denote these probability densities by P t (|C ψ |), P t (∆φ), P t (r x1,x2 ) and P t (S ψ ) respectively. Nevertheless, since the trajectories are computed for the steady-state dynamics, these probability distributions are, up to finite-size sampling errors, independent of time. Therefore, in order to reduce statistical errors, we compute their time-averaged versions from some initial time t > 0 until a final fixed time, such that |ψ(t) has sufficient time to depart from the initial state |π n sampled from π. We refer to the time-averaged probability distributions as P (|C ψ |), P (∆φ ψ ), P (r x1,x2 ) and P (S ψ ).
We find that phase-locking in the model can be detected and characterized from the shape of the probability distributions P (|C ψ |) and P (∆φ ψ ) (Fig. 2) and P (r x1,x2 ). Moreover, comparing these distributions with the entanglement probability distribution, P (S ψ ), a persistent relation between synchronization and entanglement along trajectories is observed. Synchronized trajectories tend to share a greater amount of entanglement than unsynchronized ones. That is, trajectories show high values of entanglement more often when we approach high-quality synchronization regimes (Fig. 3). The trends observed for the probability distributions of the different synchronization measures allowing the characterization of synchronization are robust, and observed for a broad range of parameters in the model (Figs. 4 and 5).
In Fig. 2(a) we show two different instances of the phase-differences probability distribution, P (∆φ ψ ), for a fixed detuning between oscillators ∆ω = γ ↑ and two different choices of the coupling strength V = {5γ ↑ , 50γ ↑ }. In Fig 2(b) we show the classical Arnold tongue (region inside the black dashed lines) together with a color map displaying |C π | in Eq. (5). In Fig. 2(c) we plot the variance of the distributions P (∆φ ψ ) and P (|C ψ |) as a function of V for same detuning. We see that for values inside the (classical) Arnold tongue, small values of V induce a phase-differences distribution smoothly peaked at ∆φ ψ = θ = 0 with a large variance. If V is increased P (∆φ ψ ) becomes sharp around ∆φ ψ = 0 and both variances Var[∆φ ψ ] and Var[|C ψ |] approach zero. The red dotted line in Fig. 2(b) correspond to the parameters used in Fig. 2(a). As can be appreciated in both Figs. 2(a) and Figs. 2(b), for small values of the coupling V a very poor phase-locking is expected even for small detunings ∆ω → 0. This is in contrast to the classical case, which predicts phase-locking inside all the region. Figure 3 shows the statistical correlations between synchronization and entanglement during single trajectories for synchronized (red circles) and unsynchronized (blue circles) regimes. Each point represents a single trajectory for which we computed (the time-averages of) ∆φ ψ and S ψ . In the top inset, the corresponding entanglement probability distributions P (S ψ ) reconstructed from the data are shown. We see that inside the good synchronization region, for V = 50γ ↑ (red circles), phase-locked trajectories show high values of entanglement (arriving up to the maximal value S ψ log 2) more frequently, as manifested in the long tail of the (red) probability distribution P (S ψ ). Instead, when synchronization is poor, V = 5γ ↑ (blue circles), this effect tends to disappear and no correlation between phase and entanglement can be inferred from the data. In this case the tail in the (blue) entanglement probability distribution is lost. This statistical correlation for the tails of the distribution provides a new link between a purely dynamical phenomenon, namely, synchronization (and in particular phase-locking) with a strong measure of quantum correlations, entanglement, along trajectories.
The shape of the probability distributions P (∆φ ψ ), P (|C ψ |), P (r x1,x2 ) and P (S ψ ) can be better appreciated in Fig. 4 where the four full probability distributions are shown for a sample of 10 3 trajectories in the regimes of nearly perfect synchronization (blue bars, V = 100γ ↑ and ∆ω = γ ↑ ), poor synchronization inside the classical Arnold tongue (red bars, V = 5γ ↑ and ∆ω = γ ↑ ), and poor synchronization outside the classical Arnold tongue (green bars, V = 20γ ↑ and ∆ω = 20γ ↑ ). In the regime of the perfect in-phase synchronization (blue bars in Fig. 4) the distribution P (|C ψ |) is highly peaked around its average value C π 0.99, the distribution of phase-differences, P (∆φ ψ ), is peaked around θ = 0, and the distribution of the Pearson indicator, P (r x1,x2 ) is peaked around the maximum value r x1,x2 1. This is accompanied by a large tail in the probability distribution of entanglement P (S ψ ). On the contrary, in the other two cases (red and green bars in Fig. 4) all the distributions for the synchronization indicators become much more flattened, spreading along all their ranges. This is a signature of a poor synchronization, even if the average values may differ in the two cases. Also in both cases P (S ψ ) becomes sharp around 0, meaning that entanglement is not produced in almost all trajectories.
In Fig. 5(a) we plot the variance of the distribution P (∆φ ψ ) (red line) as a function of ∆ω for a fixed value of the dissipative coupling strength V = 20γ ↑ . There we can see how, despite that we are still in a regime of moderate-bad synchronization, the later improves when ∆ω → 0 as expected from the classical case, since the variance of the distribution becomes small. Comparing with the variance of the entanglement probability distribution P (S ψ ) (green line), we see that it behaves in the  Fig. 4(a).
opposite way. That is, the variance of P (S ψ ) increases whenever synchronization becomes stronger, in line with the appearance of long tails in the entanglement distribution reported above. This means that the probability to see a trajectory with a high value of entanglement becomes greater when ∆ω is close to zero, that is, when we enter the region of parameters where many trajectories show good synchronization. Finally, in Fig. 5(b) we provide an additional plot showing P (∆φ ψ ) when varying the phase-locking angle θ introduced in Eq. (1). We focus on parameters leading to good synchronization (V = 100γ ↑ and ∆ω = γ ↑ ) to show the existence of phase-locking at the input angle θ also when it takes other values different from zero. Blue bars correspond to θ = 0, red bars are for θ = π/3, and green ones stand for θ = π/2. As can be seen, phase-locking is verified at the different angles θ with the probability distributions showing analogous features than for the case θ = 0.

V. DISCUSSION AND CONCLUSIONS
We have shown that synchronization can arise in quantum trajectories, here for quantum Van der Pol oscillators with dissipative coupling, providing deeper insights about the synchronization phenomenon in the quantum regime. Departures from the classical scenario are reported in the limit of almost identical weakly coupled oscillators. The monitored system also displays a clear connection between synchronization entailed by phaselocking and entanglement in quantum trajectories, spotted by the emergence of long tails in the entanglement distribution. This phenomenon is compatible with previous results showing a link between the entanglement of formation of the steady state π and synchronization in a region reminiscent of the Arnold tongue [32]. We actually find small values for the variance of P (∆φ ψ ) accompanied by high values for the variance of P (S ψ ) for similar parameters [see e.g. Fig. 5(a)]. However, the suppression of the large tails in P (S ψ ) when decreasing V inside the classical Arnold tongue region is smooth [red bars in Fig. 4(d)], in contrast to the entanglement of formation of π, which becomes suddenly zero for V ∼ 10γ ↑ [32]. A deeper comparison with the entanglement of formation in π and its reconstruction beyond the quantum limit, may be performed by extending the optimal diffusive unraveling for entanglement detection proposed in Ref. [57] to the present situation.
It would be also interesting to explore connections and possible applications to quantum control [61][62][63], quantum information processing [64][65][66] or quantum thermodynamics along trajectories [67][68][69][70]. Other possible extensions of this work include considering reactive instead of dissipative couplings as well as other systems amenable to experimental realizations, such as optomechanical systems, atomic systems or superconducting qubits. In this context, it may be relevant to extend our results to the case of imperfect detection schemes, where the finite efficiency of the detectors or the impossibility to unravel some of the Lindblad operators, leads to a description in terms of a stochastic master equation [4,42]. In such a case, the synchronization indicators introduced here could be directly applied to the corresponding stochastic density operator, but more attention should be paid to the choice of a reliable measure of entanglement along single trajectories, since it will now require a minimization procedure.
To sum up, we believe that the approach introduced in the present paper opens new possibilities of more precise determination of synchronization in the quantum regime by looking at the statistical properties of different synchronization indicators. Importantly, this approach also helps to unveil a hidden link between synchronization along single trajectories and the generation of entanglement. In addition, our results offer an operationally welldefined way to experimentally characterize quantum synchronization in systems where environmental monitoring becomes possible.  In this appendix we analytically obtain the steadystate density matrix π of the two VdP oscillators from the master equation (1) in the limit γ ↓ /γ ↑ −→ ∞. For simplicity we also assume symmetric rates in both oscillators, that is, γ ↓,↑ = γ ↓,↑ . In this limit, the VdP oscillators are restricted to their two lowest Fock states, |0 i and |1 i since any other state is annihilated by the non-linear damping term in Eq. (1) [14]. This implies that the master equation can be mapped to a dissipative spin model of the form [32]: Here the Hamiltonian reduces to H = j=1,2 ω j σ + j σ − j , and the oscillator ladder operators a and a † are transformed in spin-flip operators σ − j = |0 1| j and σ + j = |1 0| j . Importantly, in Eq. (A1) the original non-linear damping term appearing in (2), has been replaced by a linear damping with an effective rate 2γ ↑ . This can be understood from the fact that any transition |1 → |2 in the original model promoted by the pumping term (at a rate 2γ ↑ ), will immediately decay to |2 → |0 as γ ↓ → ∞, leading to an effective transition |1 → |0 .
Following Ref. [32], the steady state solution of Eq. (A1) is obtained from L(π) = 0, whose non-zero elements read: and we have 10|π|10 = 01|π|01 , and 10|π|01 * = 01|π|10 . Here we introduced The marginal states of the two VdP oscillators can be computed by partial tracing π over the complementary oscillator, π i ≡ Tr j [π] = (1 − p) |0 0| + p |1 1|, where i = j and The marginal states π i have free-phase (no off-diagonal elements) and a population ratio between ground and excited states given by 1 − p : p. We note from Eq. (A6) that p is strictly greater than zero whenever γ ↑ is finite. Therefore the oscillators never collapse to their ground states. We also have p ≤ 1/3, the maximum being reached in the limit of uncoupled oscillators, V → 0, where the population ratio becomes 2 : 1. Increasing ∆ω and V we obtain lower values of p, and in the limit V γ ↑ we have p → 1/8 independently of the detuning ∆ω. The presence of larger drops in p outside the classical Arnold's tongue region for moderate values of V > γ ↑ can be seen as a manifestation of the classical phenomenon of oscillations collapse in the quantum regime as discussed in Ref. [72].
Using Eqs. (A2)-(A5) we can now calculate the value of the complex-value correlator C introduced in Eq. (5) for the steady state π. We obtain: where ∆φ π is the phase difference of the two coupled VdP oscillators defined through: In Fig.2(b) of Sec.IV, we plot |C π | as a function of V and ∆ω and compare to the Classical Arnold tongue. As can be seen there an important region inside the Arnold tongue corresponding to small values of the detuning ∆ω and small values of V (as compared to γ ↑ ) |C π | can be far from 1. This implies a smooth transition from nosynchronized to synchronized regimes. In the transition regime, the average phase difference between the oscillators, ∆φ π , may therefore be poorly informative due to the presence of quantum fluctuations.

Supplemental Material: Synchronization along Quantum Trajectories
In this Supplemental Material we provide more technical details on the dynamical evolution of the two Van der Pol oscillators under environmental monitoring. In particular we include a derivation of the diffusive stochastic Schröndiger equation employed in the main text.

DIFFUSIVE STOCHASTIC SCHRÖDINGER EQUATION
We show how to obtain the diffusive stochastic Schrödinger equation (2) in the main text following the derivations in Refs. [S1-S3]. Our starting point is the Lindblad master equation (1) in the main text, which we will unravel by using a generalized homodyne detection scheme. For convenience we will rewrite (1) as: for the Lindblad operators ↑ a † 2 , and the collective operator , which include the corresponding rates. We notice here the following Gauge symmetry of (2), for which a double transformation L k → L k = L k + l k and . Therefore we substitute L k and H by L k and H in Eq. (S1) and unravel it using the standard direct detection scheme. When the reservoir is assumed to be made of harmonic modes, like electromagnetic radiation, adding the displacement l k to the Lindblad operators corresponds to the implementation of Homodyne detection schemes [S1]. Here we apply the same unraveling methods in a generic situation having in mind the same physical interpretation as in the Homodyne measurement of field-quadratures [S2].
The evolution is split in an infinite sequence of intervals of infinitesimal duration dt, where the dynamics is updated according to a completely positive and trace preserving (CPTP) map ρ t+dt = E(ρ t ) = n M n (dt)ρ t M † n (dt) with Kraus operators: with X k = L k e −iϕ k + L † k e iϕ k , and l k = |l k |e iϕ k . Here the operators M k correspond to the detection of a jump of type L k in the dynamical evolution, while the operator M 0 stand for the intervals where no jumps of any type are detected. Assuming that at time t the state of the system is the pure state |ψ(t) , their probabilities read where A ψ(t) ≡ ψ(t)| A |ψ(t) is the expectation value along the trajectory at time t. It can be easily verified that P 0 (dt) + k P k (dt) = 1. As can be readily appreciated from Eqs. (S4), whenever |l k | is order 1, the probability of having any jump L k is only of order dt, while the probability of having no jumps during the interval dt is of order 1. Therefore the different type of jumps correspond to Poisson processes, almost all the time no jumps of type L k will be detected, and the evolution of the system will occur according to the operator M 0 (dt). That is: (S6) which corresponds to a smooth non-unitary evolution. On the other hand, at some (rare) instant of times, where a jump k is detected, the system state changes as: The stochastic Schrödinger equation can be constructed by introducing the number of jumps of each type k detected until time t, N k (t). Whenever the probabilities P k (dt) remain of order dt the number of jumps fulfill Poisson statistics and the associated stochastic increments dN k (t) fulfill dN k (t)dN l (t) = δ kl dN k (t), with average over trajectories dN k (t) = P k (dt). The quantities dN k (t) are stochastic variables taking values either 0 (when no jumps are detected) or 1 when a jump k is detected. The infinitesimal time-evolution of the system d |ψ(t) ≡ |ψ(t + dt) − |ψ(t) can then be written in Itô form as a sum of the different pieces of the evolution: which, by noticing that dtdN k (t) ∼ O(dt 2 ), leads to the standard form of the stochastic Schrödinger equation for jumps L k = L k + l k : Here we are interested in a continuous description, where the Poissonian statistics of the jumps L k = L k +l k become a white noise. Indeed if |l k | is arbitrarily increased, we can see from Eqs. (S4) the probability of the jumps P k (dt) may become comparable to P 0 (dt). The continuous limit is achieved when the jumps become very probable, but their effect on the system is very small. We then consider a coarse-grained evolution such that many jumps are detected in every single time interval ∆t but the change in the system is still infinitesimal (see Ref. [S2]), that is ∆t 3/2 1 and |l k | −1 . In this case the central limit theorem can be applied, and the probability distribution for the number of jumps ∆N k of type k during ∆t becomes Gaussian: where ∆W k is a Wiener increment verifying ∆W k ∆W l = δ k,l ∆t. The (unnormalized) state of the system after ∆t depends on number of jumps J detected during the interval and their precise sequence {(t J , k J ), ..., (t 1 , k 1 )}, that is: where M kj are the operators in Eq. (S2) (k j = 0) and U eff (t − s) ≡ exp −i(H + i 2 k L † k L k + |l k |X k + |l k | 2 )(t − s) describes the smooth evolution periods where no jumps are detected. Here it is worth noticing that, on the relevant time-scales ( 3/2 ), the operators in Eq. (S11) commute, so that we can approximate the unnormalized state after the different ∆N k jumps of type k as In the following we will assume for simplicity that ϕ k = 0 for all k. Expanding Eq. (S12) in orders of and keeping terms up to 3/2 we obtain: where we neglected a multiplicative term l ∆N k k exp[−|l k | 2 ∆t] irrelevant for the unnormalized state |ψ(t + ∆t) . Now taking the limit l k → ∞, so that → 0, we can replace ∆t by dt and ∆W k by dW k . The stochastic Wiener increments dW k (t) represent a white noise contribution, such that dW k dW l = δ k,l dt and average over trajectories dW k = 0. Including normalization, we obtain the final form of the diffusive stochastic Schrödinger equation: which, upon identifying H eff ≡ H − i 2 k L † k L k , matches the form reported in Eq. (2) of the main text. The output currents associated to the measurements can be obtained by removing the constant displacement from the signals in Eq. (S10), and taking the continuous limit: where ξ k (t) ≡ dW k (t)/dt, corresponding to a continuous measurement of the quantity X k . Finally, from Eq. (S14) we can calculate the corresponding stochastic master equation for the conditioned density operator (t) ≡ |ψ(t) ψ(t)|. It reads: where we introduced the measurement superoperator H[L]( ) = L + L † − X k (t) . We notice that the above stochastic master equation is in the general form reported in Refs. [S4, S5] for ideal (efficient) detectors. Taking the average over trajectories, we can easily verify that since dW k (t) = 0, Eq. (S16) reduces to the standard master equation (1) of the main text.