Quantum State Reconstruction of an Oscillator Network in an Optomechanical Setting

We introduce a scheme to reconstruct an arbitrary quantum state of a mechanical oscillator network. We assume that a single element of the network is coupled to a cavity field via a linearized optomechanical interaction, whose time dependence is controlled by a classical driving field. By designing a suitable interaction profile, we show how the statistics of an arbitrary mechanical quadrature can be encoded in the cavity field, which can then be measured. We discuss the important special case of Gaussian state reconstruction, and study numerically the effectiveness of our scheme for a finite number of measurements. Finally, we speculate on possible routes to extend our ideas to the regime of single-photon optomechanics.


I. INTRODUCTION
Quantum optomechanics exploits radiation pressure to couple photons and mechanical oscillators. The field has progressed significantly in the last decade and is now entering a promising stage where the observation of quantum effects in macroscopic objects appears to be within grasp [1][2][3]. Substantial theoretical and experimental effort has been put into the preparation of mechanical systems, typically consisting of vibrating mirrors or membranes, in interesting non-classical states. In such a context, an important question arises: how can we verify that the mechanical state prepared in an experiment is indeed the desired one? The design of successful strategies to achieve such verifications requires the experimental estimation of the density operator of a mechanical system. However, it is well known that the full information encoded in the density operator cannot be accessed through the measurement of a single observable. One must instead collect the measurement statistics of several distinct observables, a task which requires access to many copies of the quantum system of interest. These could be obtained, for instance, by repeating the same experiment with the same initial conditions. By post-processing the outcomes of such measurements, an experimentalist can estimate the density operator via techniques known as quantum tomography and quantum state reconstruction [4,5]. Perhaps the best known example in this context is the reconstruction of the Wigner function of an oscillator (which brings about an amount of information equivalent to that of the density operator) through a Radon transform of the quadrature probability densities [6].
In an optomechanical setting various approaches to quantum state reconstruction have been explored in the literature [2,7,8], in particular employing weak or quan- * Electronic address: dmoore32@qub.ac.uk † Electronic address: tommaso.tufarelli@nottingham.ac.uk ‡ Electronic address: m.paternostro@qub.ac.uk § Electronic address: a.ferraro@qub.ac.uk tum non-demolition measurements of mechanical quadratures [7,[9][10][11]. Other techniques that have been put forward include the use of short laser pulses to prepare and read out the mechanical state [12][13][14], the exploitation of a detuned driving field [15], and the measurement of the phonon-number operator [16]. The possibility of a precise readout has also opened the way towards feedback cooling of a mechanical oscillator [17]. Very recently, high-efficiency state estimation of a mechanical oscillator through techniques based on Kalman filtering have been implemented experimentally, paving the way to the real-time reconstruction of mechanical state-space configurations, and their quantum-limited control [18]. Besides optomechanical systems comprising a single mechanical oscillator, one can envisage situations in which multiple mechanical oscillators interact in a small Only the mechanical mode b1 is directly coupled to the radiation field. The driving field (red arrow) can be modulated in order to realize the time-varying linearized radiationpressure coupling needed to reconstruct a selected quadrature of the oscillator network. Only the field aout leaking from the cavity is eventually measured via homodyne detection. Repeating the protocol for a sufficient set of quadratures, the state of the entire mechanical network can be reconstructed.
network and operate close to the quantum regime [19][20][21] (Fig. 1). The latter would implement systems of interacting bosons which are of paramount interest in a variety of contexts -including quantum thermodynamics [22] and quantum simulators [23], where they represent an ideal playground to test fundamental issues like equilibration [24], heat transport [25][26][27][28][29], the definition of temperature [30][31][32], and the universal scaling of ground-state entanglement [33,34]. In addition, interacting quantum oscillators have been proposed as a valid route towards quantum computation [35][36][37][38]. Quantum state reconstruction is instrumental in all these settings. Outside the optomechanical domain, the reconstruction of an oscillator network can be accomplished using a two-level system interacting with one node of the network -provided the coupling constant is time-dependent and can be controlled by the experimenter [39,40]. However, apart from measuring each oscillator individually, to the best of our knowledge no method has been proposed for the efficient readout of the quantum state of an oscillator network in an optomechanical setting. In this paper we propose a protocol of quantum state reconstruction for the mechanical portion of a generalized optomechanical system, featuring a single high-quality cavity mode coupled to a network of mechanical oscillators (see Fig. 1). Our protocol relies on the so-called linearized radiation pressure interaction, and exploits measurements on the accessible output modes of the optical cavity (rather than the mechanical modes of the oscillator network which are typically challenging to measure directly). By controlling in time the interaction strength, we show that it is possible to encode information about any mechanical quadrature in the cavity light, which can then be measured through the output fields leaking out of the system. Specifically, we discuss how an arbitrary moment of the selected quadrature can be estimated via appropriate light-quadrature measurements, followed by the inversion of a linear system of equations. Our scheme shares an important advantage with Ref. [40]: it requires minimal access to the oscillator network, in that it can be probed through interaction with just one of its elements.
The paper is organised as follows. In Section II we present our state reconstruction scheme applied to a single mechanical oscillator, a simple example which provides a gentle introduction to the more technical case of a network. We start by showing how the dynamics of interest can be solved analytically, then discuss how the interaction profile can be designed to encode a chosen mechanical quadrature in the cavity light mode. Finally, we conclude the section by explaining how to measure indirectly the light mode through collection of the cavity output field. Section III presents the main results of this paper, and generalizes our reconstruction procedure to a network of mechanical oscillators. We provide sufficient conditions under which the quantum state of the entire network can be reconstructed, as well as an explicit analytical procedure to design the interaction profile. Section IV deals with the important special case of Gaussian states, and provides a numerical simulation of our scheme for the realistic case of a finite number of measurements. In Section V we present some ideas about extending our state reconstruction scheme to the single-photon regime of optomechanics. Finally, we draw our conclusions in Section VI.

A. Optomechanical model
Our starting point is the so-called linearized optomechanical interaction, which involves a cavity mode with annihilation operator a and an oscillating mirror whose excitations are described via a second annihilation operator b. The cavity is pumped by a resonant classical field, and in a frame rotating with the cavity frequency the Hamiltonian reads where ω m is the mechanical oscillator frequency, g(t) is a time-dependent coupling constant controlled via the amplitude of the classical driving field, and X = (a+a † )/ √ 2 is a cavity quadrature operator. Notice that an alternative scheme to generate Hamiltonian (1) with tunable coupling has been recently reported in [41]. Moving to an interaction picture defined by H 0 = ω m b † b, the Hamiltonian becomes The latter satisfies the Schrödinger equationU = −iH I U , which can be solved via the ansatz U = e iφ D(Xβ) [with D(α) ≡ e αb † −α * b the displacement operator of the mechanical mode], where φ is assumed to be a time-dependent operator that commutes with the mechanical degree of freedom (namely, [φ, b] = [φ, b † ] = 0), whereas β is a time-dependent complex number. Theṅ Using the unitarity of U and the Baker-Campbell-Hausdorff formula, note thaṫ Then the Schrödinger equation implies Matching coefficients produces a set of simultaneous equations: For an interaction time τ , one has the solutions Finally one can rewrite U = e iψX 2 D(Xβ) with Therefore, the dynamics is described by a mechanical displacement operator whose amplitude depends on the in-phase quadrature operator of the optical mode modified by a quadratic term on the optical mode.
For the purposes of this reconstruction strategy, the dynamics can be greatly simplified by constructing the interaction profile in such a way as to eliminate the quadratic term (i.e., setting ψ = 0) [57].

B. State reconstruction of a single oscillator
We will now show that a link can be established between the cavity quadrature operator P = i(a † − a)/ √ 2 (namely, the canonical momentum operator conjugated to X -momentum for brevity) and an arbitrary mechanical quadrature. Then, we will use these results to illustrate how the full reconstruction of the oscillator state can be carried out. Let the initial state of light and mechanics be ρ = |0 0| ⊗ ρ 0 with |0 the vacuum state of the optical mode (in the displaced frame of reference) and ρ 0 the mechanical state to be reconstructed. For the scope of the present discussion it is convenient to switch to the Heisenberg-picture. After an interaction time τ , the cavity field's momentum evolves into where Q θ ≡ (be −iθ + b † e iθ )/ √ 2 is the arbitrarily chosen mechanical quadrature. The phase of Q θ is controlled by the parameter β, as per θ = arg(β) + π 2 . We rescale the observable as from which we can easily deduce a relationship between P r (q), the (measurable) probability distribution of Q r , and P θ (q), the probability distribution of Q θ in the state ρ 0 . Such distributions are related by the convolution integral where the Gaussian factor is due to the second term in Eq. (13) (recall that P exhibits vacuum statistics on our initial state). As the displacement paramenter |β| increases, the measured distribution P r (q) provides a better approximation to P θ (q). In the limit |β| → ∞ one would have P r (q) → P θ (q). Note however that the magnitude of |β| will be limited by physical constraints such as the maximum achievable coupling strength max t |g(t)| and the requirement to keep interaction times short enough to avoid decoherence. For finite |β| it is in principle possible to recover P θ (q) from P r (q) via standard deconvolution techniques, due to the fact that the Gaussian distribution in Eq. (14) is fully known.
As an alternative to the deconvolution approach, one may also exploit Eq. (13) to establish a relation between the statistical moments of the measured operator Q r and the mechanical moments Q k where V k indicates the statistical moments of P in the vacuum state, where Γ is the Euler-Gamma function. Eq. (15) tells us that, having experimentally determined Q j r (j = 1, . . . , n), we may calculate Q j θ (j = 1, . . . , n) from the data by inverting an n × n linear system of equations. As discussed in section IV, this second approach is particularly convenient when ρ 0 is a Gaussian state, in which case the knowledge of first and second moments of Q θ (for several values of θ) is sufficient for full state reconstruction.
Let us now illustrate the quantum state reconstruction protocol. First, the user selects a quadrature Q θ to reconstruct. This choice determines the value of arg(β) as shown above. The modulus of β along with the coupling g(t) and interaction time τ are chosen such that ψ = 0. This can be accomplished by setting and using Eqs. (9, 10) along with the choice of arg(β) and the condition ψ = 0 to solve for the coefficients A and B. After an interaction time τ , the cavity field has assimilated the information from the mechanical mode and is ready to be measured. At this point the coupling is switched off, the measurement is performed, and the result recorded. The system must be reset and the procedure repeated sufficiently many times such that the sampling of the measurement results is reliable. The probability distribution P r (q) and/or the associated moments may then be estimated from the collected data. Subsequently one may proceed to the deconvolution of P r (q), [or the inversion of Eqs. (15)] in order to estimate the distribution P θ (q) (or a finite number of its moments). The procedure must then be repeated for a sufficient number of different quadratures (i.e., different values of θ), such that the state ρ 0 may be recovered via the standard inversion techniques of quantum tomography [4][5][6]42].

