Input–output Gaussian channels: theory and application

Interfacing static and flying quantum degrees of freedom is a key step towards the realization 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⁵. In current experimental setups, a promising approach to this problem consists of coupling a quantum degree of freedom with cavity light and then adopting the light leaking out of the cavity as the flying quantum degree of freedom [1–5]. The latter may then be used to entangle distant trapped degrees of freedom by mixing on a beam splitter [6] or to achieve entanglement in a local site by the detection of the emerging light [7]. From the theoretical standpoint, this paradigm is well described by the seminal input–output formalism developed in the 1980s by Yurke, Collett and Gardiner [8–10].

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 [11–13] 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 characterization of their state will hence be provided in such a case. Let us emphasize 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 characterize. Furthermore, such a treatment is wholly independent of 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 optomechanical setups) 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.

We consider a localized system composed of n bosonic modes, whose canonical operators are grouped together in the column vector $\hat {\boldsymbol {v}}(t)=({a}_1,{a}^{\dagger }_1,\ldots ,{a}_n,{a}^{\dagger }_n)^\intercal$ and subject to the most general affine quantum Langevin equation

$$\begin{equation} \frac{\mathrm{d}\hat{\boldsymbol{v}}}{\mathrm{d}t} = A \hat{\boldsymbol{v}} + \hat{\boldsymbol{v}}_{\mathrm{in}}(t) , \end{equation} \tag{ 1 }$$

where A is the 'drift matrix', determined by the system's own Hamiltonian and linear coupling to the environment (see below), and

$$\begin{equation} \hat{\boldsymbol{v}}_{\mathrm{in}}(t) = K (a_{\mathrm{in},1}(t),a^{\dagger}_{\mathrm{in},1}(t),\ldots,{a}_{\mathrm{in,n}}(t),{a}^{\dagger}_{\mathrm{in,n}}(t))^\intercal + \bar{\boldsymbol{v}}_{\mathrm{in}}, \end{equation} \tag{ 2 }$$

$$\begin{equation} K=\mathrm{diag}(\sqrt{\kappa_1},\sqrt{\kappa_1},\ldots,\sqrt{\kappa_n},\sqrt{\kappa_n}). \end{equation} \tag{ 3 }$$

Here, κ_j quantifies the loss rate of the jth mode of the system, which interacts with the input operators $\hat {a}_{\mathrm {in},\mathrm {j}}$ and a^†_in,j, while $\bar {\boldsymbol {v}}_{\mathrm {in}}\in {\mathbbm C}^{2n}$ is just a constant vector, satisfying $\bar {v}_{\mathrm {in},2j-1}=\bar {v}_{\mathrm {in},2j}^{*}$, accounting for the possible linear driving of the system. Given a quadratic Hamiltonian

$$\begin{equation} \hat H=1/2\sum_{jk}H_{jk}\hat{v}_j\hat{v}_k+\sum_jV_j\hat{v}_j, \end{equation} \tag{ 4 }$$

with $H=H^\intercal$,V_2j−1 = V*_2j, A and $\bar {\boldsymbol {v}}_{\mathrm {in}}$ can be expressed as

$$\begin{equation} A=-\mathrm{i}\Sigma H-K^2/2, \end{equation} \tag{ 5 }$$

$$\begin{equation} \bar{\boldsymbol{v}}_{\mathrm{in}}=-\mathrm{i}\Sigma\boldsymbol{V}, \end{equation} \tag{ 6 }$$

where Σ is the commutation matrix $\Sigma _{jk}=[\hat {v}_j,\hat {v}_k]$. Throughout the paper, we shall only consider white Gaussian input noise, characterized by $\langle \hat {\boldsymbol {v}}_{\mathrm {in}}\rangle =\bar {\boldsymbol {v}}_{\mathrm {in}}$ and by the following correlation functions:

