Reconfigurable longrange phonon dynamics in optomechanical arrays
Abstract
We investigate periodic optomechanical arrays as reconfigurable platforms for engineering the coupling between multiple mechanical and electromagnetic modes and for exploring manybody phonon dynamics. Exploiting structural resonances in the coupling between light fields and collective motional modes of the array, we show that tunable effective longrange interactions between mechanical modes can be achieved. This paves the way towards the implementation of controlled phononic walks and heat transfer on denselyconnected graphs as well as the coherent transfer of excitations between distant elements of optomechanical arrays.
pacs:
42.50.Wk, 07.10.Cm, 07.60.Ly, 42.79.GnI Introduction
Optomechanical systems (OMS), naturally lying in the intersection between optical technologies and electronics, play a major role in communication and informationprocessing sciences Roukes (2001). Recent advances in the fabrication of highquality mechanical resonators and their integration with electromagnetic fields have allowed to bring the control of mechanical motion to, or close to, the quantum regime, with potential applications in areas as different as metrology and sensing, quantum information processing, and tests of the fundamental laws of physics Kippenberg and Vahala (2008); Meystre (2013); Aspelmeyer et al. (2013). While these investigations have principally focused on the interplay between electromagnetic radiation and single mechanical resonators, multielement OMS are beginning to be actively studied theoretically Eisert et al. (2004); Bhattacharya et al. (2008); Hartmann and Plenio (2008); Ludwig et al. (2010); Dobrindt and Kippenberg (2010); Stannigel et al. (2010); Heinrich et al. (2011); Chang et al. (2011); Stannigel et al. (2012); Xuereb et al. (2012); Seok et al. (2012); Ludwig and Marquardt (2013); Xuereb et al. (2013); Chesi et al. (2014), as well as experimentally Lin et al. (2010); Mahboob et al. (2011, 2012); Massel et al. (2011); Zhang et al. (2012); Botter et al. (2013). The motivations for exploring their potential are manifold. First, their multimode nature makes them well suited for applications in communication technology Roukes (2001); Stannigel et al. (2010). In addition, they hold the promise for enhanced performance in quantum optomechanics and metrology Xuereb et al. (2012); Seok et al. (2012). Finally, the common interaction of several mechanical elements with one or more electromagnetic fields allows, in principle, for the engineering of complex longrange interactions among the mechanical components, paving the way to the investigation of quantum manybody phenomena with macroscopic elements Ludwig et al. (2010); Xuereb et al. (2012); Ludwig and Marquardt (2013); Tomadin et al. (2012); Xuereb et al. (2013). A key challenge in OMS is to engineer reconfigurable systems, in which the interactions are not predetermined by the bulk properties of the system but can be tailored and switched on or off. This would open the way towards, e.g., efficient and controlled manipulation of heat transfer and single excitations in optomechanical arrays.
In this paper we propose to use periodic optomechanical arrays as reconfigurable platforms for engineering the coupling between multiple mechanical and electromagnetic modes. Such a device operates in a regime where the array is transmissive and light permeates through the structure. This allows both for the enhancement of the optomechanical response Xuereb et al. (2012) and the coupling to specific collective motional modes of the array Xuereb et al. (2013). We show that effective longrange phonon–phonon interactions can be achieved by addressing these transmissive modes. Arising from structural resonances defined by the light fields, these interactions are naturally tunable and reconfigurable. We provide two illustrations of controlled manybody dynamics made possible in this setting: (i) In the badcavity regime of optomechanics, the structure acts as a beamsplitter array for phonons with effective longrange mode coupling, enabling the investigation of phononic random walks on highlyconnected graphs and controlled transfer of heat between distant elements in the array. (ii) In the goodcavity regime, coherent and reconfigurable transfer of single excitations is shown to be possible between distant array elements.
These results should enable the investigation of, e.g., nonstandard heat transport and thermodynamics as well as excitation and information transfer in a wide range of periodically ordered OMS, e.g., nanoelectromechanical resonators Buks and Roukes (2002); Teufel et al. (2011), microtoroids Kippenberg et al. (2005); Arbabi et al. (2011), dielectric membranes Thompson et al. (2008) or particles RomeroIsart et al. (2011), optomechanical crystals Eichenfield et al. (2009), or cold atoms Botter et al. (2013). The engineering of genuine quantum manybody effects in such an array of mesoscopic systems will provide an additional element into the “mechanical quantum simulator” that we propose here. This will allow for addressing, e.g., fundamental issues related to the persistence of quantum features in multielement systems with comparatively large masses, dimensions, and at finite temperature. While these conditions would normally imply Newtonian mechanics, the results presented here suggest that clear signatures of nonclassical behavior can persist even in such a mesoscopic simulator.
Ii Generic model
We consider the dynamics of an externally driven optomechanical system composed of identical mechanical elements, here dubbed ‘mirrors,’ and optical cavity modes. The mirrors are modeled as harmonic oscillators with annihilation operators , vibrational frequency , and decay rate . The ^{th} optical mode is detuned by with respect to its driving field, has a decay rate , and is described by the annihilation operator . Here, and . We treat the mechanical oscillators as a periodic array of lossless 1D scatterers operating in the Lamb–Dicke regime. Such an array displays optical resonances for which the effective reflectivity vanishes Deutsch et al. (1995); Xuereb et al. (2012) and for which the ‘transmissive’ light modes strongly couple to collective motional modes of the array Xuereb et al. (2012, 2013). The Hamiltonian of the system reads (see Appendix) ()
(1) 
where the ‘coupling vectors’ are dimensionless, have unit norm, and are determined mainly by the optical properties of the system. In the case of a periodic array of identical scatterers, at the frequencies where the array is transparent these vectors take the sinusoidal form Xuereb et al. (2013). The optomechanical couplings of the elements thus have a longranged sinusoidal profile spanning the whole array (Fig. 1). Each complex frequency is determined by the mean field amplitude of the respective mode () and the overall optomechanical coupling strength multiplying .
Hamiltonian (1) allows for the engineering of a flexible toolbox for the manipulation of phonon dynamics in an optomechanical array. In the following we shall investigate two regimes: (i) In the ‘badcavity’ regime () we derive an effective Hamiltonian for the mechanics and investigate phonon diffusion and heat transfer through the array. (ii) In the ‘goodcavity’ regime () we derive an analytical expression for the matrix describing the unitary evolution, which allows for the engineering of controlled coherent phonon dynamics.
Iii Badcavity limit
By generalizing the standard procedure Gardiner (1984), the optical fields can be eliminated from the dynamics of the optomechanical system provided that . This yields the effective linearcoupling Hamiltonian (see Appendix) with . The matrices , illustrated in Fig. 2 for , and the tuning of , performed by adjusting and , determine how excitations spread through the array. The resulting system is the phononic analog of the random walks explored in Refs. Broome et al. (2010); Peruzzo et al. (2011); Schreiber et al. (2011); Sansoni et al. (2012); Tillmann et al. (2013); Spring et al. (2013); Crespi et al. (2013); Jeong et al. (2013).
Using the vectors to build an orthonormal ‘similarity’ matrix , we can cast the evolution of the operators describing the mechanical modes as (see Appendix), where . In a similar fashion to Ref. Reck et al. (1994), the matrix can be decomposed into linear optics components (cf. Fig. 3), allowing a general and physically transparent description of the dynamics, and illustrating the way phonons flow through the array.
Figure 4 illustrates a situation where phonons are initially prepared in a coherent state localized at one element of the optomechanical array [panel (a)]. As expected, panel (b) shows that imposes a final population distribution with a sinusoidal shape mimicking that of . [It can be demonstrated numerically that if the coherences between the different modes are set to zero after each step in the interferometer, the resulting “classical” distribution does not bear any resemblance to .] Furthermore, we find that the population distribution can be modified by a proper design of the effective beam splitters and phase shifters described above Jeong et al. (2004). As examples of this flexibility in the manipulation of phonon dynamics, we impose two kinds of randomness on the system: (i) a random phase offset to the phase shifters making up , which can be generated by adding noise to the optical parameters; and (ii) a randomization of the transmission of the beamsplitters in the decomposition of , which corresponds to perturbing , i.e., changing the properties of the mechanical elements Xuereb et al. (2013). In the former case, Fig. 4(c) shows that averaging over many realizations of random phase distributions yields almost uniform phonon populations. In the latter case, panel (d) shows that disorder has instead the opposite effect: the probability distribution collapses into a highlylocalized one with significant population only at the element where the excitation was injected. Combinations of these possibilities can be realized, resulting in a flexible control of the type of phonon walk to be implemented.
Let us now explore the flow of heat throughout the array. The Hamiltonian described above is quadratic and therefore preserves the Gaussian nature of any input state of this kind. With this in mind, we constrain the present analysis to the set of Gaussian states. Each of the mechanical elements is coupled to a Markovian bath characterized by a temperature giving rise to a mean number of excitations in element . We choose for . Therefore, each phonon bath has a mean number of excitations except for that of element , which has . The heat dynamics in the array is then analyzed by solving the differential equation governing the evolution of the covariance matrix of the partite system Ferraro et al. (2005).
The adiabatic elimination discussed previously yields a system of harmonic oscillators coupled not only to the aforementioned thermal baths, but also to shared reservoirs. These reservoirs, which arise through the coupling of each optical field to collective mechanical modes Xuereb et al. (2013), complicate the picture and prevent the standard identification of ‘heat flowing through an element’ in the spirit of Ref. Asadian et al. (2013). The alternative we explore in this section is to calculate the occupation number of the mechanical elements and infer from this the effective heat flow through the array.
The results, illustrated in Fig. 5, exhibit two ‘nonstandard’ behaviors that are intimately tied with the properties of the optomechanical system under consideration. First, this system exhibits static reconfigurability, i.e., the form of the steadystate phonon population distribution can be chosen by changing which of the optical fields is used to induce the indirect couplings between them. Whilst it is not possible to choose an arbitrary distribution, owing to the symmetry of the effective Hamiltonian, this choice still admits considerable flexibility. Second, the flow of energy from one mechanical element to another takes place indirectly, through the cavity field. It therefore proceeds at a similar rate throughout the entire array, governed not by the distance between the source element and the element in question but by the coupling constant of the latter to the optical field. A corollary of this is: if for some , one can speak of heat flow from the element to another element without necessitating any form of heat conduction through element itself. This situation occurs, e.g., for and whenever is odd Xuereb et al. (2013). For even , the elements closest to the center of the array are the least affected. What distinguishes opticallymediated from direct coupling is thus (i) reconfigurability; and (ii) timescales, as excitations flow to every element simultaneously in the optical case, rather than sequentially. These studied interactions enable the study of heat transfer and thermodynamics in nonstandard settings Dubi and Di
Ventra (2011); Bermudez et al. (2013). The parameters selected for plotting Fig. 5 were such that for the steadystate occupation numbers were, to a good approximation, all equal to , regardless of the cooling effects of the optomechanical coupling.
Iv Goodcavity limit
We now turn to the ‘goodcavity’ regime, for which . This allows us to neglect the nonunitary dynamics in an approximate picture if we confine ourselves to times . Upon setting , moving into a rotating frame with respect to the free Hamiltonian (i.e., with ), and neglecting rapidlyoscillating terms in the Hamiltonian, we obtain the evolution operator
(2) 
Our interest lies in the coherent shuttling of single excitations around the system. We therefore constrain ourselves to the singleexcitation subspace and express the state vector as a dimensional vector with the first (last ) entries representing the probability amplitude of the excitation to be found in the respective optical (mechanical) mode. To simplify the notation, let us define the matrix ( columns, rows). It can then be shown that the unitary evolution matrix can be written in the blockmatrix form (see Appendix)
(3) 
where , , and ^{1}^{1}1This is a multimode generalization, restricted to the singleexcitation subspace, of the Jaynes–Cummings evolution operator Stenholm (1973)..
In principle, this evolution can even be dynamically reconfigurable if we allow for the possibility that the amplitudes of the optical modes can be changed on a timescale , and therefore significantly shorter than any other timescale of the problem. The implementation of this is discussed in detail in the Appendix; we note that it is crucial that this switching occurs when the mechanical and optical subsystems are uncorrelated and no excitations reside in the optical subsystem. With this in mind, we can therefore string together sequences of , between which the amplitudes are changed instantaneously.
The result of this procedure is a set of linear equations that allow us to engineer the route of an excitation through the array. As an example, we illustrate the case for , where the fact that independently of and allows for particularly simple protocols to be devised. We demonstrate our ideas by means of the two different examples shown in Fig. 6: (a) By switching the amplitudes of two fields, we transport a phonon from mirror to mirror , and (b) starting from an initial superposition of the phonon on mirrors and , we end up with a polariton oscillating between mirrors and and the light fields.
V Discussion and outlook
We have investigated collective dynamics in multimode OMS with the goal of simulating manybody effects. The dynamical regimes considered in our analysis showcases distinctive possibilities, ranging from diffusionlike propagation of phononic excitations across the array to the controlled transfer of phonons between targeted elements of the mechanical system. Other regimes of interest could be similarly explored. For instance, operating with bluedetuned cavity fields would allow for investigating collective selfoscillations and synchronization Ludwig and Marquardt (2013) in such systems; exploiting the intrinsic nonlinearity of the optomechanical coupling could enable simulation of manybody models (e.g., the Bose–Hubbard Hamiltonian) Jacobs (2012) or quantum information processing Rips and Hartmann (2013) with mechanical systems; and using ring cavities would allow exploring geometric phases Chesi et al. (2014). Such studies are promising for engineering nontrivial manybody dynamics, a possibility we plan to pursue in future works addressing dissipative quantum state engineering, dynamical phase transitions, and fluctuation theorems of thermodynamics origin Mazzola et al. (2013); Dorner et al. (2013); Mazzola et al. (2014).
Vi Acknowledgements
A.X. would like to thank C. Di Franco for interesting discussions and the Royal Commission for the Exhibition of 1851 for financial support. C.G. acknowledges support from the Austrian Science Fund (FWF): P24968N27, and G.P. from the European Commission through ERC StGrant “COLDSIM” (No. 307688), AFOSR, and UdS Labex. M.P. thanks the UK EPSRC for a Career Acceleration Fellowship and a grant awarded under the “New Directions for Research Leaders” initiative (EP/G004579/1), the Alexander von Humboldt Stiftung, the John Templeton Foundation (grant 43467), and the EC Collaborative project TherMiQ (No. 618074). A.D. acknowledges funding from the EU (CCQED project), the Institut Francais du Danemark (IFD2013 program) and the Danish Council for Independent Research (Sapere Aude program). Some calculations were carried out using computational facilities funded by the European Regional Development Fund, Project ERDF080.
Appendix A Choice of working point: Linearization
Throughout the main text, we worked exclusively in the linearized domain, assuming that the singlephoton coupling strength between each of the mechanical elements and the optical fields is small. As is wellknown Meystre (2013); Aspelmeyer et al. (2013), upon application of a strong driving field this means that the effective optomechanical coupling strength is equal to the singlephoton coupling strength multiplied by the squareroot of the expectation value of the photon number inside the relevant optical mode. This photon number can be chosen arbitrarily through the strength of the driving field and, therefore, so can the effective coupling strength. This means that we do not need to optimize the parameters of the system (in the optical case, this would mean having to carefully choose highlyreflective mirrors and a small cavity) as done when trying to obtain strong singlephoton coupling: What we require is simply that the physical configuration is chosen such that the system is in the transmissive configuration defined in Ref. Xuereb et al. (2012). We note, however, that experimental limitations (including absorption and limited laser power) and optomechanical instabilities act to limit the maximum effective coupling strength achievable in any one setup. Under such conditions, it may still be advantageous to design the system to optimise the singlephoton coupling strength.
This choice of working point is crucial in distinguishing the family of systems we explore from that proposed in Ref. Seok et al. (2012) and others in the literature. Indeed, our protocols require:

at least two optical modes interacting with different mechanical oscillators, and

that the resulting collective mechanical oscillators being different superpositions of the same physical oscillators.
Appendix B Adiabatic elimination of the optical field in the optomechanical master equation
Here we briefly go through the steps for adiabatically eliminating the cavity field from the master equation describing an optomechanical system consisting of one optical field (annihilation operator ) interacting with one mechanical mode (). Our treatment follows that in Ref. Jaehne et al. (2008), the Supplementary Information of Ref. Stannigel et al. (2012), and §5.1.2 in Ref. Gardiner and Zoller (2004). We begin by splitting the optomechanical Hamiltonian Law (1995) that describes linear coupling to the optical field into the free () and interaction () parts:
(4) 
together with the cavitydriving Hamiltonian , where the (coherent) driving field amplitude is taken to be real for simplicity. To account for nonunitary dissipative processes, we consider the density matrix for the composite system and introduce the superoperators
(5) 
and
(6) 
where represents the mean number of excitations in the bath that is coupled to. Given the decay rates and for and , respectively, we can now write the master equation of our system as
(7) 
Our first step is to define two unitary transformations and , which act on and , respectively, via the relations
(8) 
These are used to transform our master equation to (i) eliminate , and (ii) shift and to have zero mean value. The resulting Hamiltonian can be truncated to be quadratic in the operators, and can be written
(9) 
after a redefinition of to absorb . Let us call these three terms , , and , respectively. At this point we define the system (‘s’), bath (‘b’), and interaction (‘i’) Liouvillians as follows:
(10)  
(11) 
and
(12) 
and our master equation reads . We are interested in an equation of motion for the density matrix with eliminated from it. Our treatment will only be valid to lowest order in . We define an operator such that
(13) 
which projects onto a tensor product of the vacuum state for the cavity (), which is the steadystate solution for the bath master equation, and the reduced density matrix where the cavity field has been traced out. It is easy to see that is a projection operator () and that, if we define as the identity operator, so is . Two properties of these operators are particularly useful to us (Gardiner and Zoller, 2004, §5.2.1):
(14)  
(15) 
By projecting our master equation, we therefore obtain
(16)  
(17) 
The adiabatic (‘weak coupling’, in the sense of being lowest order in ) approximation allows us to formally solve the second equation:
(18) 
We now take the initial time (which allows us to ignore the first term) and substitute, in the integrand, the zeroorder approximation
(19) 
The in the exponent can be ignored to this order of approximation in :
(20) 
or, since commutes with ,
(21)  
(22) 
Finally, we trace over the cavity field to obtain the approximate master equation for the reduced density matrix ,
(23) 
Substitution of into this expression and application of the rotatingwave approximation, after some algebra, give the final result
(24) 
where the spectral density of the cavity field to zeroth order in is defined by
(25) 
We can therefore identify an effective Hamiltonian
(26) 
as well as the effective cooling and heating Liouvillians
(27) 
and
(28) 
This calculation is readily generalized to many fields at different frequencies, which do not interfere, and multiple mechanical elements.
Appendix C Basis vectors diagonalizing the Hamiltonian
Let us start from the effective Hamiltonian
(29) 
As shown in previous work Xuereb et al. (2013), we note that the vectors and are identical for symmetry reasons. Therefore, the set of vectors does not span the dimensional space necessary to describe the mechanical oscillators. Indeed, in general this set has linearlyindependent vectors. For this reason, we can absorb the into the , and set . To proceed, let us complete basis by adding necessary unit vectors, forming the set of orthonormal basis vectors , where for . The choice of these vectors is completely arbitrary so long as forms an orthonormal basis for the space. For the most part, different choices of the vectors with lead to the same dynamics. However, this does not hold true when we add randomness to the system (as we did to generate Fig. 4 of the main text).
Let us now define the matrix whose rows are the , as well as the vector of operators and the diagonal matrix . This leads us to a simple diagonal decomposition of the Hamiltonian: . In the Heisenberg picture, we therefore simply have .
Appendix D Choice of parameters for Fig. 5
Standard optomechanical theory applies to the system in Fig. 5 in the main text. This means that the light field can itself create or destroy collective mechanical excitations, even if all the mechanical baths are at the same temperature. To present the clearest case, we chose the parameters for this figure to nullify this optical effect: If the baths are held at the same temperature, the resulting populations would then be the same for all the oscillators. For clarity, only one optical mode was activated (i.e., each was set to zero except for one). We chose to use the mode because it has the longest spatial period and is therefore clearest to depict in figures, but any other mode would have done; we chose , but any other value would have given rise to qualitatively similar data. Next, was selected such that the resulting quality factor of the mechanical oscillator, , well within experimental reach Karuza et al. (2013). Following this, was chosen to lie on the red sideband of the oscillator. Finally, was selected to minimize the opticalonly effect.
Appendix E Evolution operator matrix
Suppose we have mechanical elements and optical fields. Working in the singlephonon subspace, we will work with kets of the form and , where and . The former indicates a photon in field , the latter a phonon on element . Our free Hamiltonian is
(30) 
and the interaction Hamiltonian
(31) 
is the detuning of the ^{th} optical field, the mechanical frequency of the ^{th} mechanical oscillator, the optomechanical coupling frequency between the ^{th} field and ^{th} oscillator, and the cnumber component of the ^{th} field. In the interaction picture with respect to , and taking we obtain the Hamiltonian
(32)  
(33) 
We now define the evolution operator
(34)  
(35) 
noting that this expression may be generalized appropriately for . To aid analysis, define
(36)  
(37) 
such that . Thus, for example,
(38)  
(39)  
(40)  
(41) 
etc. Define the matrix ( columns, rows):
(42) 
which allows us to compactly write quantities of the type
(43) 
Now, with each we associate a element vector
(44) 
With this notation, and generalizing the above work, we obtain as in the main text, where we also made use of the definitions
(45)  
(46) 
When is singular, the form of the resulting equation allows us to use the Moore–Penrose pseudoinverse to compute .
Appendix F Nonunitary evolution in the goodcavity limit
When considering the singleexcitation subspace in the goodcavity limit we ignored decay, both optical and mechanical, and took the rotatingwave approximation (RWA). In this section, we present some data to justify these approximations. First, in Fig. 7, we simulated the master equation using the full linearized Hamiltonian and the standard dissipation channels. For the chosen parameters (, , , ), where both the optical and mechanical decay rates are very small, no decay is evident and the unitary approximation, together with the RWA, can be freely used. Breakdown of the RWA could be detected through appearance of higherfrequency oscillations in the populations, which are absent in the presented data. We are presently exploring the positivedetuning case, where the Hamiltonian approximates the “twomode squeezingââ Hamiltonian Meystre (2013); Aspelmeyer et al. (2013). Under these conditions, one can no longer restrict the analysis to the singleexcitation subspace and a concise analysis of the type we explore in the main text is not available.
Appendix G Dynamical reconfigurability
In the main text we make use of the fact that the system is dynamically reconfigurable. Our protocols necessitate the ability to change the amplitudes of the optical modes on a timescale , and therefore significantly shorter than any other timescale of the problem. This can be achieved by a rapid change of the cavity linewidth Xu et al. (2007); Kondo et al. (2013); Liu et al. (2013).
Furthermore, a straightforward calculation shows that if is a function of time, the input optical field amplitude can be changed simultaneously to yield an intracavity field amplitude that is constant in time. It is simplest to illustrate this in the case when changes instantaneously in time, and for a single optical mode; we therefore drop the index labelling the optical modes. Suppose that there exists some time for which and . The equation of motion for the mean amplitude of the cavity field can be written
(47) 
where includes corrections due to the mean amplitude of the mechanical motion. The solution for this equation for a time is
(48)  
(49) 
where the first line refers to the case where and the second . Let us now suppose that for and for . Then,
(50)  
(51) 
We now choose such that
(52) 
in which case does not depend on time.
With this in mind, and considering that for the optical mode , one can switch quasiinstantaneously by using the following protocol:

Start from the steadystate at time , where the intracavity mean amplitude is timeindepenent

Increase the cavity linewidth, quasiinstantaneously, simultaneously adjusting the input field strength according to the above prescription such that the intracavity field amplitude stays constant

Change the input amplitude; the cavity field readjusts over a timescale

Switch the cavity linewidth back to its original value , once again adjusting the input field amplitude to keep the intracavity field amplitude constant
The above protocol ‘opens up and closes’ the cavity, such that any excitations in the optical field are lost. Therefore, we require the excitation to be found entirely in the mechanical subspace during this manipulation, such that the state of the system is not affected by this procedure.
The net effect of this protocol is that for the dynamics proceeds as before, albeit with a stepchange in the value of . With this in mind, we can therefore string together sequences of , between which (at times ) the amplitudes are changed instantaneously, to form the evolution matrix
(53) 
We note that there are two constraints one must keep in mind when implementing this protocol. First, as already discussed, the switching must occur at a time when the mechanical and optical subsystems are uncorrelated, and when the excitation resides entirely in the mechanical subsystem. Second, the switching time must be much slower than the timescale set by the inverse of the optical frequency in order to avoid creating photons through mechanisms such as the dynamic Casimir effect.
Appendix H Practical implementations
The generic features described in the main text are in principle applicable to a wide range of systems, whether in the optical or the microwave domain. While we used a language pertaining mostly to the optical domain language in the main text, micromechanical elements coupled to superconducting microwave resonator fields represent a very promising system for observing these effects. Indeed, resolved sideband cooling to the ground state and operation in the strong coupling regime Teufel et al. (2011) have been achieved, and even scalable experiments with multiple elements have been recently carried out Massel et al. (2012). Initialization and readout of the mechanics could be performed either by coupling the mechanical elements to additional microwave resonators or to artifical atoms, as has also been demonstrated in Refs. O’Connell et al. (2010); Pirkkalainen et al. (2013), amongst others. In the optical domain potential candidates could be dielectric membranes Thompson et al. (2008); Karuza et al. (2013); Purdy et al. (2013), forming a periodic array inside an optical resonator, as studied in e.g., Refs. Xuereb et al. (2012); Seok et al. (2012). Using parameters close to the experiments of Ref. Karuza et al. (2013), for instance, one could operate in the good cavity regime discussed in the paper. Initialization and readout could be performed either optically using standard techniques or electrically using functionalized membranes Lee et al. (2013); Adiga et al. (2013); Bagci et al. (2013); Andrews et al. (2013). In the same fashion toroidal cavities with indentations Arbabi et al. (2011) could also be used. Another promising candidate for the realization of such a mechanical simulator would be ensembles of cold atoms in optical cavities StamperKurn (2012); Brahms et al. (2012), for which the required regimes and techniques have all been experimentally demonstrated.
References
 Roukes (2001) M. Roukes, Phys. World 14, 25 (2001).
 Kippenberg and Vahala (2008) T. J. Kippenberg and K. J. Vahala, Science 321, 1172 (2008).
 Meystre (2013) P. Meystre, Ann. Phys. 525, 215 (2013).
 Aspelmeyer et al. (2013) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, eprint (2013), arXiv:1303.0733 .
 Eisert et al. (2004) J. Eisert, M. B. Plenio, S. Bose, and J. Hartley, Phys. Rev. Lett. 93, 190402 (2004).
 Bhattacharya et al. (2008) M. Bhattacharya, H. Uys, and P. Meystre, Phys. Rev. A 77, 033819 (2008).
 Hartmann and Plenio (2008) M. J. Hartmann and M. B. Plenio, Phys. Rev. Lett. 101, 200503 (2008).
 Ludwig et al. (2010) M. Ludwig, K. Hammerer, and F. Marquardt, Phys. Rev. A 82, 012333 (2010).
 Dobrindt and Kippenberg (2010) J. M. Dobrindt and T. J. Kippenberg, Phys. Rev. Lett. 104, 033901 (2010).
 Stannigel et al. (2010) K. Stannigel, P. Rabl, A. S. Sørensen, P. Zoller, and M. D. Lukin, Phys. Rev. Lett. 105, 220501 (2010).
 Heinrich et al. (2011) G. Heinrich, M. Ludwig, J. Qian, B. Kubala, and F. Marquardt, Phys. Rev. Lett. 107, 043603 (2011).
 Chang et al. (2011) D. E. Chang, A. H. SafaviNaeini, M. Hafezi, and O. Painter, New J. Phys. 13, 023003 (2011).
 Stannigel et al. (2012) K. Stannigel, P. Komar, S. J. M. Habraken, S. D. Bennett, M. D. Lukin, P. Zoller, and P. Rabl, Phys. Rev. Lett. 109, 013603 (2012).
 Xuereb et al. (2012) A. Xuereb, C. Genes, and A. Dantan, Phys. Rev. Lett. 109, 223601 (2012).
 Seok et al. (2012) H. Seok, L. F. Buchmann, S. Singh, and P. Meystre, Phys. Rev. A 86, 063829 (2012).
 Ludwig and Marquardt (2013) M. Ludwig and F. Marquardt, Phys. Rev. Lett. 111, 073603 (2013).
 Xuereb et al. (2013) A. Xuereb, C. Genes, and A. Dantan, Phys. Rev. A 88, 053803 (2013).
 Chesi et al. (2014) S. Chesi, Y.D. Wang, and J. Twamley, eprint (2014), arXiv:1402.0926 .
 Lin et al. (2010) Q. Lin, J. Rosenberg, D. Chang, R. Camacho, M. Eichenfield, K. J. Vahala, and O. Painter, Nat. Photonics 4, 236 (2010).
 Mahboob et al. (2011) I. Mahboob, E. Flurin, K. Nishiguchi, A. Fujiwara, and H. Yamaguchi, Nat. Commun. 2, 198 (2011).
 Mahboob et al. (2012) I. Mahboob, K. Nishiguchi, H. Okamoto, and H. Yamaguchi, Nat. Phys. 8, 387 (2012).
 Massel et al. (2011) F. Massel, T. T. Heikkilä, J.M. Pirkkalainen, S. U. Cho, H. Saloniemi, P. J. Hakonen, and M. A. Sillanpää, Nature 480, 351 (2011).
 Zhang et al. (2012) K. Zhang, P. Meystre, and W. Zhang, Phys. Rev. Lett. 108, 240405 (2012).
 Botter et al. (2013) T. Botter, D. W. C. Brooks, S. Schreppler, N. Brahms, and D. M. StamperKurn, Phys. Rev. Lett. 110, 153001 (2013).
 Tomadin et al. (2012) A. Tomadin, S. Diehl, M. D. Lukin, P. Rabl, and P. Zoller, Phys. Rev. A 86, 033821 (2012).
 Buks and Roukes (2002) E. Buks and M. L. Roukes, J. Microelectromech. S. 11, 802 (2002).
 Teufel et al. (2011) J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds, Nature 475, 359 (2011).
 Kippenberg et al. (2005) T. J. Kippenberg, H. Rokhsari, T. Carmon, A. Scherer, and K. J. Vahala, Phys. Rev. Lett. 95, 033901 (2005).
 Arbabi et al. (2011) A. Arbabi, Y. M. Kang, C.Y. Lu, E. Chow, and L. L. Goddard, Appl. Phys. Lett. 99, 091105 (2011).
 Thompson et al. (2008) J. D. Thompson, B. M. Zwickl, A. M. Jayich, F. Marquardt, S. M. Girvin, and J. G. E. Harris, Nature 452, 72 (2008).
 RomeroIsart et al. (2011) O. RomeroIsart, A. C. Pflanzer, M. L. Juan, R. Quidant, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, Phys. Rev. A 83, 013803 (2011).
 Eichenfield et al. (2009) M. Eichenfield, J. Chan, R. M. Camacho, K. J. Vahala, and O. Painter, Nature 462, 78 (2009).
 Deutsch et al. (1995) I. H. Deutsch, R. J. C. Spreeuw, S. L. Rolston, and W. D. Phillips, Phys. Rev. A 52, 1394 (1995).
 Gardiner (1984) C. W. Gardiner, Phys. Rev. A 29, 2814 (1984).
 Broome et al. (2010) M. A. Broome, A. Fedrizzi, B. P. Lanyon, I. Kassal, A. AspuruGuzik, and A. G. White, Phys. Rev. Lett. 104, 153602 (2010).
 Peruzzo et al. (2011) A. Peruzzo, A. Laing, A. Politi, T. Rudolph, and J. L. O’Brien, Nat. Commun. 2, 224 (2011).
 Schreiber et al. (2011) A. Schreiber, K. N. Cassemiro, V. Potoček, A. Gábris, I. Jex, and C. Silberhorn, Phys. Rev. Lett. 106, 180403 (2011).
 Sansoni et al. (2012) L. Sansoni, F. Sciarrino, G. Vallone, P. Mataloni, A. Crespi, R. Ramponi, and R. Osellame, Phys. Rev. Lett. 108, 010502 (2012).
 Tillmann et al. (2013) M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit, and P. Walther, Nat. Photonics 7, 540 (2013).
 Spring et al. (2013) J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Kolthammer, X.M. Jin, M. Barbieri, A. Datta, N. ThomasPeter, N. K. Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R. Smith, and I. A. Walmsley, Science 339, 798 (2013).
 Crespi et al. (2013) A. Crespi, R. Osellame, R. Ramponi, V. Giovannetti, R. Fazio, L. Sansoni, F. De Nicola, F. Sciarrino, and P. Mataloni, Nat. Photonics 7, 322 (2013).
 Jeong et al. (2013) Y.C. Jeong, C. Di Franco, H.T. Lim, M. S. Kim, and Y.H. Kim, Nat. Commun. 4, 2471 (2013).
 Reck et al. (1994) M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, Phys. Rev. Lett. 73, 58 (1994).
 Jeong et al. (2004) H. Jeong, M. Paternostro, and M. S. Kim, Phys. Rev. A 69, 012310 (2004).
 Ferraro et al. (2005) A. Ferraro, S. Olivares, and M. G. A. Paris, Gaussian states in quantum information, Napoli series on Physics and Astrophysics (Bibliopolis, Napoli, 2005).
 Asadian et al. (2013) A. Asadian, D. Manzano, M. Tiersch, and H.J. Briegel, Phys. Rev. E 87, 012109 (2013).
 Dubi and Di Ventra (2011) Y. Dubi and M. Di Ventra, Rev. Mod. Phys. 83, 131 (2011).
 Bermudez et al. (2013) A. Bermudez, M. Bruderer, and M. B. Plenio, Phys. Rev. Lett. 111, 040601 (2013).
 (49) This is a multimode generalization, restricted to the singleexcitation subspace, of the Jaynes–Cummings evolution operator Stenholm (1973).
 Stenholm (1973) S. Stenholm, Phys. Rep. 6, 1 (1973).
 Jacobs (2012) K. Jacobs, eprint (2012), arXiv:1209.2499 .
 Rips and Hartmann (2013) S. Rips and M. J. Hartmann, Phys. Rev. Lett. 110, 120503 (2013).
 Mazzola et al. (2013) L. Mazzola, G. De Chiara, and M. Paternostro, Phys. Rev. Lett. 110, 230602 (2013).
 Dorner et al. (2013) R. Dorner, S. R. Clark, L. Heaney, R. Fazio, J. Goold, and V. Vedral, Phys. Rev. Lett. 110, 230601 (2013).
 Mazzola et al. (2014) L. Mazzola, G. De Chiara, and M. Paternostro, eprint (2014), arXiv:1401.0566 .
 Jaehne et al. (2008) K. Jaehne, K. Hammerer, and M. Wallquist, New J. Phys. 10, 095019 (2008).
 Gardiner and Zoller (2004) C. W. Gardiner and P. Zoller, Quantum Noise, 3rd ed. (Springer, 2004).
 Law (1995) C. Law, Phys. Rev. A 51, 2537 (1995).
 Karuza et al. (2013) M. Karuza, C. Biancofiore, M. Bawaj, C. Molinelli, M. Galassi, R. Natali, P. Tombesi, G. Di Giuseppe, and D. Vitali, Phys. Rev. A 88, 013804 (2013).
 Xu et al. (2007) Q. Xu, P. Dong, and M. Lipson, Nat. Phys. 3, 406 (2007).
 Kondo et al. (2013) K. Kondo, M. Shinkawa, Y. Hamachi, Y. Saito, Y. Arita, and T. Baba, Phys. Rev. Lett. 110, 053902 (2013).
 Liu et al. (2013) Y.C. Liu, Y.F. Xiao, X. Luan, and C. W. Wong, Phys. Rev. Lett. 110, 153606 (2013).
 Massel et al. (2012) F. Massel, S. U. Cho, J.M. Pirkkalainen, P. J. Hakonen, T. T. Heikkilä, and M. A. Sillanpää, Nat. Commun. 3, 987 (2012).
 O’Connell et al. (2010) A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang, M. Weides, J. Wenner, J. M. Martinis, and A. N. Cleland, Nature 464, 697 (2010).
 Pirkkalainen et al. (2013) J.M. Pirkkalainen, S. U. Cho, J. Li, G. S. Paraoanu, P. J. Hakonen, and M. A. Sillanpää, Nature 494, 211 (2013).
 Purdy et al. (2013) T. P. Purdy, R. W. Peterson, and C. A. Regal, Science 339, 801 (2013).
 Lee et al. (2013) S. Lee, V. P. Adiga, R. A. Barton, A. M. van der Zande, G.H. Lee, B. R. Ilic, A. Gondarenko, J. M. Parpia, H. G. Craighead, and J. Hone, Nano Lett. 13, 4275 (2013).
 Adiga et al. (2013) V. P. Adiga, R. De Alba, I. R. Storch, P. A. Yu, B. Ilic, R. A. Barton, S. Lee, J. Hone, P. L. McEuen, J. M. Parpia, and H. G. Craighead, Appl. Phys. Lett. 103, 143103 (2013).
 Bagci et al. (2013) T. Bagci, A. Simonsen, S. Schmid, L. G. Villanueva, E. Zeuthen, J. Appel, J. M. Taylor, A. Sørensen, K. Usami, A. Schliesser, and E. S. Polzik, eprint (2013), arXiv:1307.3467 .
 Andrews et al. (2013) R. W. Andrews, R. W. Peterson, T. P. Purdy, K. Cicak, R. W. Simmonds, C. A. Regal, and K. W. Lehnert, eprint (2013), arXiv:1310.5276 .
 StamperKurn (2012) D. M. StamperKurn, eprint (2012), arXiv:1204.4351 .
 Brahms et al. (2012) N. Brahms, T. Botter, S. Schreppler, D. W. C. Brooks, and D. M. StamperKurn, Phys. Rev. Lett. 108, 133601 (2012).