C. Measuring the cavity through its output field
Note that the scheme presented so far relies on the measurement of the intra-cavity field, which is typically not directly accessible. However, by making use of the small but inevitable transmittance of the cavity mirrors, one may measure the associated output fields and infer the intra-cavity field properties via input-output theory [43]. For the case under consideration, assume that the emission rate of the cavity κ is small enough as to be negligible during the reconstruction protocol described so far (i.e., κτ 1). It is also convenient to assume that such emission occurs through only one of the two cavity mirrors, so that it may be more easily collected. At time τ , we assume that the optomechanical coupling has been switched off, and that the cavity obeys the standard quantum Langevin equation [43] a whose formal solution can be arranged as This is a Heisenberg-picture relation indicating that the full information about the cavity field (at the time τ of interest) is shared between the output field modes and the cavity field at the later time t f . Observe that for t f κ −1 the desired information is fully encoded in the output field. Formally, this amounts to the expression where the bosonic operator f out represents an appropriate combination of output field modes that can be measured directly. Note that, so far, we have considered an idealized cavity in which all internal losses are associated with emission into detectable modes. In the presence of genuine optical losses, which irreversibly deteriorate the amount of accessible information, Eq. (20) must be modified as follows [44,45].
where 0 < < 1 is the probability of single-photon loss, while a vac is a bosonic mode accounting for the associated added noise: it commutes with f out and is assumed to display vacuum statistics. Constructing the usual quadrature operators we obtain the relation P out = √ 1 − P (τ ) + √ P vac , which can be exploited to obtain relationships analogous to Eqs. (14) and (15), linking the measurement statistics of P out to those of P (τ ). In particular, the measured moments read The moments P n out can be estimated via homodyne detection, so that an inversion of the equations above allows one to retrieve the moments P n (τ ) . As outlined in Section II B, the latter can then be rescaled to obtain the moments of Q r , which in turn allows to retrieve the desired mechanical quadratures. Alternatively, the measured probability distribution of P out may be deconvoluted to obtain that of P (τ ), if the noise parameter is known. Appendix A details examples of the reconstruction when the noise parameter is significant. The primary effect is to increase the number of measurements required for a good representation of P n out . The next significant noise factor is that of damping of the mechanical oscillator at rate Γ. In the regime in which the reconstruction protocol takes place the mechanical damping rate is very small compared to the cavity decay rate. To clarify, we operate in the resolved sideband regime which requires that Γ κ ω m . We thus conclude that the effect of the mechanical damping over the timescale of the interaction with the cavity field is negligible. Furthermore, the explored examples of g(s) (Figs. 4 and 5) required for reconstruction show that the reconstruction is effective over the course of a few mechanical periods.