$$\begin{equation} \left\langle\{\hat{v}_{\mathrm{in},j}(t),\hat{v}^{\dagger}_{\mathrm{in},k}(t')\}\right\rangle - 2\bar{v}_{\mathrm{in},j}\bar{v}_{\mathrm{in},k}^{*}= ({\boldsymbol{\sigma}}_{\mathrm{in}})_{jk}\delta(t-t'), \end{equation} \tag{ 7 }$$

$$\begin{equation} \langle \hat{{v}}_{\mathrm{in},j}(t) \hat{{v}}^\dagger_k(0)\rangle-\bar{{v}}_{\mathrm{in},j}\langle\hat{{v}}^\dagger_k(0)\rangle=0, \end{equation} \tag{ 8 }$$

where σ_in is related to a physical CM, as it obeys K⁻¹σ_inK⁻¹ + iΣ ⩾ 0. The output field is determined by the boundary condition

$$\begin{equation} \hat{\boldsymbol{v}}_{\mathrm{out}}=K\hat{\boldsymbol{v}}-K^{-1}\hat{\boldsymbol{v}}_{\mathrm{in}} , \end{equation} \tag{ 9 }$$

where we deviate slightly from the standard definitions (see, e.g., [10]) to keep the formulae that will follow more compact. One can then insert the formal solution of equation (1) into (9) to obtain

$$\begin{equation} \hat{\boldsymbol{v}}_{\mathrm{out}}(t) = K \,\mathrm{e}^{At} \hat{\boldsymbol{v}}(0) + K \int_{0}^{t} \mathrm{e}^{A(t-s)} \hat{\boldsymbol{v}}_{\mathrm{in}}(s) \mathrm{d}s - K^{-1}\hat{\boldsymbol{v}}_{\mathrm{in}}(t) \, . \end{equation} \tag{ 10 }$$

Note that, in general, not all the components of the output field $\hat {\boldsymbol {v}}_{\mathrm {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.

In many applications, one is interested in output modes characterized by an exponential temporal profile (see figure 1 for an example). In general, this amounts to considering output modes of the form

$$\begin{equation} \skew6\hat{\boldsymbol{f}} = N \int_{0}^{\infty} \mathrm{e}^{\Lambda t} \hat{\boldsymbol{v}}_{\mathrm{out}}(t) \mathrm{d}t , \end{equation} \tag{ 11 }$$

where Λ = diag(λ₁,λ₂,...,λ_2n) = diag(μ₁,μ*₁,...,μ_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=\left |2\,\mathrm {Re}(\Lambda)\right |^{1/2}$ is a (diagonal) normalization matrix which guarantees that the $\skew6\hat{f}_j$ form a set of bosonic modes: $[\skew6\hat{f}_j,\skew6\hat{f}_k]=\Sigma _{jk}$.

Zoom In Zoom Out Reset image size

The matrix integrations involved in deriving the explicit solution for the vector of modes $\skew6\hat {\boldsymbol {f}}$ are best dealt with by considering each component $\skew6\hat{f}_j$ separately: inserting equation (10) into (11) yields (see appendix A.1)

$$\begin{equation} \fl \skew6\hat{f}_j=-\left[NK(A+\lambda_j \mathbbm{1})^{-1}\hat{\boldsymbol{v}}(0)+N(K(A+\lambda_j \mathbbm{1})^{-1}+K^{-1})\int_{0}^{\infty}\mathrm{e}^{\lambda_j t} \hat{\boldsymbol{v}}_{\mathrm{in}} \,\mathrm{d}t\right]_j, \end{equation} \tag{ 12 }$$

where $\mathbbm {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\rightarrow \infty }\exp (A t)=0$. The vector $\skew6\hat {\boldsymbol {f}}$ can be more compactly expressed by defining the elements of the matrix X and the vector $\hat {\boldsymbol {u}}$ as per

$$\begin{equation} X_{jk}=\left[NK(A+\lambda_j \mathbbm{1})^{-1}\right]_{jk} , \end{equation} \tag{ 13 }$$

$$\begin{equation} \hat{u}_j=\left[N(K(A+\lambda_j \mathbbm{1})^{-1}+K^{-1})\int_{0}^{\infty}\mathrm{e}^{\lambda_j t} \hat{\boldsymbol{v}}_{\mathrm{in}}\, \mathrm{d}t\right]_j . \end{equation} \tag{ 14 }$$

Note that the statistics of the vector $\hat {\boldsymbol {u}}$ are known, in that they can be reconstructed from equation (7). In particular, one has

$$\begin{equation} Y_{jk}=\langle\{\hat{u}_j,\hat{u}^{\dagger}_k\}\rangle-2\langle\hat{u}_j\rangle\langle\hat{u}^{\dagger}_k\rangle = -\frac{\left(X'\boldsymbol{\sigma}_{\mathrm{in}}{X'}^{\dagger}\right)_{jk}}{\lambda_{j}+\lambda_k^*} \, , \end{equation} \tag{ 15 }$$

where X' = X + NK⁻¹. Then, if σ(0) is the symmetrized covariance matrix (CM) of the initial state of the system ($\sigma _{jk}(0)=\langle \{\hat {v}_j(0),\hat {v}^{\dagger }_k(0)\}\rangle -2\langle \hat {v}_j(0)\rangle \langle \hat {v}^{\dagger }_k(0)\rangle$), one can determine the output CM $\boldsymbol {\sigma }_{\mathrm {out}}=\langle \{\skew6\hat{f}_j,\skew6\hat{f}^{\dagger }_k\}\rangle -2\langle \skew6\hat{f}_j\rangle \langle \skew6\hat{f}^{\dagger }_k\rangle$ from equation (12), obtaining

$$\begin{equation} \boldsymbol{\sigma}_{\mathrm{out}} = X \boldsymbol{\sigma}(0) X^{\dagger} + Y . \end{equation} \tag{ 16 }$$

Equation (16) defines a Gaussian quantum channel relating the system's CM at time t = 0 to the output CM σ_out [11, 13]. 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 [13], one has simply to transform all the equations by similarity with the unitary

$$\begin{eqnarray*} U=\frac{1}{\sqrt{2}}\bigoplus_{j=1}^{n}\left(\begin{array}{@{}cc@{}} 1 & 1\\ -i & i \end{array}\right) . \end{eqnarray*}$$

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 [14]

$$\begin{equation} \hat{\boldsymbol{g}}(t) = N_t \int_{0}^{t} \mathrm{e}^{\Lambda (t-t')} \hat{\boldsymbol{v}}_{\mathrm{out}}(t') \mathrm{d}t' , \end{equation} \tag{ 17 }$$

where the diagonal matrix Λ has the same form as above, and $N_t=\left |2\,\mathrm {Re}\,\Lambda \right |^{1/2}[\mathbbm {1}-\mathrm {e}^{2\,\mathrm {Re}(\Lambda) t}]^{-1/2}$ guarantees that the modes $\hat g_j(t)$ are bosonic. As anticipated, the modes of equation (17) have two main interpretations. On the one hand, they may describe continuous detection of the output fields with 'slow' detectors. In this case, |2 Re(λ_j)| quantifies the bandwidth (inverse of the response time) of the jth detector. On the other hand, they may be used to model spectral measurements with finite resolution, so that |2 Re(λ_j)| quantifies the spectral resolution of the jth frequency filter (figure 2).

Zoom In Zoom Out Reset image size

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)

$$\begin{equation} \boldsymbol{\sigma}_{\mathrm{out}} = Z \boldsymbol{\sigma}(0) Z^{\dagger} + T , \end{equation} \tag{ 18 }$$

with

$$\begin{equation} Z_{jk} = [NK(A-\lambda_j \mathbbm{1})^{-1}(\mathrm{e}^{At}-\mathbbm{1}\mathrm{e}^{\lambda_jt})]_{jk} , \end{equation} \tag{ 19 }$$

$$\begin{equation} T = FIF^{\dagger} - G + H + (J+J^{\dagger}) - (L+L^{\dagger}), \end{equation} \tag{ 20 }$$

and

$$\begin{equation} F_{jk} = [NK(A-\lambda_j \mathbbm{1})^{-1}]_{jk} , \end{equation} \tag{ 21 }$$

$$\begin{equation} G_{jk} = \frac{\left[(F+NK^{-1})\boldsymbol{\sigma}_{\mathrm{in}}(F+NK^{-1})^{\dagger}\right]_{jk}}{\lambda_j+\lambda_k^{*}} , \end{equation} \tag{ 22 }$$

$$\begin{equation} H_{jk} = {\mathrm{e}^{(\lambda_j+\lambda_k^{*})t}}G_{jk} , \end{equation} \tag{ 23 }$$

$$\begin{equation} J_{jk} = \left[F(A+\lambda_k^{*}\mathbbm{1})^{-1}\boldsymbol{\sigma}_{\mathrm{in}}(F+NK^{-1})^{\dagger}\right]_{jk} , \end{equation} \tag{ 24 }$$

$$\begin{equation} L_{jk} = \left[F\mathrm{e}^{(A+\lambda_k^{*}\mathbbm{1})t}(A+\lambda_k^{*}\mathbbm{1})^{-1}\boldsymbol{\sigma}_{\mathrm{in}}(F+NK^{-1})^{\dagger}\right]_{jk} , \end{equation} \tag{ 25 }$$

$$\begin{equation} I = \int_{0}^{t} \mathrm{e}^{At'}\boldsymbol{\sigma}_{\mathrm{in}}\mathrm{e}^{A^{\dagger}t'} \mathrm{d}t' . \end{equation} \tag{ 26 }$$

Note 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 [15].

Clearly, the matrices Z and T defining the Gaussian channel are now time dependent, as one should expect since the monitored output modes themselves depend on time in this approach.

4.1. Stationary spectrum

Let us briefly mention that, if one assumes the stability condition $\lim _{t\rightarrow \infty }\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:

$$\begin{equation} \boldsymbol{\sigma}_{\mathrm{out}} =FIF^\dagger - G + (J+J^{\dagger}) , \end{equation} \tag{ 27 }$$

where all quantities are intended in the limit t → ∞.

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 [1]. 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 [16, 18, 19]. 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 'phonon–photon 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 the 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 non-stationary regime.

We investigate in detail the case when 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 equation (29) can be taken as the starting point). Coupled ion–cavity systems have been realized experimentally in [35], and cavity cooling of a trapped ion was proposed in [17].

Let us consider an ion with motional trapping frequency ν along a chosen axis, corresponding to the mode a₁.⁶ The ion is trapped inside an optical cavity that sustains a mode a₂ of frequency ω_c, in turn coupled via a Jaynes–Cummings interaction of strength g₀ to an internal transition of the ion, described by a mode a₃ of frequency ω₀. 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)

$$\begin{eqnarray} &&\skew3\hat{H} = \nu a_1^{\dagger}a_1 + \delta a_2^{\dagger}a_2 + \Delta a_3^{\dagger}a_3 + g_0(a_2^{\dagger}a_3+a_2a_3^{\dagger}) -\eta\Omega (a_1+a_1^{\dagger}) (a_3+a_3^{\dagger}) , \end{eqnarray} \tag{ 28 }$$

where δ = (ω_c − ω_L), Δ = (ω₀ − ω_L) and η is the Lamb–Dicke parameter. As we shall see shortly, we are only interested in inducing virtual transitions of the ion, which will remain in its internal ground state with good approximation. Therefore, we only introduce a small error if we treat the mode a₃ 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 equation (3), where now n = 3), where κ₁ is the motional heating rate per phonon, κ₂ the cavity loss rate and κ₃ the spontaneous emission rate of the ion's internal transition. The thermal excitations of each environment are instead included in the matrix σ_in (see equation (7)). In particular, only the phononic environment is thermally excited in our problem, and we shall indicate its thermal occupancy by $\bar n_{\mathrm {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₀,Ω) [20, 21]. If the ion is initially in its internal ground state, one obtains the Hamiltonian

