Input-output Gaussian channels: theory and application

Setting off from the classic input-output formalism, we develop a theoretical framework to characterise the Gaussian quantum channels relating the initial correlations of an open bosonic system to those of properly identified output modes. We then proceed to apply our formalism to the case of quantum harmonic oscillators, such as the motional degrees of freedom of trapped ions or nanomechanical oscillators, interacting with travelling electromagnetic modes through cavity fields and subject to external white noise. Thus, we determine the degree of squeezing that can be transferred from an intra-cavity oscillator to light, and also show that the intra-cavity squeezing can be transformed into distributed optical entanglement if one can access both output fields of a two-sided cavity.


Introduction
Interfacing static and flying quantum degrees of freedom is a key step towards the realisation of operational quantum technologies. In fact, such a requirement figures as one of the additional 'networkability' Di Vincenzo criteria, central to quantum communication, distributed computation and any sort of quantum networking [1]. In current experimental set-ups, a promising approach to this problem consists in coupling a quantum degree of freedom with cavity light, and then to adopt the light leaking out of the cavity as the flying quantum degree of freedom [2,3,4,5,6]. The latter may then be used to entangle distant trapped degrees of freedom by mixing on a beam-splitter [7] or achieve entanglement in a local site by detection of the emerging light [8]. From the theoretical standpoint, this paradigm is well described by the seminal input-output formalism developed in the eighties by Yurke, Collett and Gardiner [9,10,11].
In this paper, we consider general confined bosonic degrees of freedom (cavity light fields, trapped atoms, ions, optomechanical systems, or combinations thereof), adopt the input-output formalism to include interaction with a bath and a set of (possibly) accessible output modes, and develop a very general theoretical framework to study the correlations of the output fields. Central to our analysis is the identification of suitable output modes, which are relevant to experimental detections in practical cases and will also lead us to the analytic determination of Gaussian channels [12,13,14] relating the initial correlations of the trapped system to the output correlations of the modes. If the initial state of the system is Gaussian, which is often the case when one considers bilinear Hamiltonians, the output fields will also be in a Gaussian state and the complete characterisation of their state will hence be provided in such a case. Let us emphasise that the standard input-output approach would imply the numerical solution of rather convoluted integrations. Our treatment, instead, reduces the reconstruction of the output fields to an algebraic problem (involving only the computation of matrix inverses, matrix exponentials and the solution of Sylvester equations), which we completely characterise. Furthermore, such a treatment is wholly independent from the details of the quadratic intra-cavity dynamics and extends, up to a quadratic scaling in the computational resources, to any number of modes, even into regimes where brute-force integration would be impractical.
We then proceed to apply our theory to the relevant case of mechanical oscillators (such as trapped ions, particles, or vibrating mirrors in opto-mechanical set-ups) coupled to cavity light, and show that, given some initial degree of squeezing in the mechanical element, one can obtain squeezed travelling output light, or even distributed entanglement if one can access both output modes of a double-sided cavity. We also show how our formalism can provide a rigorous treatment of the continuous monitoring of mechanical motion, when performed via realistic (non-instantaneous) detectors and outside stationary regimes.

General Formalism
We consider a localised system comprised of n bosonic modes, whose canonical operators are grouped together in the column vectorv(t) = (a 1 , a † 1 , . . . , a n , a † n ) and subject to the most general affine quantum Langevin equation: where A is the 'drift matrix', determined by the system's own Hamiltonian and linear coupling to the environment (see below), and Here, κ j quantifies the loss rate of the j-th mode of the system, which interacts with the input operatorsâ in,j and a † in,j , whilev in ∈ C 2n is just a constant vector, satisfyinḡ v in,2j−1 =v * in,2j , accounting for the possible linear driving of the system. Given a quadratic Hamiltonian with H = H ,V 2j−1 = V * 2j , A andv in can be expressed as: where Σ is the commutation matrix Σ jk = [v j ,v k ]. Throughout the paper, we shall only consider white Gaussian input noise, characterised by v in =v in and by the following correlation functions: where σ in is related to a physical covariance matrix, as it obeys The output field is determined by the boundary condition: where we deviate slightly from the standard definitions (see, e.g., [11]) to keep the formulae that will follow more compact. One can then insert the formal solution of Eq. (1) into Eq.
Note that, in general, not all the components of the output fieldv out will correspond to experimentally accessible modes. For example, the light leaking out of an optical cavity might be collected and further processed, while the 'quantum information' that a mechanical mode dissipates into its phononic environment is lost for all practical purposes. In what follows, we shall assume that the accessible output fields do not suffer further losses before their experimental manipulation. However, as shown in Appendix D, the treatment of the lossy case simply amounts to combining the Gaussian channels described below with appropriate beam-splitter transformations.