III. A NETWORK OF OSCILLATORS
A. Hamiltonian and time evolution The protocol described above can be generalized to the case of a network of N harmonically coupled oscillators, whose mechanical excitations are described by a set of bosonic operators Only one mechanical oscillator, say b 1 , is coupled to the optical mode a (see Fig.1). The Hamiltonian for such a system reads H = H 0 + H int , where ω n , J nm and K nm are the bare frequencies and coupling constants characterizing the network -which are assumed to be known in advance [46,47]. By Williamson's theorem, the mechanical portion of our system can be brought into diagonal form by a symplectic transformation S, which has the general structure In terms of the mechanical normal modes defined by S, the Hamiltonian reads where G n = (S 1 − S 2 ) * n1 , d n are the annihilation operators for the normal modes ([d n , d † m ] = δ nm ), and ν n the associated eigenfrequencies. Moving to the interaction picture defined by H 0 , one has where These operators have the property that [h j (t), h j (t )] = 0 ∀j = j . This allows the unitary for the system to be written as U = ⊗ j u j (t) , with each u j satisfying the equatioṅ with the initial condition u j (0) = I . Following Section II, these equations can be solved by the ansatz with φ j commuting with all involved mechanical modes. There are two coupled equations associated with each j which have solutions Then where Ψ = j ψ j , β = β 1 β 2 . . . β N and The evolution of this more general system bears a clear resemblance to the single oscillator case. First, it comprises a quadratic term on the optical mode that depends on the global phase Ψ. We will set this to zero, extending the method used for the reconstruction of one oscillator (see below). Second, we can recognise an X-conditioned multimode displacement on the mechanical modes. In order to proceed with the state reconstruction, it is necessary to introduce two assumptions on the properties of the network. These are • (A1) G n = 0 ∀ n.
These assumptions embody the ability of the cavity field to interact with, and distinguish, all the normal modes of the network [see Eq. (27)]. A further practical requirement is that the interaction time should be sufficiently long to allow the resolution of modes that vibrate with similar frequencies. The dynamics described here is reminiscent of the one derived in Ref. [40], where an oscillator network is probed with an auxiliary two-level system rather than with a cavity field.