$$\begin{equation} \skew3\hat{H}' = \nu a_1^{\dagger}a_1 + \delta' a_2^{\dagger}a_2 + J(a_1+a_1^{\dagger})(a_2+a_2^{\dagger}) , \end{equation} \tag{ 29 }$$

where δ' = δ − g²₀/Δ and J = ηg₀Ω/Δ. As anticipated, a Hamiltonian of this form is also readily available in optomechanical setups, where it provides a reliable description of the dynamics under strong driving of the cavity field [16].

The direct coupling between the oscillator and the cavity input–output fields is finally obtained by choosing δ' = ν, and considering the limit $\nu \gg \kappa _2\gg J\gg \bar n_{\mathrm {th}}\kappa _1$. Then, neglecting the counter-rotating terms in equation (29), along with the adiabatic elimination of a₂ (see appendix C), yields the following input–output equations for the mechanical mode:

$$\begin{equation} \frac{\mathrm{d}a_1}{\mathrm{d}t}\simeq-\mathrm{i}\nu a_1-\frac{\kappa}{2}a_1-\mathrm{i}\sqrt{\kappa}a_{\mathrm{in},2}, \end{equation} \tag{ 30 }$$

$$\begin{equation} a_{\mathrm{out},\mathrm{2}}\simeq-\mathrm{i}\sqrt{\kappa}a_1+a_{\mathrm{in},2}, \end{equation} \tag{ 31 }$$