Exponential Pulses
In many applications, one is interested in output modes characterized by an exponetial temporal profile (see Fig. 1 for an example). In general, this amounts to considering output modes of the form: where Λ = diag(λ 1 , λ 2 , ..., λ 2n ) = diag(µ 1 , µ * 1 , ..., µ n , µ * n ) is a 2n × 2n diagonal matrix with Re(Λ) < 0, such that the mode profiles decay exponentially in time, with rates |Reλ j |, and no mixing between the output modes is allowed. N = |2Re(Λ)| 1/2 is a (diagonal) normalisation matrix which guarantees that thef j form a set of bosonic modes: [f j ,f k ] = Σ jk . The matrix integrations involved in deriving the explicit solution for the vector of modesf are best dealt with by considering each componentf j separately: inserting Eq. (10) into Eq. (11) yields (see Appendix A.1) where 1 is the identity matrix in dimension 2n. In this derivation, we have assumed that the system is stable, in the sense that lim t→∞ exp(At) = 0. The vectorf can be more compactly expressed by defining the elements of the matrix X and the vectorû as per Notice that the statistics of the vectorû are known, in that they can be reconstructed from Eq. (7). In particular, one has: where X = X + N K −1 . Then, if σ(0) is the symmetrised covariance matrix (CM) of the initial state of the system (σ , obtaining: Eq. (16) defines a Gaussian quantum channel relating the system's CM at time t = 0 to the output CM σ out [12,14]. Note that it was convenient for us to write the correlation matrices in terms of the field operators a j and a † j . In order to switch to a description in terms of quadratures, and retrieve the notation more commonly used in the recent literature [14], one has simply to transform all the equations by similarity with the unitary

Slow detectors and time-dependent spectral measurements.
Let us now extend our treatment to the case of detectors with a finite response time, and to spectral analysers with finite resolution. In both cases, assuming a Lorentzian response in the frequency domain, the relevant output modes are of the form [15]: where the diagonal matrix Λ has the same form as above, and N t = |2ReΛ| 1/2 [1 − e 2Re(Λ)t ] −1/2 guarantees that the modesĝ j (t) are bosonic. As anticipated, the modes of Eq. (17) have two main interpretations. On one hand, they may describe continuous detection of the output fields with 'slow' detectors. In this case, |2Re(λ j )| quantifies the bandwidth (inverse of the response time) of the j−th detector. On the other hand, they may be used to model spectral measurements with finite resolution, so that |2Re(λ j )| quantifies the spectral resolution of the j−th frequency filter. The exact treatment of this case, while more convoluted, runs along the same lines as the previous one, and the Gaussian channel can still be obtained in closed form as (see Appendix A.2) The quantity |2Re(λ j )| may represent either the bandwidth of a realistic detector, or the frequency resolution of a spectral analyzer. and Notice that even though the matrix I cannot be given in a simple explicit form, it is the solution of the time-dependent Sylvester equation AI + IA † = e At σ in e A † t − σ in , and can be solved by standard linear algebra methods [16]. Clearly, the matrices Z and T defining the Gaussian channel are now timedependent, as one should expect since the monitored output modes themselves depend on time in this approach.

Stationary spectrum.
Let us briefly mention that, if one assumes the stability condition lim t→∞ exp(At) = 0, as well as Re(Λ) < 0, the output field admits a stationary limit reflecting the status of the asymptotic steady state inside the cavity, with CM: where all quantities are intended in the limit t → ∞.

Applications.
Let us now apply our general formalism to instances of practical interest. We will deal with mechanical oscillators interacting with cavity fields and, typically, concern ourselves with regimes where the cavity field can be adiabatically eliminated, yielding a direct coupling of the mechanical degree of freedom to the input and output fields of the cavity [2]. Exploiting this mechanism has become a standard practice to achieve cavity cooling of a mechanical oscillator, and to perform continuous detections of fluctuations in the oscillator's position [17,19,20]. Here, we consider a further application which is readily available in the same setting. We are going to consider the possibility of high fidelity transfer of the quantum state of the motional degree of freedom, typically difficult to access, to a travelling pulse of light emitted by the cavity, which can be processed via standard optical tools (this process is sometimes referred to as 'phononphoton conversion'). In particular, we shall focus on the transfer of squeezing between matter and light and on the creation of distributed Gaussian entanglement between two output modes, given initial mechanical squeezing. In addition, we show how our formalism provides a rigorous and quantitative description of the continuous detection of the mechanical motion, performed via a realistic detector and in a manifestly nonstationary regime. We investigate in detail the case in which the oscillator is provided by the motion of a trapped ion, keeping in mind that very similar results hold for optomechanical systems [in which case, Eq. (29) below can be taken as the starting point]. Coupled ion-cavity systems have been realized experimentally in [36], and cavity cooling of a trapped ion was proposed in [18].
Let us consider an ion with motional trapping frequency ν along a chosen axis, corresponding to the mode a 1 . ‡ The ion is trapped inside an optical cavity that sustains a mode a 2 of frequency ω c , in turn coupled via a Jaynes-Cummings interaction of strength g 0 to an internal transition of the ion, described by a mode a 3 of frequency ω 0 . In order to couple the ion motion to the cavity field, we drive the ion internal transition with a detuned classical laser with Rabi frequency Ω. In a frame rotating with the laser frequency ω L , operating in the Lamb-Dicke regime, and considering only quadratic terms in the bosonic operators §, we have the Hamiltonian (see Appendix B) where δ = (ω c − ω L ), ∆ = (ω 0 − ω L ) and η is the Lamb-Dicke parameter. As we shall see shortly, we are only interested in inducing virtual transitions of the ion, which ‡ We shall assume that the motion along the remaining axes is either negligible, or it can be factored out.
§ We recall that linear terms in the Hamiltonian can be accounter for by adding a constant to the vectorv in .
will remain in its internal ground state with good approximation. Therefore, we only introduce a small error if we treat the mode a 3 as bosonic, which allows us to include the effects of spontaneous emission in the two-level system, all while remaining within the scope of our formalism and yielding an input-output Gaussian channel. In all our studies, we have simply verified a posteriori that the population of excited states above the first stays negligible at all times during the dynamics. As outlined in our general theory, we include loss rates in the problem via the matrix K [see Eq. (3), where now n = 3], where κ 1 is the motional heating rate per phonon, κ 2 the cavity loss rate and κ 3 the spontaneous emission rate of the ion's internal transition. The thermal excitations of each environment are instead included in the matrix σ in [see Eq. (7)]. In particular, only the phononic environment is thermally excited in our problem, and we shall indicate its thermal occupancy byn th . Note also that in this system only the output mode a out,2 is accessible experimentally. A direct coupling between ion motion and cavity field can be achieved if the internal levels of the ion are eliminated by time-averaging the Hamiltonian (28) in the large detuning regime ∆ (ν, δ, g 0 , Ω) [21,22]. If the ion is initially in its internal ground state, one obtains the Hamiltonian where δ = δ − g 2 0 /∆ and J = ηg 0 Ω/∆. As anticipated, a Hamiltonian of this form is also readily available in optomechanical set-ups, where it provides a reliable description of the dynamics under strong driving of the cavity field [17].
The direct coupling between the oscillator and the cavity input-output fields is finally obtained by choosing δ = ν, and considering the limit ν κ 2 J n th κ 1 . Then, neglecting the counter-rotating terms in Eq. (29), along with the adiabatic elimination of a 2 (see Appendix C), yields the following input-output equations for the mechanical mode: where κ = 4J 2 /κ 2 and we have also neglected the effects of the phononic environment (n th κ 1 κ has to be checked a posteriori). As the photonic environment corresponding to a in,2 is effectively at zero temperature, Eq. (30) describes red-sideband cooling of the mechanical mode, with cooling rate κ. Moreover, the boundary condition (31) provides us with the desired link between the cavity output field (the only detectable mode) and the motional degree of freedom.
While a useful guidance, the two relationships (30) and (31) are only approximations, whose reliability does not extend to the full dynamics we will be interested to study. However, the general analysis of Sections 3 and 4 allows us to identify and treat exactly output modes that are suitable to study the applications we are interested in, without the need to resort to any approximation. We shall hence apply our general theory to those output modes, construct the input-output channels defined in Eqs. (16) and (18), and obtain exact quantitative results which are reliable even when Motional frequency of the ion: ν = 2π · 5MHz Jaynes-Cummings coupling strength: g 0 = 2π · 0.62MHz Laser Rabi frequency: Ω = 2π · 1MHz Laser-internal transition detuning: ∆ = 2π · 10MHz Lamb-Dicke parameter: η = 0.08 Cavity decay rate: κ 2 = 2π · 53kHz Spontaneous emission rate: κ 3 = 2π · 360kHz Phononic heating raten th κ 1 = 2π · 24Hz Effective cooling rate κ 2π · 19kHz Detector's bandwidth: Γ = 2π · 50MHz Table 1. Adopted values of dynamical parameters. These are obtained by considering the 230nm transition of an Indium ion [27], trapped at the center of a Fabry-Perot cavity of length L 1cm, waist w 6µm and average photon lifetime (κ 2 ) −1 = 3µs. We have considered a stabilized trap with average heating time (n th κ 1 ) −1 = 6.6ms at room temperature (T = 300K). The effective cooling rate shown above corresponds to a temporal length (κ) −1 83µs for the modef . The 50MHz bandwidth of the detector corresponds to an average response time Γ −1 3.2ns. For further details, see Appendix E.
the approximations leading to Eqs. (30) and (31) are only partially justified. In what follows, we shall consider an ion-cavity system characterised by typical experimental parameters, which we summarise in Table 1.

Transfer of squeezing between matter and light
Let us start by considering the possibility of transferring the initial state of the mechanical oscillator to a travelling pulse of light. This follows from standard inputoutput theory being applied to Eqs. (30) and (31), yielding [26] which means that the cavity emits a pulse of effective duration τ ∼ 1/κ, whose quantum state is approximately the same as the oscillator initial state. Notice that, as we shall verify in the following, Eq. (32) is only an approximation that will serve us as a guideline in the choice of the output modes, which we will however treat exactly through our general framework. In our language, the first equality in Eq. (32) corresponds to modes of the form (11), with µ 2 = iν − κ/2. Note that the parameters µ 1 , µ 3 can be set to any value, as we are only interested in the output field a out,2 . In particular, we shall focus here on the case in which the mode a 1 is initially in a pure squeezed state, and investigate quantitatively the amount of squeezing transferred to the travelling mode f . The reliable realisation of such squeezing transfer constitutes a basic requirement for   Table 1. continuous variable quantum communication and information processing. Moreover, in the medium term this set-up might be employed to generate squeezed light, with the crucial advantage of exploiting the stronger nonlinear terms allowed by material degrees of freedom with respect to light, where nonlinear interactions giving rise to squeezing are always comparatively weak. More specifically, the squeezing of trapped ions could be previously obtained either by manipulating the trapping potentials [23,24,29] or through the internal degrees of freedom [30,31,32], while in optomechanical systems, an high degree of squeezing could be obtained via indirect position measurements [33].
The squeezing of a Gaussian state with CM σ has been quantified, in dB, as max[0, −10 log 10 (σ ↑ 1 )], where σ ↑ 1 is the smallest eigenvalue of σ (note that in our notation the vacuum state has eigenvalue 1). In the following example, we have considered a trapped Indium ion, characterised by the experimental parameters shown in Table 1 (see also Appendix E). Our quantitative findings are summarised in Fig. 3(a) showing that, even though the state transfer is not accurate for the considered experimental parameters (i.e., the curve differs from a straight line of unit slope), a substantial amount of squeezing can be transferred from the trapped ion to the output travelling mode. However, the output squeezing seems to hit a plateau at about 9dB after an initial linear rise, indicating a complete breakdown of the approximated relation f −ia 1 (0) for high input squeezing. This is also reflected in the purity Tr( 2 ) of the output state which, being our state Gaussian, can be determined as 1/(detσ out ) 1/2 [34]. Fig. 3(b) shows clearly that the purity of the output state is a decreasing function of the initial intra-cavity squeezing (in dB), proving that the conversion is not perfect since part of such squeezing is not coherently transferred but contributes instead to phase-insensitive noise in the output fields.

Entanglement generator.
By considering, once again, output modes defined as in Eq. (11), we can also study the leakage of both cavity mirrors in a double-sided cavity, for an initial squeezed state of the mechanical oscillator, and consider the state of the two output modes leaking out of each mirror. Entanglement between those modes can be expected in this scenario, since the adiabatic elimination of the cavity field yields an effective coupling of the two output modes to the same motional degree of freedom. Indeed, it can be shown that the formal treatment of this case is identical to that of a single-sided cavity, up to mixing its output field with the vacuum at a beam-splitter (see Appendix D). For the Gaussian squeezed state of Section 5.1, this operation results in entanglement between the two output modes, which can be quantified in terms of the logarithmic negativity N [35].
For a two-mode Gaussian state with CM σ, the logarithmic negativity can be evaluated as max[0, − log 2 (ν − )], whereν − is the smallest eigenvalue of the matrix |ΣT σT |, with T = 1 2 ⊕ σ x (σ x being the x-Pauli matrix). Fig. 3(c) shows that a substantial amount of entanglement can be generated, and hence distributed, for realistic amounts of internal squeezing. Our findings clearly point at the potential held by the motional degrees of freedom of massive particles as competitors of nonlinear crystals for the generation of optical entanglement.

Continuous detection of the mechanical motion.
Eq. (31) suggests that, by measuring the cavity output light with a fast detector, one should be able to monitor the mechanical motion in real time. In particular we shall focus here on the continuous detection of the mechanical population, keeping in mind that a similar analysis can be performed if one is interested in measuring the mechanical quadratures (see Appendix C). We stress that here we are concerned with quantities averaged over many experimental runs, hence our treatment does not deal with the partial collapses of the quantum state that typically occur in a single realization of a continuous measurement [26].
To formulate the problem in terms of an input-output Gaussian channel, we consider modes of the form (17), with µ 2 = iν − Γ/2, Γ being the bandwidth of the detector. Again, the parameters µ 1 , µ 3 can be set to any value, as we are only interested in the detection of the mode associated to a out,2 , which we shall denote as g(t). From Eqs. (17) and (31), assuming that Γ is large compared to the other frequencies in the problem, we get (see Appendix C) Therefore, by rescaling the average intensity registered by the detector, one should be able to infer the average number of phonons in the oscillator at time t. Using our general formalism, we can now investigate quantitatively the validity of this prediction. In Fig. 3(d), we compare inferred and actual values of a † 1 a 1 , as a function of time, for an Indium ion with an initial motional squeezing of 20dB. It can be seen that, after an initial transient time, the rescaled intensity at the detector follows faithfully and monotonically the mechanical population. Note however that the approximate relationship (33) does not account for the detailed behaviour of the output field intensity, which is instead captured via our techniques.

Conclusions
The current developments in the control of harmonic oscillators at the micro-, nano-and atomic scale and in the technologies to couple them to travelling electromagnetic degrees of freedom [17,36], as well as the wealth of applications quantum information theory has envisaged for such interfaces, call for compact and general frameworks to handle inputoutput processes in a variety of applied settings. This work responds to such a need by delivering an algebraic description, in terms of Gaussian channels, of input-output relationships for general white noise and quadratic interactions. As demonstrated, the applicability of the method is broad while its predictions are directly observable, and we hence believe it to hold potential for further applications in the context of input-output quantum interfaces.

Appendix A. Derivation of the Gaussian Channels
Appendix A.1. Exponential Pulses.
Let us start with the explicit derivation of the channel associated to Eq. (11) . Combining with Eq. (10), and considering the j-th component, one haŝ where we have used that the primitive of a matrix exponential e Bt , in a domain where B is invertible, is given by B −1 e Bt , and that all the exponentials involved in our calculations vanish in the limit t → ∞. Integration by parts has been used to simplify the second term in Eq. (A.1); more specifically, we have integrated the term e (A+λ j 1)t and differentiated the term t 0 e −Asv in (s)ds. Eqs. (12)-(15) follow by calculating explicitly the second moments of the modesf j . Combining Eq. (A.2) with Eqs. (7) and (8), and defining δf j =f j − f j we have:

. Slow detectors and time-dependent spectral measurements
We now move on to the derivation of the Gaussian channel associated to Eq. (17) . Combining Eqs. (17) and (10) Considering the j-th component, we havê where N t,j = |2Reλ j | 1/2 [1 − e 2Re(λ j )t ] −1/2 , and we have again performed integration by parts on the second term in Eq. (A.4). The expression for the matrix Z of Eq. (19) follows easily from the first term in Eq. (A.5). To show how the matrix T [Eq. (20)] is derived, we shall now compute explicitly the second moments of the modesĝ j . Using Eqs. (7) and (8), performing explicitly the integrals of the form t 0 dt e Bt = B −1 (e Bt −1), we have From Eq. (A.6), one can derive easily Eqs. (20)- (26). The only integral which cannot be performed straightforwardly is given by As emphasized in the main text, one can verify explicitly that I verifies a Sylvester equation It is known that Eq. (A.8) has a unique solution if and only if A and −A † have no common eigenvalues [16]. Note that in our system this condition is automatically satisfied. In fact, the hypothesis of stability for the matrix A (that is, lim t→∞ e At = 0), implies that the eigenvalues of A have strictly negative real parts, hence those of −A † must have positive real parts. Consequently, the unique algebraic solution to Eq. (A.8) provides the result of the integral in Eq. (A.7).

Appendix B. Ion-cavity Hamiltonian
Our system is composed of a two-level ion, trapped inside an optical cavity and driven by an external laser beam. The internal levels of the ion are denoted as |e , |g , have frequency splitting ω 0 and their dynamics is conveniently described via a set of Pauli operators σ z , σ ± (In our notation, we have σ z = |e e| − |g g|, σ + = |e g|, σ − = |g e|). We assume that it is sufficient to consider a single mode of the electromagnetic field inside the optical cavity, with annihilation operator a 2 and frequency ω c . We also assume that the motion of the center of mass of the ion is relevant only along one axis, which we may choose as x without loss of generality. The total Hamiltonian, in a frame rotating with the frequency ω L of the driving laser, reads where a 1 is the annihilation operator for the ion motion along the x axis [x = x 0 (a 1 +a † 1 ), where x 0 is the ground state spread], ν the motional trapping frequency, δ = ω c − ω L , ∆ = ω 0 − ω L are respectively the detuning of the cavity and the ion internal transition with respect to the laser frequency, g 0 is the strength of the Jaynes-Cummings interaction between the internal levels of the ion and the cavity field, while Ω is the Rabi frequency of the driving laser, k being the projection of the light wavevector on the x axis (k ≤ ω L /c). Note that, as a consequence of its motion, the ion experiences variations in the phase of the driving field, while the strength of the Jaynes-Cummings interaction remains approximately constant. This is the case if, for example, the equilibrium position of the ion coincides with a maximum of the cavity field intensity. In the Lamb-Dicke regime k x 2 1, we can approximate e ikx 1 + ikx, hence where η = kx 0 is the Lamb-Dicke Parameter. As we have emphasized in the main text, we are interested in a large detuning regime in which the internal transition of the ion is only virtually excited (that is, ∆ g 0 , Ω, ν, δ). Therefore, we shall introduce only a small error by substituting where a 3 is a bosonic mode. In this way, the problem becomes tractable in our general language of Inputoutput Gaussian channels. The resulting Hamiltonian, up to a constant, is then given by Now, the quadratic part of the Hamiltonian in Eq. (B.3) coincides with the expression found in Eq. (28) , while the linear part can be included in the equations of motion by adding a constant to the vectorv in (see main text). In order for the bosonic treatment of the two-level system to be consistent, one has to check a posteriori that, at all times, the population of the mode a 3 is confined to the ground and first excited state with good approximation. Below we show estimates for the probability that the mode a 3 occupies excited states above the first [indicated as P (a † 3 a 3 > 1)], as a function of rescaled time, for the applications described in the main text (the parameters used are given in Table 1). We have fixed the initial motional squeezing to 20dB, which yields the largest occupation for the mode a 3 in the considered parameter range. Note how the results are consistent with the bosonic approximation of the two level system. This quantity has been calculated by explicitly integrating Eq. (1), using the full Hamiltonian of Eq. (B.3), and assuming a thermal distribution for the population of the mode a 3 , so that the only relevant parameter is a † 3 a 3 , which can calculated via standard input-output theory. One then has the estimate

Appendix C. Adiabatic elimination of the cavity field
Let us now take Eq. (29) as a starting point. Taking δ = ν, and assuming J ν, we can use a Rotating Wave Approximation [21,22] to neglect counter-rotating terms. Hence, we havê Correspondingly, the equations of motion for the bosonic modes reaḋ If the cavity is initially empty, and the coupling of the cavity field to the mechanical mode is weak, i.e. J κ 2 , we can expect that the mode a 2 will remain in its ground state with good approximation. In the Heisenberg picture, this amounts to takinġ a 2 −iνa 2 . Substituting in Eq. (C.3), we get Defining κ = 4J 2 /κ 2 , substituting Eq. (C.4) in Eq. (C.2) and considering the limit κ n th κ 1 (so that the environment associated to a in,1 can be neglected), we get Eq. (30). To get Eq. (31), one simply has to combine Eq. (C.4) with the boundary condition a out,2 = √ κ 2 a 2 − a in,2 .
Let us now show how to obtain Eq. (33). The bosonic mode of interest (see main text) has the explicit form where we have already taken into account the adiabatic elimination of the mode a 2 . In the limit where Γ is large compared to any other frequency in the problem, me may assume that the main contribution to the integral is achieved for t ∼ t. Hence, we may replace a 1 (t ) a 1 (t) and perform the integration explicitly, yielding where ξ in is a noise operator resulting from the integration of a in,2 , and we have approximated e −Γ/2t 0. By construction, the bosonic mode associated to ξ in is in the vacuum state. Therefore, Eq. (33) can be derived easily from Eq. (C.6). In addition, one has Combining appropriately Eqs. (33) and (C.7), the relationship between the statistics of arbitrary quadratures of g(t) and a 1 (t) can be derived.

Appendix D. Two-sided cavities, and losses in the output fields
Consider a cavity mode a, such that photons can leak out of both mirrors (left or right) with rates κ L , κ R respectively. In addition, the mode is subject to some Hamiltonian H. The corresponding Quantum Langevin equation readṡ where a in,L , a in,R are two uncorrelated and independent bosonic noise operators, with the usual properties. The output fields on the two sides of the cavity are defined via the boundary conditions: Now, let us consider the beam-splitter transformation where we have defined κ tot = κ L + κ R . We can now define the transformed input-output operators where we see that the cavity mode now couples only to the noise operator a in , while the mode ξ in is simply reflected. Our general theory as described in the main text can be applied to the effective output mode a out , however in this case the physically meaningful modes are associated to the outputs of the left and right mirror, according to a relation of the type a out,R (t) , (D.10) where for simplicity we have assumed the same temporal profile ϑ(t). Different temporal profiles can also be accounted for, but they are outside the scopes of this paper. Let us now show the relation between the physical modes f L , f R , and the mode that can be calculated with our theory. By inverting Eq. (D.6), and substituting in Eq. (D.10), one can easily see that where we see explicitly that a two-sided cavity is equivalent to a one-sided cavity followed by a beam-splitter (with the parameters given above). Note that the statistics of f can be calculated via the methods given in the main text, while that of ξ is uniquely determined by the noise operators a in,L , a in,R (Typically, ξ will be in the vacuum or in a thermal state). As a final remark, we note that the mathematics of a two-sided cavity described here can be also used to describe partial losses in the output field. Indeed, one can take a out,L to describe the accessible portion of the output field, while a out,R can be interpreted as the portion that is lost (e.g. due to scattering and absorption on the mirrors). Then, only f L describes the physically accessible bosonic mode, while the mode f R can be traced out.

Appendix E. Experimental parameters
The parameters of Table 1 can be obtained by considering the λ = 230nm transition of an Indium ion [27]. For this transition, the decay rate is κ 3 = 2π · 360kHz. Assuming that the ion equilibrium position coincides with a maximum of the cavity field, the Jaynes-Cummings coupling strength is given by [28] g 0 = 1 2 where c is the speed of light in vacuum, and V is the effective volume of the cavity field. For a Gaussian TEM 00 mode, we have [28] V L × π w 0 2 2 , where L is the length of the cavity and w 0 the waist of the mode at the cavity center. Assuming L = 1cm and w 0 = 6µm we get Let us briefly discuss the remaining parameters. The value κ 2 = 2π · 53kHz, as given in the main text, corresponds to a cavity finesse Finally, we assume that the ion trap is stabilized to provide an average heating time (n th κ 1 ) −1 = 6.6ms at T = 300K. This implies an average phonon number n th 1.2 · 10 6 , (E. 6) and an heating rate per phonon κ 1 2π · 2 · 10 −5 Hz. (E.7)