B. Quantum state reconstruction
Similarly to the single oscillator case, let the initial state of the system be the factorised state ρ = |0 0| ⊗ ρ 0 where ρ 0 now indicates an arbitrary state of the oscillator network. Reconstruction proceeds as before, with the following modification: there are now multiple mechanical quadratures, defined by Q θj = (d j e −iθj + d † j e iθj )/ √ 2. Note that these quadratures are defined in terms of the normal mode operators, so that the state is reconstructed in the normal-mode basis. However, a representation in terms of the original modes b 1 , ..., b N can be obtained through the inverse symplectic transformation S −1 , corresponding to a reshaping of the reconstructed Wigner function. As before, let us work in the Heisenberg picture, and let us indicate by τ the interaction time in which the controlled displacement is implemented. After the interaction, the cavity quadrature P evolves into where θ j = arg(β j )+ π 2 . By proceeding as for the case of a single mechanical oscillator, one may write a convolution integral connecting the probability distribution of P (τ ) to that of the mechanical quadratures. Since the choice of each quadrature Q θj is determined only by the phase of β j , by varying |β j | in Eq. (36) it is possible to measure a sufficient number of linearly independent observables to enable a deconvolution, and hence estimate the multivariate probability distribution of (Q θ1 , ..., Q θ N ). As we will be primarily interested in Gaussian states, however, here we focus on the reconstruction of arbitrary moments of the mechanical quadratures, rather than their probability distribution. From Eq. (36), it follows that the moments of P are linked to the mechanical quadrature moments by with the sum over all permutations of integers (k 1 , ..., k N ) such that k 1 + k 2 + ... + k N = n, and we recall that V k0 ≡ 0| P k0 |0 is given in Eq. (16). This system of simultaneous equations is under-determined, however we can exploit the freedom in |β j | to generate as many independent extra equations as necessary, involving the same variables but different coefficients. This results in a linear system of equations that, once inverted, provides an arbitrary moment of the mechanical network quadratures. For Gaussian states, we recall that first and second moments will suffice to fully characterize the mechanical quantum state.
As in Section II B, in order to (i) reconstruct an arbitrary quadrature Q θj and (ii) reduce the dynamics to a multimode displacement (i.e. set Ψ = 0)[58], it is crucial to properly select the interaction profile g(s). Let us define the latter as where the additional term outside the sum does not depend on any frequency of the system (i.e. ω is arbitrary). From Eq (32), it follows that where with matrix elements G is the vector of interaction profile coefficients, and where with S being a rectangular matrix. For long interaction times the matrix R is invertible. It can be shown that This implies that there exists an interaction time such that det(R) > 0 and hence R is invertible. Then, as β is chosen, the coefficients G are determined via Eqs. (38) and (48) set the interaction profile as a function of β and h. The additional constraint to be taken into account is that the quadratic parameter Ψ should vanish. As shown in Appendix B, this amounts to a quadratic equation in {h, h * } which can be solved numerically.