where κ = 4J²/κ₂ and we have also neglected the effects of the phononic environment ($\bar n_{\mathrm {th}}\kappa _1 \ll \kappa$ has to be checked a posteriori). As the photonic environment corresponding to a_in,2 is effectively at zero temperature, equation (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 being 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 in studying. However, the general analysis of sections 3 and 4 allows us to identify and treat exactly output modes that are suitable for studying 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 equations (16) and (18) and obtain exact quantitative results that are reliable even when the approximations leading to equations (30) and (31) are justified only partially. In what follows, we shall consider an ion–cavity system characterized by typical experimental parameters, which we summarize in table 1.

Table 1. Adopted values of dynamical parameters. These are obtained by considering the 230 nm transition of an indium ion [26], trapped at the centre of a Fabry–Perot cavity of length L ≃ 1 cm, waist w ≃ 6 μm and average photon lifetime (κ₂)⁻¹ = 3 μs. We have considered a stabilized trap with average heating time $(\bar n_{\mathrm {th}}\kappa _1)^{-1}=6.6\,\mathrm {ms}$ at room temperature (T = 300 K). The effective cooling rate shown above corresponds to a temporal length (κ)⁻¹ ≃  83 μs for the mode $\skew6\hat f$. The 50 MHz bandwidth of the detector corresponds to an average response time Γ⁻¹ ≃ 3.2 ns. For further details, see appendix E.

Motional frequency of the ion	ν=2π×5 MHz
Jaynes–Cummings coupling strength	g0=2π×0.62 MHz
Laser Rabi frequency	Ω=2π×1 MHz
Laser-internal transition detuning	Δ=2π×10 MHz
Lamb–Dicke parameter	η=0.08
Cavity decay rate	κ2=2π×53 kHz
Spontaneous emission rate	κ3=2π×360 kHz
Phononic heating rate	$\bar n_{\mathrm {th}}\kappa _1=2\pi \times 24\,\mathrm {Hz}$
Effective cooling rate	κ≃2π×19 kHz
Detector's bandwidth	Γ=2π×50 MHz

5.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 input–output theory being applied to equations (30) and (31), yielding [25]

$$\begin{equation} f=\sqrt{\kappa}\int_0^\infty\mathrm{d}t\,\mathrm{e}^{(\mathrm{i}\nu-\kappa/2)t}a_{\mathrm{out},2}(t)\simeq-\mathrm{i}a_1(0), \end{equation} \tag{ 32 }$$

which means that the cavity emits a pulse of effective duration τ ∼ 1/κ, whose quantum state is approximately the same as the oscillator initial state. Note that, as we shall verify in the following, equation (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 equation (32) corresponds to modes of the form (11), with μ₂ = iν − κ/2. Note that the parameters μ₁,μ₃ 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 when the mode a₁ is initially in a pure squeezed state, and investigate quantitatively the amount of squeezing transferred to the travelling mode f. The reliable realization of such squeezing transfer constitutes a basic requirement for continuous variable quantum communication and information processing. Moreover, in the medium term this setup 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 [22, 23, 28] or through the internal degrees of freedom [29–31], while in optomechanical systems, a high degree of squeezing could be obtained via indirect position measurements [32].

The squeezing of a Gaussian state with CM σ has been quantified, in dB, as $\max [0,-10\log _{10}(\sigma _{1}^{\uparrow })]$, where σ^↑₁ 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, characterized by the experimental parameters shown in table 1 (see also appendix E). Our quantitative findings are summarized in figure 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 9 dB after an initial linear rise, indicating a complete breakdown of the approximated relation f ≃ − ia₁(0) for high input squeezing.

Zoom In Zoom Out Reset image size

This is also reflected in the purity Tr(ϱ²) of the output state ϱ which, being our state Gaussian, can be determined as 1/(det σ_out)^1/2 [33]. Figure 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.

5.2. Entanglement generator

5.3. Continuous detection of the mechanical motion

Equation (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 [25].

To formulate the problem in terms of an input–output Gaussian channel, we consider modes of the form (17), with μ₂ = iν − Γ/2, Γ being the bandwidth of the detector. Again, the parameters μ₁,μ₃ can be set to any value, as we are only interested in the detection of the mode associated with a_out,2, which we shall denote by g(t). From equations (17) and (31), assuming that Γ is large compared to the other frequencies in the problem, we obtain (see appendix C)

$$\begin{equation} \langle{g(t)^\dagger g(t)}\rangle\simeq\frac{4\kappa}{\Gamma}\langle{a_1(t)^\dagger a_1(t)}\rangle. \end{equation} \tag{ 33 }$$

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 figure 3(d), we compare inferred and actual values of 〈a^†₁a₁〉, as a function of time, for an indium ion with an initial motional squeezing of 20 dB. 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.

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 [16, 35], as well as the wealth of applications that quantum information theory has envisaged for such interfaces, call for compact and general frameworks to handle input–output 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.

In fond memory of Wolfgang Lange. We thank V Giovannetti, G Milburn, A Bayat, S Bose, M S Kim and M Tame for useful discussions. We acknowledge financial support from the Royal Society, Wolfson Foundation, UK EPSRC, Qatar National Research Fund (NPRP 4-554-1-084), EU STREP HIP and PICC, EU integrated project QESSENCE and the Alexander von Humboldt Foundation.

A.1. Exponential pulses

Let us start with the explicit derivation of the channel associated with equation (11). Combining with equation (10) and considering the jth component, one has

$$\begin{eqnarray} &&\fl \skew6\hat f_j=|2\,\mathrm{Re}\,\lambda_j|^{1/2}\left[K \int_0^\infty\mathrm{d}t\;\mathrm{e}^{(A+\lambda_j\mathbbm{1})t} \hat{\boldsymbol{v}}(0) + K\int_0^\infty\mathrm{d}t\;\mathrm{e}^{(A+\lambda_j\mathbbm{1})t} \int_{0}^{t} \mathrm{e}^{-As} \hat{\boldsymbol{v}}_{\mathrm{in}}(s) \mathrm{d}s \right.\nonumber\\ &&\left.- K^{-1}\int_0^\infty\mathrm{d}t\;\mathrm{e}^{\lambda_j t}\hat{\boldsymbol{v}}_{\mathrm{in}}(t) \right]_j \end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray} &=&-|2\,\mathrm{Re}\,\lambda_j|^{1/2}\left[\vphantom{\frac{\int}{1}}K(A+\lambda_j\mathbbm{1})^{-1} \hat{\boldsymbol{v}}(0)\right.\nonumber\\ & &\left.+(K(A+\lambda_j\mathbbm{1})^{-1}+K^{-1})\int_0^\infty\mathrm{d}t\;\mathrm{e}^{\lambda_j t}\hat{\boldsymbol{v}}_{\mathrm{in}}(t)\right]_j, \end{eqnarray} \noindent \tag{ A.2 }$$

where we have used that the primitive of a matrix exponential e^Bt, in a domain where B is invertible, is given by B⁻¹ 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 equation (A.1); more specifically, we have integrated the term $\mathrm {e}^{(A+\lambda _j\mathbbm {1})t}$ and differentiated the term $\int _{0}^{t} \mathrm {e}^{-As} \hat {\boldsymbol {v}}_{\mathrm {in}}(s) \mathrm {d}s$. Equations (12)–(15) follow by calculating explicitly the second moments of the modes $\skew6\hat f_j$. Combining equation (A.2) with equations (7) and (8) and defining $\delta \skew6\hat f_j=\skew6\hat f_j-\langle {\skew6\hat f_j}\rangle$, we have

$$\begin{eqnarray} &&\fl \langle\{\delta\skew6\hat{f}_j,\delta\skew6\hat{f}^{\dagger}_k\}\rangle=|2\,\mathrm{Re}\,\lambda_j|^{1/2}|2\,\mathrm{Re}\,\lambda_k|^{1/2}\left[\vphantom{\frac{\int}{1}}K(A+\lambda_j\mathbbm{1})^{-1}\boldsymbol{\sigma}(0)(K(A+\lambda_k\mathbbm{1})^{-1})^\dagger\right.\nonumber\\ &&\left.-(K(A+\lambda_j\mathbbm{1})^{-1}+K^{-1})\frac{\boldsymbol{\sigma}_{\mathrm{in}}}{\lambda_j+\lambda_k^*}(K(A+\lambda_k\mathbbm{1})^{-1}+K^{-1})^\dagger\right]_{jk}. \end{eqnarray} \tag{ A.3 }$$

A.2. Slow detectors and time-dependent spectral measurements

We now move on to the derivation of the Gaussian channel associated with equation (17). Combining equations (17) and (10) considering the jth component, we have

$$\begin{eqnarray} &&\fl \hat g_j(t)=N_{t,j}\,\mathrm{e}^{\lambda_j t}\left[K\int_0^t\mathrm{d}t'\;\mathrm{e}^{(A-\lambda_j\mathbbm{1})t'}\hat{\boldsymbol{v}}(0)+K\int_0^t\mathrm{d}t'\;\mathrm{e}^{(A-\lambda_j\mathbbm{1})t'}\int_0^{t'}\mathrm{d}s\;\mathrm{e}^{-As}\hat{\boldsymbol{v}}_{\mathrm{in}}(s)\right.\nonumber\\ &&\left.-K^{-1}\int_0^t\mathrm{d}t'\;\mathrm{e}^{-\lambda_jt'}\hat{\boldsymbol{v}}_{\mathrm{in}}(t')\vphantom{K\int_0^t\mathrm{d}t'\;\mathrm{e}^{(A-\lambda_j\mathbbm{1})t'}\hat{\boldsymbol{v}}(0)+K\int_0^t\mathrm{d}t'\;\mathrm{e}^{(A-\lambda_j\mathbbm{1})t'}\int_0^{t'}\mathrm{d}s\;\mathrm{e}^{-As}\hat{\boldsymbol{v}}_{\mathrm{in}}(s)}\right]_j \end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray} &&\fl\hphantom{\hat g_j(t)} =N_{t,j}\left[\vphantom{\frac{\int}{1}}K(A-\lambda_j\mathbbm{1})^{-1}(\mathrm{e}^{At}-\mathrm{e}^{\lambda_jt})\hat{\boldsymbol{v}}(0)+K(A-\lambda_j\mathbbm{1})^{-1}\right.\nonumber\\ &&\times \left.\int_0^t\mathrm{d}t'\;\mathrm{e}^{A(t-t')}\hat{\boldsymbol{v}}_{\mathrm{in}}(t')-(K(A-\lambda_j\mathbbm{1})^{-1}+K^{-1})\int_0^t\mathrm{d}t'\mathrm{e}^{\lambda_j(t-t')}\hat{\boldsymbol{v}}_{\mathrm{in}}(t')\right]_j, \end{eqnarray} \noindent \tag{ A.5 }$$

where $N_{t,j}= |2\,\mathrm {Re}\,\lambda _j |^{1/2}[\mathbbm {1}-\mathrm {e}^{2\,\mathrm {Re}(\lambda _j) t}]^{-1/2}$, and we have again performed integration by parts on the second term in equation (A.4). The expression for the matrix Z of equation (19) follows easily from the first term in equation (A.5). To show how the matrix T (equation (20)) is derived, we shall now compute explicitly the second moments of the modes $\hat g_j$. Using equations (7) and (8), performing explicitly the integrals of the form $\int _0^t\mathrm {d}t'\,\mathrm {e}^{Bt'}=B^{-1}(\mathrm {e}^{Bt}-\mathbbm {1})$, we have

$$\begin{eqnarray} &&\fl \langle\{ \delta \hat g_j,\delta \hat g_k^\dagger\}\rangle =[Z\boldsymbol{\sigma}(0)Z^\dagger]_{jk}+N_{t,j}N_{t,k}\left[\vphantom{\sum}K(A-\lambda_j\mathbbm{1})^{-1}\right.\nonumber\\ &&\left.\times \int_0^t\mathrm{d}t'\mathrm{e}^{At'}\boldsymbol{\sigma}_{\mathrm{in}}\mathrm{e}^{A^\dagger t'}(K(A-\lambda_k\mathbbm{1})^{-1})^\dagger\right.\nonumber\\ &&+\left.\frac{\mathrm{e}^{(\lambda_j+\lambda_k^*)t}-1}{\lambda_j+\lambda_k^*}(K(A-\lambda_j\mathbbm{1})^{-1}+K^{-1})\boldsymbol{\sigma}_{\mathrm{in}}(K(A-\lambda_k\mathbbm{1})^{-1}+K^{-1})^\dagger\right.\nonumber\\ &&\left.-K(A+\lambda_j\mathbbm{1})^{-1}(A+\lambda_k^*\mathbbm{1})^{-1}(\mathrm{e}^{(A+\lambda_k^*\mathbbm{1})t}-\mathbbm{1})\right.\nonumber\\ &&\times\left. \boldsymbol{\sigma}_{\mathrm{in}}(K(A-\lambda_k\mathbbm{1})^{-1}+K^{-1})^\dagger-(K(A-\lambda_j\mathbbm{1})^{-1}+K^{-1})\right.\nonumber\\ &&\left.\times \boldsymbol{\sigma}_{\mathrm{in}}(A^\dagger+\lambda_j\mathbbm{1})^{-1}(\mathrm{e}^{(A^\dagger+\lambda_j\mathbbm{1})t}-\mathbbm{1})(K(A-\lambda_k\mathbbm{1})^{-1})^\dagger\right]_{jk}. \end{eqnarray} \tag{ A.6 }$$

From equation (A.6), one can derive easily equations (20)–(26). The only integral which cannot be performed straightforwardly is given by

$$\begin{equation} I=\int_0^t\mathrm{d}t'\,\mathrm{e}^{At'}\boldsymbol{\sigma}_{\mathrm{in}}\,\mathrm{e}^{A^\dagger t'}. \end{equation} \noindent \tag{ A.7 }$$

As emphasized in the main text, one can verify explicitly that I verifies a Sylvester equation

$$\begin{equation} AI+IA^\dagger=\int_0^t\mathrm{d}t'\;\frac{\mathrm{d}}{\mathrm{d}t'}\left(\mathrm{e}^{At'}\boldsymbol{\sigma}_{\mathrm{in}}\mathrm{e}^{A^\dagger t'}\right)=\mathrm{e}^{At}\boldsymbol{\sigma}_{\mathrm{in}}\mathrm{e}^{A^{\dagger}t}-\boldsymbol{\sigma}_{\mathrm{in}}. \end{equation} \noindent \tag{ A.8 }$$

It is known that equation (A.8) has a unique solution if and only if A and −A^† have no common eigenvalues [15]. 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 equation (A.8) provides the result of the integral in equation (A.7).

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 by $\left \vert e \right \rangle , \left \vert g \right \rangle$, have frequency splitting ω₀ 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₂ and frequency ω_c. We also assume that the motion of the centre 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

$$\begin{eqnarray} &&\fl H_{\mathrm{tot}}=\nu a_1^\dagger a_1+\delta a_2^\dagger a_2+\frac{\Delta}{2}\sigma_z+g_0(a_2^\dagger \sigma^-+a_2\sigma^+)+\mathrm{i}\Omega(\sigma^+\mathrm{e}^{\mathrm{i}kx}-\sigma^-\mathrm{e}^{-\mathrm{i}kx}), \end{eqnarray} \tag{ B.1 }$$

where a₁ is the annihilation operator for the ion motion along the x-axis (x = x₀(a₁ + a^†₁), where x₀ is the ground state spread), ν the motional trapping frequency, δ = ω_c − ω_L, Δ = ω₀ − ω_L are, respectively, the detuning of the cavity and the ion internal transition with respect to the laser frequency, g₀ 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\sqrt {\langle {x^2}\rangle }\ll 1$, we can approximate e^ikx ≃ 1 + ikx, hence

$$\begin{eqnarray} &&\fl H_{\mathrm{tot}} \simeq \nu a_1^\dagger a_1+\delta a_2^\dagger a_2+\frac{\Delta}{2}\sigma_z+g_0(a_2^\dagger \sigma^-+a_2\sigma^+)+\mathrm{i}\Omega(\sigma^+-\sigma^-) -\eta\Omega(\sigma^++\sigma^-)(a_1+a_1^\dagger), \end{eqnarray} \tag{ B.2 }$$

where η = kx₀ 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₀,Ω,ν,δ). Therefore, we shall introduce only a small error by substituting σ_z ≃ 2a^†₃a₃ − 1,σ⁺ ≃ a^†₃,σ⁻ ≃ a₃, where a₃ is a bosonic mode. In this way, the problem becomes tractable in our general language of input–output Gaussian channels. The resulting Hamiltonian, up to a constant, is then given by

$$\begin{eqnarray} &&\fl H_{\mathrm{tot}}\simeq \nu a_1^\dagger a_1+\delta a_2^\dagger a_2+\Delta a_3^\dagger a_3+g_0(a_2^\dagger a_3+a_2a_3^\dagger)+\mathrm{i}\Omega(a_3-a_3^\dagger) -\eta\Omega(a_3+a_3^\dagger)(a_1+a_1^\dagger). \end{eqnarray} \tag{ B.3 }$$

Now, the quadratic part of the Hamiltonian in equation (B.3) coincides with the expression found in equation (28), while the linear part can be included in the equations of motion by adding a constant to the vector $\bar {\boldsymbol {v}}_{\mathrm {in}}$ (see the 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₃ is confined to the ground and first excited states with a good approximation. Below we show estimates for the probability that the mode a₃ occupies excited states above the first (indicated as P(a^†₃a₃ > 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 20 dB, which yields the largest occupation for the mode a₃ 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 equation (1), using the full Hamiltonian of equation (B.3), and assuming a thermal distribution for the population of the mode a₃, so that the only relevant parameter is 〈a^†₃a₃〉, which can calculated via standard input–output theory. One then has the estimate

$$\begin{equation} P(a_3^\dagger a_3>1)\sim\left(\frac{\langle a_3^\dagger a_3\rangle}{1+\langle a_3^\dagger a_3\rangle}\right)^2. \end{equation} \tag{ B.4 }$$

Let us now take equation (29) as a starting point. Taking δ' = ν and assuming J ≪ ν, we can use a rotating wave approximation [20, 21] to neglect counter-rotating terms. Hence, we have

$$\begin{equation} \skew3\hat{H}' \simeq \nu a_1^{\dagger}a_1 + \nu a_2^{\dagger}a_2 + J(a_1a_2^\dagger+a_2a_1^{\dagger}) . \end{equation} \tag{ C.1 }$$

Correspondingly, the equations of motion for the bosonic modes read

$$\begin{equation} \dot a_1=-\mathrm{i}\nu a_1-\frac{\kappa_1}{2}a_1-\mathrm{i}Ja_2+\sqrt{\kappa_1}a_{\mathrm{in},1}, \end{equation} \tag{ C.2 }$$

$$\begin{equation} \dot a_2=-\mathrm{i}\nu a_2-\frac{\kappa_2}{2}a_2-\mathrm{i}Ja_1+\sqrt{\kappa_2}a_{\mathrm{in},2}. \end{equation} \tag{ C.3 }$$

If the cavity is initially empty and the coupling of the cavity field to the mechanical mode is weak, i.e. J ≪ κ₂, we can expect that the mode a₂ will remain in its ground state with good approximation. In the Heisenberg picture, this amounts to taking $\dot a_2\simeq -\mathrm {i}\nu a_2$. Substituting in equation (C.3), we obtain

$$\begin{equation} a_2\simeq-\mathrm{i}\frac{2J}{\kappa_2}a_1+\frac{2}{\sqrt{\kappa_2}}a_{\mathrm{in},2}. \end{equation} \tag{ C.4 }$$

Defining κ = 4J²/κ₂, substituting equation (C.4) in equation (C.2) and considering the limit $\kappa \gg \bar n_{\mathrm {th}}\kappa _1$ (so that the environment associated with a_in,1 can be neglected), we obtain equation (30). To obtain equation (31), one simply has to combine equation (C.4) with the boundary condition $a_{\mathrm {out},2}=\sqrt {\kappa _2}a_2-a_{\mathrm {in},2}$.

Let us now show how to obtain equation (33). The bosonic mode of interest (see the main text) has the explicit form

$$\begin{equation} g(t)\simeq\frac{\Gamma^{1/2}}{(1-\mathrm{e}^{-\Gamma t})}\int_0^t\mathrm{d}t'\;\mathrm{e}^{(\mathrm{i}\nu-\Gamma/2)(t-t')}(-\mathrm{i}\sqrt{\kappa}a_1(t')-a_{\mathrm{in},2}(t')), \end{equation} \tag{ C.5 }$$

where we have already taken into account the adiabatic elimination of the mode a₂. In the limit where Γ is large compared to any other frequency in the problem, we may assume that the main contribution to the integral is achieved for t' ∼ t. Hence, we may replace a₁(t') ≃ a₁(t) and perform the integration explicitly, yielding