IV. GAUSSIAN STATE RECONSTRUCTION
A special case of the reconstruction strategy applies when the state ρ 0 to be reconstructed is Gaussian. Such states are characterised completely by the first and second order moments of two conjugate quadratures (per mode) and the correlations between them. The collection of these terms is directly accessible to the reconstruction scheme outlined here. In other words, if one knows already that the state is Gaussian, one does not have to reconstruct higher order moments.
These relevant cases give the opportunity to show explicitly the functioning of the protocol here introduced and to exemplify its requirements. In particular, we will provide examples of the required time-dependent coupling g(t) in units of the mechanical frequency (or smallest eigenfrequency as appropriate), together with estimates of the protocol's performance in terms of fidelity and of the number of measurements necessary to determine the moments associated with each mechanical quadrature. Figs. 2 and 3 show the behaviour of the fidelity between the reconstructed state and the original state ρ 0 . In order to obtain these plots we proceeded as follows. For the single mode case of Fig. 2, we considered squeezed thermal states with temperature T and squeezing parameter r [48]. We set the total interaction time τ = 2π/ω m and fixed the mechanical quadratures to be reconstructed (Q θ with θ = 0, ±π/4, π/2 suffice in this case). The generic interaction profile is given by Eq. (17). The selection of the mechanical quadrature -together with the additional requirement of deleting the quadratic term e iψX 2 -determine as per Eqs. (9,11) the specific interaction profile g(t). The latter is reported in Fig. 4 for the four cases of interest θ = 0, ±π/4, π/2. Notice that the required tuneability in time is of the order of the mechanical frequency [see also Eq. (17)] and that the profiles are clearly distinguishable, thus indicating robustness against small perturbations. This range of interaction strengths are typically available experimentally. However, due to experimental limitations, it is also possible that, given a fixed β, the required magnitude of g(s) is too high. This obstacle can be circumvented by allowing for a longer interaction time, given that the magnitude of g(s) is inversely proportional to it [see Eq. (38)] For a fixed number of measurements N , we numerically sampled Q m and Q 2 m for the four choices of Q θ . Then, an inversion of Eq. (15) allowed for the reconstruction of the first and second moments of Q θ from which the covariance matrix of the original state can be reconstructed. We then used the latter and the covariance matrix of ρ 0 in order to obtain the fidelity F between the reconstructed state and the original one [49,50]. Fig. 2 shows that, regardless of the state to be reconstructed, a few hundreds of measurements N are sufficient to achieve high fidelity.
For the two-mode case, we considered two-mode squeezed thermal states.
Again, for a fixed number of measurements N and for each choice of {(θ 1 , θ 2 )}, the quadrature P is sampled in order to estimate the second order moments of the mechanical quadratures. In this case, however, extra equations must be generated to make the system in Eq. (37) solvable. Each extra equation costs an additional N measurements. As can be expected, by increasing the number of oscillators to reconstruct the number of required measurements increases. However, the latter is not significantly affected by the state to be reconstructed.
Clearly, tests of non-Gaussianity are also possible within this scheme. In fact, apart from a full reconstruction of the state as explained in Sections II and III, one could check the non-Gaussian character of the state ρ 0 by reconstructing only few higher-order moments and comparing them with the first and second moments. This is a general feature of this scheme, that it allows direct access to partial information of the state without full tomography.