$$\begin{equation} g(t)\simeq2\mathrm{i}\sqrt{\frac{\kappa}{\Gamma}}a_1(t)+\xi_{\mathrm{in}}, \end{equation} \tag{ C.6 }$$

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 with ξ_in is in the vacuum state. Therefore, equation (33) can be derived easily from equation (C.6). In addition, one has

$$\begin{equation} \langle g(t)^2\rangle\simeq\frac{4\kappa}{\Gamma}\langle a_1(t)^2\rangle. \end{equation} \tag{ C.7 }$$

Combining appropriately equations (33) and (C.7), the relationship between the statistics of arbitrary quadratures of g(t) and a₁(t) can be derived.

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 reads

$$\begin{equation} \dot a=\mathrm{i}[H,a]-\frac{\kappa_{\mathrm{L}}+\kappa_{\mathrm{R}}}{2}a+\sqrt{\kappa_{\mathrm{L}}}a_{\mathrm{in}, \mathrm{L}}+\sqrt{\kappa_{\mathrm{R}}}a_{\mathrm{in},\mathrm{R}}, \end{equation} \noindent \tag{ D.1 }$$

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

$$\begin{equation} a_{\mathrm{ out},\mathrm{L}}=\sqrt{\kappa_{\mathrm{L}}}a-a_{\mathrm{in}, \mathrm{L}}, \end{equation} \tag{ D.2 }$$

$$\begin{equation} a_{\mathrm{ out},\mathrm{R}}=\sqrt{\kappa_{\mathrm{R}}}a-a_{\mathrm{in},\mathrm{R}}. \end{equation} \tag{ D.3 }$$

Now, let us consider the beam-splitter transformation

$$\begin{eqnarray} BS=\frac{1}{\sqrt{\kappa_{\mathrm{tot}}}}\left(\begin{array}{@{}cc@{}} \sqrt{\kappa_{\mathrm{L}}} & \sqrt{\kappa_{\mathrm{R}}}\\ -\sqrt{\kappa_{\mathrm{R}}} & \sqrt{\kappa_{\mathrm{L}}} \end{array}\right), \end{eqnarray} \tag{ D.4 }$$

where we have defined κ_tot = κ_L + κ_R. We can now define the transformed input–output operators

$$\begin{equation} \left(\begin{array}{@{}c@{}} a_{\mathrm{in}}\\ \xi_{\mathrm{in}} \end{array}\right)=BS\left(\begin{array}{@{}c@{}} a_{\mathrm{in}, \mathrm{L}}\\ a_{\mathrm{in},\mathrm{R}} \end{array}\right), \end{equation} \tag{ D.5 }$$

$$\begin{equation} \left(\begin{array}{@{}c@{}} a_{\mathrm{out}}\\ \xi_{\mathrm{out}} \end{array}\right)=BS\left(\begin{array}{@{}c@{}} a_{\mathrm{ out},\mathrm{L}}\\ a_{\mathrm{ out},\mathrm{R}} \end{array}\right). \end{equation} \tag{ D.6 }$$