V. SINGLE PHOTON-PHONON COUPLING
In the previous Sections, we have considered a system in which the radiation-pressure interaction has been treated in the linearized regime [see Eqs. (1) and (24)]. This is certainly the situation that has been explored most in experiments to date -both in the opto-and electro-mechanical settings -giving us a clear motivation to focus on it. However, it is worthwhile to briefly outline how our protocol can be modified to the case in which the radiation-pressure coupling retains its nonlinear character. Remarkably, the protocol modifies substantially in this case, and it gives access directly to the characteristic function of the network, as we will now see.
For brevity, we will only mention here the case of one non-linear interacting mechanical resonator b. Eq (1) is which, in an interaction picture defined by the free terms, becomes: The dynamics is solved using the same techniques as above, producing a unitary operator where ψ and β retain their definitions from before and N = a † a. We can see that now dynamics is described by a number-operator conditioned displacement of the mechanical mode modified by a Kerr-like term on the optical mode. Given a factorised initial state, ρ = |α α| ⊗ ρ 0 , where |α denotes a coherent state, the characteristic function of the mechanical state, χ(β) = tr{D(β)ρ 0 }, may be recovered from the first moments of the cavity position and momentum operators. Evolving the cavity's position and momentum operators for a time τ under this displacement operator produces the following relations (53) whereβ = βb † − β * b and X and P are defined as before. In this regime of nonlinear coupling we do not have the tuneable coupling required to set ψ = 0 and therefore the Kerr term cannot be avoided. Including this term, the first moments of X and P are X = α| e −iψN 2 Xe iψN 2 |α tr(coshβ)+ i α| e −iψN 2 P e iψN 2 |α tr(sinhβ) (54) P = α| e −iψN 2 P e iψN 2 |α tr(coshβ)+ i α| e −iψN 2 Xe iψN 2 |α tr(sinhβ) (55) Taking the sum of these produces an expression involving the characteristic function of the mechanical state, Eq. (57) above shows a direct link between the expectation values of X and P and the value of the mechanical characteristic function at point β. By exploring enough points in the phase space it is then in principle possible to reconstruct directly the characteristic function of the mechanical oscillator, which in turn gives full information about its state. This feature is very different with respect to the reconstruction procedure outlined in the previous Sections, and shares a much stricter resemblance with the protocol of [39,40]. In general, the direct reconstruction of the characteristic function entails a set of useful features which have been already outlined in the literature. We refer the reader to Refs. [39,40] (and references therein) for a detailed account.
We would like to emphasise here that the freedom to explore phase space lies in the definition of β [Eq. (9)]. There are two parameters that one could in principle control: the interaction time τ at which the measurement must be performed and the optomechanical coupling strength g 0 . A functional dependence on time for g 0 would grant the greatest control, but this is difficult to achieve experimentally. Its definition is g 0 = ωc L 2mωm , with m the mass of the oscillator and L the length of the cavity. Since these parameters are usually fixed, the only real freedom is in the interaction time. As is clear from Eq. (9), changing the interaction time allows exploration of only a ring in phase space, and not the entire space. However, the partial information on the state from this ring may still provide valuable details on the mechanical state.