In terms of these, the equation of motion for a and the boundary conditions for the input–output fields can be rewritten as

$$\begin{equation} \dot a=\mathrm{i}[H,a]-\frac{\kappa_{\mathrm{tot}}}{2}a+\sqrt{\kappa_{\mathrm{tot}}}a_{\mathrm{in}}, \end{equation} \tag{ D.7 }$$

$$\begin{equation} a_{\mathrm{out}}=\sqrt{\kappa_{\mathrm{tot}}}a-a_{\mathrm{in}}, \end{equation} \tag{ D.8 }$$

$$\begin{equation} \xi_{\mathrm{out}}=-\xi_{\mathrm{in}}, \end{equation} \noindent \tag{ D.9 }$$

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 with the outputs of the left and right mirrors, according to a relation of the type

$$\begin{equation} \left(\begin{array}{@{}c@{}} f_{\mathrm{L}}\\ f_{\mathrm{R}} \end{array}\right)=\int\mathrm{d}t\; \vartheta(t)\left(\begin{array}{@{}c@{}} a_{\mathrm{ out},\mathrm{L}}(t)\\ a_{\mathrm{ out},\mathrm{R}}(t) \end{array}\right), \end{equation} \noindent \tag{ 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 scope 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 equation (D.6) and substituting in equation (D.10), one can easily see that

$$\begin{equation} \left(\!\begin{array}{@{}c@{}} f_{\mathrm{L}}\\ f_{\mathrm{R}} \end{array}\!\right)=BS^\dagger \left(\!\begin{array}{@{}c@{}} f\\ \xi \end{array}\!\right), \end{equation} \tag{ D.11 }$$

$$\begin{equation} \left(\!\begin{array}{@{}c@{}} f\\ \xi \end{array}\right)=\int\mathrm{d}t\; \vartheta(t) \left(\begin{array}{@{}c@{}} a_{\mathrm{out}}(t)\\ \xi_{\mathrm{out}}(t) \end{array}\!\right), \end{equation} \tag{ D.12 }$$

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.

The parameters of table 1 can be obtained by considering the λ = 230 nm transition of an indium ion [26]. For this transition, the decay rate is κ₃ = 2π × 360 kHz. Assuming that the ion equilibrium position coincides with a maximum of the cavity field, the Jaynes–Cummings coupling strength is given as [27]

$$\begin{equation} g_0=\frac{1}{2}\sqrt{\frac{3\lambda^2 c\kappa_3}{2\pi V}}, \end{equation} \noindent \tag{ E.1 }$$

where c is the speed of light in vacuum, and V is the effective volume of the cavity field. For a Gaussian TEM₀₀ mode, we have [27]

$$\begin{equation} V\simeq L\times\pi\left(\frac{w_0}{2}\right)^2, \end{equation} \noindent \tag{ E.2 }$$

where L is the length of the cavity and w₀ the waist of the mode at the cavity centre. Assuming that L = 1 cm and w₀ = 6 μm, we obtain

$$\begin{equation} g_0\simeq2\pi\times 0.62\,\mathrm{MHz}. \end{equation} \noindent \tag{ E.3 }$$

Let us briefly discuss the remaining parameters. The value κ₂ = 2π × 53 kHz, as given in the main text, corresponds to a cavity finesse

$$\begin{equation} \mathcal F=\frac{\pi c}{L\kappa_2}\simeq 2.8\times 10^5. \end{equation} \noindent \tag{ E.4 }$$

A motional trap frequency of ν = 2π × 5 MHz, combined with a mass of m ≃ 114 au for the indium ion, yields the Lamb–Dicke parameter (in our units, ℏ = 1)

$$\begin{equation} \eta=\frac{2\pi x_0}{\lambda}=\frac{2\pi}{\lambda}\sqrt{\frac{1}{2m\nu}}\simeq0.081. \end{equation} \tag{ E.5 }$$

Finally, we assume that the ion trap is stabilized to provide an average heating time $(\bar n_{\mathrm {th}}\kappa _1)^{-1}=6.6\,\mathrm {ms}$ at T = 300 K. This implies an average phonon number

$$\begin{equation} \bar n_{\mathrm{th}}\simeq 1.2\times 10^6 \end{equation} \tag{ E.6 }$$

and a heating rate per phonon

$$\begin{equation} \kappa_1\simeq 2\pi\times 2\times 10^{-5}\,\mathrm{Hz}. \end{equation} \tag{ E.7 }$$