VI. CONCLUSIONS
We have introduced a method to reconstruct an arbitrary state of a harmonic network of mechanical oscillators. The reconstruction strategy applies to any setting in which a distinguished mechanical oscillator of the network is coupled to a bosonic probe via a linearized interaction. Then, the main feature of our reconstruction protocol is that by measuring a single system (the probe) the state of the entire mechanical network can be recovered. Given that the probe interacts with one mechanical oscillator only, suitable counter-measures can in principle be envisaged in order to screen the rest of the network from sources of noise that are typically unavoidable whenever a system is coupled to a probe. In this sense, our method provides a minimally invasive configuration to monitor a network of oscillators, contrary to a more standard strategy in which each oscillator of the network is individually measured. This is reminiscent of the approach reported in Ref. [40], where a finite dimensional probe was considered. However, in many settings it is more convenient to use an infinite dimensional probe instead. As said, suitable experimental platforms include opto-and electro-mechanical settings, where linearized coupling at the quantum level has recently been demonstrated between optical or microwave radiation and a single mechanical oscillator [1][2][3]. The main feature that differentiates our setting from the latter is that we consider, rather than only one oscillator, a network of them. However, first implementations of such systems have been reported recently [19][20][21], thus providing a promising route towards the realization of small opto-and electro-mechanical networks. In addition, one can show that our protocol can be adapted to configurations in which the mechanical oscillators, rather then being directly coupled, interact only indirectly via a common cavity mode [51][52][53][54][55][56].
In order to assess the performance of our method, we have considered the relevant case in which the state to be reconstructed is Gaussian. In particular, by giving a detailed analysis of one and two-mode cases, we have shown that the quality of the method is oblivious to the details of the reconstructed state. As one could expect, in order for the method to succeed with high fidelity, the number of required measurements increases with the number of modes. We have explored this feature in some detail by numerically evaluating the fidelity for the realistic case of a finite number of measurements, rather than the limit of an asymptotically large number of measurements as per our analytical results. In addition, we have also shown how the detrimental effect of non-ideal measurements can be taken into account, by considering losses in the coupling between the radiation probe and the mode that is actually measured. Other noise mechanisms could certainly be at work in an actual implementation of our protocol, however they would depend specifically on the platform under consideration and an exhaustive analysis is outside the scope of the present work.
In view of the rapid progress in the development of opto-and electro-mechanical technologies, we believe that the method here introduced could prove useful in assessing the generation of non-classical states of a network of mechanical oscillator, as well as its dynamics. This is of importance in many aspects relevant to the develop-ment of quantum technologies, where the quantumness of a system needs to be assessed in detail while minimally compromising its state and its coherent dynamics. In the plots of Fig. 2 it is assumed that the all optical information leaking from the cavity is measurable, or equivalently, that the intracavity field is accessible to direct measurement. In a physical scenario, the losses described in Eq. (21) must be taken into account. If the loss coefficient is known, then the statistical moments of the intracavity field can be recovered via Eq. (22). However, the effect of the losses is to increase the number of measurements (sampling size) required to accurately represent these moments. Fig. 6 shows the same reconstruction as in Fig. 2 with additional curves demonstrating this effect for various values of .
Appendix B: Eliminating the quadratic term e iΨX 2 To write down explicitly the constraint Ψ = 0, we start by recalling Ψ = Im j τ 0 β jβj * dt . (B1) Using the notation of Section III B, we may perform each integral explicitly Finally, we can exploit Eq. (48) to substitute the explicit expression for G as a linear function of β and h. Thus it is evident that the constraint Ψ = 0 amounts to a quadratic equation in {h, h * }. We have not been able to prove that a solution exists in any instance, however the freedom in the modulus of β j should provide alternatives in case of possible pathological cases.