Out-of-equilibrium thermodynamics of quantum optomechanical systems

We address the out-of-equilibrium thermodynamics of an isolated quantum system consisting of a cavity optomechanical device. We explore the dynamical response of the system when driven out of equilibrium by a sudden quench of the coupling parameter and compute analytically the full distribution of the work generated by the process. We consider linear and quadratic optomechanical coupling, where the cavity field is parametrically coupled to either the position or the square of the position of a mechanical oscillator, respectively. In the former case we find that the average work generated by the quench is zero, whilst the latter leads to a non-zero average value. Through fluctuations theorems we access the most relevant thermodynamical figures of merit, such as the free energy difference and the amount of irreversible work generated. We thus provide a full characterization of the out-of-equilibrium thermodynamics in the quantum regime for nonlinearly coupled bosonic modes. Our study is the first due step towards the construction and full quantum analysis of an optomechanical machine working fully out of equilibrium.


I. INTRODUCTION
As a result of several decades of efforts stemming from different communities, the classical scientific body of thermodynamics have been experiencing a true renaissance.The reasons of this revival can mainly be traced back to the release of two constraints: on the one hand the departure from the thermodynamic limit, motivated by investigation of increasingly smaller systems, enabled fluctuations to be incorporated; on the other hand the tight requirement of quasistatic processes has been relaxed, in favor of generic finite-time transformations connecting non-equilibrium states.The overall picture is an exact, non-perturbative extension of thermodynamics to mesoscopic systems lying arbitrarily far from equilibrium; stochastic thermodynamics [1] is now a mature field which addresses thermodynamical quantities such as work, free energy and entropy at the level of single trajectories and fluctuation theorems relate the value that these quantities assume at equilibrium to outof-equilibrium finite-time dynamics [2,3].
Furthermore, given the ever-increasing control achievable over microscopic systems and the technological quest for devices miniaturization, one would eventually reach a point where quantum fluctuations, besides thermal ones, start playing a non-negligible role [4,5].The former scenario must then be amended with a full quantum treatment.Performances of thermal machines working in the quantum regime have recently been investigated in a plethora of different physical systems [6], and the statistics of relevant figures of merit such as work and entropy generated during time-dependent protocols inquired for different models [7].
Another motivation to achieve a better understanding of thermodynamics in the quantum regime, somehow complementary with respect to the perspective of scaling thermal machines down to the nanoscale, comes from the exploration of macroscopic quantum systems.The extension of quantum-limited control over objects in the mesoscopic-and possibly macroscopic-domain, is of primary interest both for fundamental problems, e.g. the comprehension of the mechanism of decoherence, and for quantum technology.In particular, optomechanical systems provide an ideal platform where to investigate macroscopic quantum phenomena: mechanical oscillators made of 10 15 particles are now approaching the quantum regime, offering unprecedented levels of tunability and control [8].For that reason they are among the most promising candidates to shed light on the interplay between quantum theory and thermodynamics.
In this work we try to merge these scenarios: We explore and characterize the thermodynamical behavior of an optomechanical system driven out of equilibrium by a time-dependent transformation.We address an isolated quantum system, consisting of an optical mode confined in a cavity and parametrically coupled to a mechanical oscillator, evolving according to a time-dependent Hamiltonian and undergoing a two-step measurement protocol.Specifically, we will be concerned with a sudden quench of the interaction, realized by suddenly switching on the coupling between the two, initially uncoupled, modes.We derive analytic expressions for the characteristic function of the work distribution and analyze the full statistics of the work generated.Two different interaction Hamiltonians, both of relevance for present quantum technology, will be considered.We shall first discuss the more common case where radiation-pressure interaction couples the cavity field to the position of the oscillator, followed by the case of a quadratic optomechanical interaction, where the optical field couples to the square of the position of the mechanical resonator.The starting point for most analyses of optomechanical devices is a n m FIG. 1. Graphical depiction of the two-step protocol for the work distribution.At t < 0 a system is in contact with a bath until thermal equilibrium is reached [panel (a)].At t = 0 + , system and bath are detached, while the energy of the system is measured.Let the outcome of such measurement be E 0 n , which projects the state of the system onto the energy eigenstates . The system's Hamiltonian is then changed following to a given protocol and the system evolves according to the unitary evolution operator U (τ, 0) for a time τ [panel (c)], at which time it is measured (over the eigenbasis of the new Hamiltonian).Outcome E τ m is achieved, which gives the new state . By repeating this protocol many times a distribution of values E τ m − E 0 n is achieved, which embodies the probability distribution of the work done by/on the system as a result of the protocol that has been implemented.linearization of the interaction, where the Hamiltonian is cast into a quadratic form that is more amenable to analysis.Here, we eschew this simplification, which is formally valid when the cavity field is strongly driven [9], and address the full nonlinear optomechanical Hamiltonian.We note at this point that the thermodynamical properties of the equivalent linearized model were recently explored by some us in Ref. [10].By retaining the full optomechanical coupling, our work therefore aims to address the out-of-equilibrium thermodynamical behavior of nonlinearly coupled bosonic modes in the quantum regime, and thus go beyond the results reported in literature so far.
The remainder of this work is organized as follows: In Sec.II we introduce the two-measurement protocol necessary to extract the work distribution, and review the quantum fluctuation relations.Sec.III contains a detailed analysis of the dynamical features of an optomechanical system subject to a sudden quench of the coupling parameter and assesses its thermodynamical behavior, first in the case of linear optomechanical coupling and then in the quadratically-coupled case.Finally, in Sec.IV we summarize our findings and discuss new perspectives opened up by this work.

II. WORK DISTRIBUTION AND QUANTUM FLUCTUATION THEOREMS
Let us consider a system described by a timedependent Hamiltonian Ĥ(G t ), whose dependence on time is realized via the externally tunable parameter G t .This parameter, which we refer to as the driving parameter, determines the configuration of the system at any time.Moreover, let us assume that at t = 0 the system is in thermal equilibrium with a bath at inverse temper-ature β, and is hence described by the Gibbs state where Z(G 0 ) = Tr e −β Ĥ(G0) is the canonical partition function of the system.This system is taken out of equilibrium by applying a chosen transformation that modifies G t in time.Here we are concerned with the statistics of the work done on or by the system when applying such a protocol.We thus proceed as follows (cf.Fig. 1 for a graphical depiction of the the process): At time t = 0 + the system is detached from the reservoir and a projective energy measurement is performed on the system in the energy eigenbasis of Ĥ(G 0 ), yielding an eigenstate which we label E 0 n .The driving parameter is changed according to the aforementioned transformation until a final time τ .During this period, the state of the system evolves as dictated by the action of the unitary evolution operator Ûτ,0 on the post-measurement state.Finally, a second projective energy measurement is made on the system, this time in the eigenbasis of Ĥ(G τ ) and yielding eigenstate |E τ m .Given the spectral decompositions of the initial and final Hamiltonians, Ĥ(G , respectively, the energy difference between the two outcomes E τ m − E 0 n may be interpreted as the work performed by the external driving in a single realization of the protocol.This particular value of the work occurs with probability p 0 n p τ m|n , where p 0 n = e −βE 0 n /Z(G 0 ) keeps track of the initial thermal statistics, while p τ m|n = | E τ m | Ûτ,0 E 0 n | 2 embodies the transition probability arising from the change of basis.The work performed due to the protocol described above can be characterized by a stochastic variable W following the probability distribution Instead of dealing directly with Eq. ( 2), it is often useful to work with its Fourier transform χ(u, τ ) = dW e i uW P (W ), which is referred to as the characteristic function of the work distribution and can be cast in the form (3) The utility of the characteristic function becomes apparent when calculating the moments of the work probability distribution explicitly.Indeed, the k th moment of P (W ) can be obtained from the characteristic function as For the special cases of k = 1, 2 it can be shown that this relation acquires the simple form In what follows we will be concerned with a specific driving protocol, known as sudden quench, where G t is abruptly changed from its initial value to the final one.
In this case, Ûτ,0 = 1 and any dependence on τ disappears.We will thus refer to the characteristic function simply as χ(u).
Work fluctuation theorems relate the probability distribution of a given process [cf.Eq. ( 2)] with its timereversed counterpart, and account for the emergence of irreversibility in isolated systems.In the time-reversed (or backward) process the system is initially in a Gibbs state of the final Hamiltonian Ĥ(G τ ), and the transformation acting on the driving parameter is reversed in time as G t → G τ −t .Expressed in terms of the characteristic functions for the forward [χ(u)] and backward [ χ(u)] processes, the Tasaki-Crooks fluctuation relation [12] reads where ∆F = −β −1 log[Z(G τ )/Z(G 0 )] is the free energy difference between the initial states for the forward and backward processes.The main implication of this relation is that the probability to extract an amount of work W from the system during the backward process is exponentially suppressed with respect to the probability that the same amount of work is done on the system during the forward process.Linked to such relation is the celebrated Jarzynski equality [13] which links the average of a quantity arbitrarily far from equilibrium with the state function ∆F .From Eq. ( 7) ∆F ≤ W follows immediately, which embodies a statement of the second principle of thermodynamics.The difference between the two quantities, which we denote by W irr ≡ W − ∆F , is referred to as the irreversible work generated during the transformation.

III. WORK DISTRIBUTION OF QUENCHED OPTOMECHANICAL SYSTEMS
Let us consider the optomechanical interaction between a field mode within a single-mode electromagnetic cavity of resonance frequency ω c and a mechanical resonator characterized by its mass M and oscillation frequency ω m .These two subsystems will be associated to bosonic annihilation operators, denoted by â ([â, â † ] = 1) and b ([ b, b † ] = 1), respectively.The cavity frequency is modulated by, and couples parametrically to, the mechanical displacement x, so that it can be expanded as If the leading term in the expansion is the linear one, the two oscillators interact via radiation-pressure and the much-studied linear optomechanical regime is recovered.On the contrary, if this term vanishes only the position-squared term contributes so that the so-called quadratic optomechanical regime is accessed; examples of physical systems where the latter coupling is achievable are "membrane-in-the-middle" setup [14], levitating nano-beads [15,16], trapped ions or atoms [17].Note that the adjectives 'linear' and 'quadratic' here refer to the power of the mechanical displacement coupled to the field; we stress, however, that the interaction is inherently nonlinear in the field modes, involving three-or four-wave mixing processes.In order to proceed, we assume to be able to control the optomechanical coupling strength, and suddenly turn it on at t = 0 + .As a function of the mechanical position and momentum variables where k = 1 leads to the linear regime and k = 2 to the quadratic one, G is the coupling parameter, and Θ(t) is the Heaviside step function.Since we set G 0 = 0, both systems are initially uncorrelated and prepared in a global thermal state at inverse temperature , and N α = (e β ωα − 1) −1 being the average number of thermal excitations in mode α = c, m.Our main goal is to evaluate the characteristic function of the work distribution Eq. ( 3), which encompasses all the thermodynamically relevant information.Using the above notation, we have Before moving to the calculation of χ(u), P (W ), and ∆F for both linear and quadratic coupling cases, let us make a remark about the implementation of the quench.The somehow contrasting requirements of having an initial equilibrium state of the cavity-mirror system and turning on the optomechanical interaction at a desired time can be reconciled in the following way (here illustrated for the linear coupling case).Let us consider a perfectly reflecting mirror coupled on each side to the field mode âj of cavity c j , j = 1, 2, with equal strength, so that G c1 = −G c2 = G and the interaction Hamiltonian will be given by Ĥint = G (â † 1 â1 − â † 2 â2 )x.If we assume the tripartite system to equilibrate and consider the reduced state of one cavity mode and the mirror we have We can see that, unless the thermal states of the two cavities are perfectly correlated (in a classical way), this state does not reduce to ˆ β , namely the initial state required by the protocol.However, we computed the Kullback-Leibler divergence of the diagonal part ˆ (c1m) (the only entering the protocol) with respect to thermal statistics p k , and we found that in the range of parameters explored in this work it never exceeds values of the order of 10 −4 .Therefore, this configuration may provide a viable method for approximating the initial state of the protocol.The quench would then consist in the sudden shut-off of the auxiliary mode â2 .A detailed feasibility analysis of the whole protocol is however beyond the scope of this work and it is left for future investigations.

A. Quenched linear optomechanical interaction
For the case of a Fabry-Pérot cavity of length L and oscillating mirror of mass M the coupling can be shown to be equal to G (1) t>0 = ω c /L ≡ g/x zpf , where g is referred as the single-photon coupling strength and quantifies the shift in the equilibrium position of the mechanical resonator induced by a single photon.In order to keep the notation as simple as possible, we will explicitly denote by ĤI the (initial) uncoupled Hamiltonian and by ĤF the (final) interacting one It is straightforward to prove that where η = (1 − e −iωmu ) [18].Expression (13) provides us with physical insight into the dynamical evolution induced by radiation-pressure interaction: Apart from two free-rotating terms (the first and last in the above product), the propagator reduces to a displacement of the mechanical mode conditioned on the number of cavity photons, followed by an evolution generated by a Kerrlike term.The characteristic function in Eq. ( 10) can then be explicitly worked out.The form of the interaction suggests taking the trace over the number states {|n c } for mode â and over the coherent states {|α m } for b (we reserve Latin letters for Fock-state labels and Greek letters for coherent-state labels throughout), i.e., where P (m) (α) = exp (−|α| 2 /N m )/(πN m ) is the Glauber-Sudarshan P -representation of an equilibrium thermal state in the coherent-state basis and the compound kets are defined as |n, α ≡ |n c ⊗ |α m .It is possible to gather the following analytical expression for the characteristic function (1 + N c ) n+1 (15) which cannot be summed analytically.We can however appreciate a few significant features of such expression: First, we recognize the thermal statistics of the cavity field modulated by an exponential whose argument keeps track of the average number of phonons N m .Second, the characteristic function is periodic in u.
To proceed further, since the Fourier transform of Eq. ( 15) cannot be explicitly worked out, we evaluate the probability distribution of the work by calculating Eq. ( 2) directly.To do this, we require the energy eigenvalues and eigenstates of ĤI and ĤF .As ĤI is the free Hamiltonian of the uncoupled system, it satisfies the eigenvalue equation ĤI 2 )− g 2 ωm n 2 .A pictorial view of pre-and post-quench eigenstates in the subspace at fixed number n of photons is sketched in Fig. 2. As stated by Eq. ( 2), the transitions from a set of eigenstates to another are responsible-at the microscopic level-for the work performed on or by the system.The probability distribution of the work is thus given by ) where L b a (x) are the generalized Laguerre polynomials coming from the evaluation of the overlap between preand post-quench mechanical oscillator eigenstates [20].A comparison with Eq. ( 2) enables to unambiguously discriminate the contribution of the first projective measurement (which consist of a sampling from the joint thermal distribution of the cavity and the mirror) from the quantum transition probability, and explicitly provides an analytical expression for the latter.The probability distri- bution of the work, together with real and imaginary parts of the characteristic function, is shown in Fig. 3, for different values of N c , N m , and coupling strength.By differentiating the expression of characteristic func-tion Eq. ( 15) and evaluating it in the origin, according to the prescription in Eq. ( 4), one can see that each term of the series identically vanishes, so that the average work generated by quenching the optomechanical coupling is in fact zero.This is in agreement with the behavior of the imaginary part of χ(u), shown in the inset of Fig. 3, which approaches u = 0 with zero derivative; the distribution of the work values is therefore centered around W = 0. Having access to the characteristic function also gives us information about the statistical moments of P (W ); e.g., the variance of the distribution is given by As expected, this quantity increases both with respect to the intensity of the quench, as quantified by g/ω m , and the average number of thermal excitations.This feature is apparent by comparing the topmost distribution, relative to N c = 0.001, N m = 1 and g/ω m = 0.2, to the other two, both obtained for N c = 0.1 and N m = 1-thus varying the ratio ω c /ω m -but corresponding to g/ω m = 0.1 and g/ω m = 0.8 respectively, i.e., increasing both the temperature and the coupling strength.
Let us first analyze P (W ) as illustrated for a few representative cases in Fig. 3, where we consider small values of g/ω m 1.In such conditions and for relatively small values for N c , the probability distribution appears to be dominated by peaks occurring close to multiple values of ω m .These peaks originate from different initiallypopulated Fock states of the mechanical subsystem.Indeed, the number of peaks with appreciable amplitude increases strongly with N m .In Fig. 3 (b) we notice that the sparse peak-distribution associated with very low values of N c changes into a "clustered" one, where groups of peaks develop close to multiples of ω m and are biased towards less positive values of W .This is directly caused by the Kerr-like term in ĤF , whose contribution to the overall energy is always negative.A natural question to ask at this point is why the average work done is zero when each of these fine structures is biased in the same direction.The answer to this lies in the positive skewness of the distribution, which is given by and is more apparent in the low-temperature regime; indeed, by simply looking at the distribution shown in Fig. 3 (b), it is possible to appreciate the positive skewness of the distribution.
Shifting our attention from Fig. 3 to Fig. 4, we can appreciate the effects of increasing the temperature significantly.The two effects we discussed above, namely the increasing number of peaks upon increasing N m and the fine structure that appears more and more prominently when increasing N c , work together to turn P (W ) from a distribution consisting of well-separated peaks to a dense forest of points.It is readily apparent from the latter figure that the tails of the distribution decay exponentially with increasing |W |.In order to investigate this effect more thoroughly, we show in Fig. 4 a coarse-graining of the probability distributions.This coarse-graining was performed by convolving P (W ) with a Gaussian of appropriate width (0.5 ω m in this case).The resulting distributions, drawn as solid curves in this figure, display clearly a tripartite structure.First, around W = 0, a prominent peak is apparent whose width in this figure is entirely due to the convolved Gaussian.Second, a quadratic decay is appreciated for slightly larger values of W .The probability distribution in this region is thus Gaussian in nature.Third, the tails of the distribution have a manifestly exponential character: the coarse-grained curve displays a prominent kink where the exponential tail meets the Gaussian part of the distribution.
It is worth discussing the validity of our coarse-graining approach.We have verified that the discussion above is not modified significantly when the function used to coarse-grain is changed from a Gaussian or a Lorentzian, or when the width of this function is changed within reason.A final check we performed was to construct the cumulative distribution function W −∞ dw P (w).This function was interpolated and smoothed, and then differentiated to give a continuous version of P (W ).Once again, the conclusions we drew above were left unmodified.It is possible to attach a physical meaning to the coarse-graining of P (W ) as follows.Should the probability distribution be measured using any realistic apparatus, the measurement results will not be infinitely sharp, and will be distributed according to some distribution, usually assumed to be Gaussian.Such an experiment would directly yield the coarse-grained distribution we calculate and display in Fig. 4.
We have thus shown, analytically and numerically, that despite turning on a nonlinear interaction between the two modes, on average there is no net production of work.This is perhaps a surprising fact, given that it has been established that either by quenching the frequency of the harmonic potential of a single oscillator [21], or the linear interaction between two bosonic modes [10], net work is produced on average.We shall return to this point in the next subsection, where we discuss the physical origin of this fact and demonstrate a method for producing nonzero average work.
Using Eq. ( 13) we can easily compute the evolution of the initial Gibbs state, as defined by ˆ (t) = e − i t ĤF ˆ (c) β e i t ĤF .In our case, it is easily seen that this always leads to a separable state, where any correlations between the optical and mechanical modes are fully classical.The dynamics is periodic in time: At t = 2πr/ω m (r ∈ Z), the system goes back to the initially factorized state, while for t = (2r + 1)π/ω m (r ∈ Z), one gets the maximally (classically) correlated state.
Eq. ( 13) also allows us to compute the partition function of the system, via a suitable Wick rotation of the argument, i.e., u → −i β, which effectively identifies the imaginary time as an inverse temperature.For the initial state of the system the partition function factorizes in two canonical contributions while for the coupled system we obtain The free energy difference is correspondingly given by e − βωcn e β g 2 n 2 ωm , which, as can be verified, agrees with the Jarzynski equality ∆F = − 1 β ln χ(iβ).Upon close inspection, it is readily apparent that the series involved in the latter expression is actually divergent.Indeed, for every finite value of β, g/ω m , and ω c /ω m , there exists n = n(g, r) such that ∀n > n, we have that g 2 n > r.This causes the sum to diverge exponentially, such that ∆F is formally undefined.This divergent term can be traced back to the part of ĤF that reads ω c â † â − g 2 /ω m (â † â) 2 .As is apparent, the spectrum of this Hamiltonian is not bounded from below.Occupation of levels with n ≥ n, which occurs naturally for any non-zero β, can thus be mapped into a negative temperature with respect to ĤF .To resolve this issue, we impose a cutoff on the number of terms in the series; When g/ω m approaches or even exceeds unity, with the system entering the interesting strong-coupling regime of optomechanics, we must truncate the series to correspondingly small photon numbers in order to prevent dynamical instability, and the ensuing divergence of ∆F , upon quenching the system.For the rest of this work, we will therefore restrict ourselves to the physical domain in which the series does converge.
An explicit calculation of ∆F , as illustrated in Fig. 5, shows that the free energy difference is negative, in agreement with the statement of the second law ∆F ≤ W ≡ 0.Moreover, the irreversible work reduces to W irr = −∆F .Upon moving towards lower temperatures, both the evolved state and the reference thermal state tend to collapse onto the ground state, leading to vanishing values of the irreversible work, as is apparent from the figure.On the other hand, upon increasing the coupling g/ω m , the free energy difference grows in modulus.

B. Initial displacement of the mechanical oscillator
In the previous subsection we observed how W = 0 for an initial thermal state of the Hamiltonian H I , independently of the strength of the quench.The fact can be seen as a direct consequence of the symmetry of the interaction which, being proportional to x, is an odd function in the mechanical field operators, such that In other words, the average work generated by this kind of quench will be zero.In order to remedy this, we now add an initial displacement of amplitude E ω m ∈ R to the mechanical mode b of the Hamiltonian (9) so that the initial and final Hamiltonians will now read ĤI,F,E = ĤI,F + which differs from Eq. ( 15) by a phase factor.This extra factor is actually responsible for positive derivative of the imaginary part Im[χ(u, E)] at the origin and hence to a non-zero value of the average work.Indeed, applying Eq. ( 4), one finds that the average work done by quenching the optomechanical interaction is given by which depends linearly on the displacement E, on the number of thermal photons populating the cavity, and on the quenching parameter.Finally, the free energy difference for this model is given by The behavior of the irreversible work W irr is reported in Fig. 6, with respect to the inverse temperature and the magnitude of the displacement.

C. Quenched quadratic optomechanical interaction
We will consider now the case where the photon number operator of the cavity field is coupled to the square of the position operator of the mirror.As before, we will concentrate on the single-photon regime where the interaction of a single photon with the mechanical mode is enough to appreciably change its frequency and also squeeze its state.In this instance, we can introduce the single-photon coupling strength κ through the relation G (2) = κ/x 2 zpf , in analogy with the linear case.The initial Hamiltonian H I is unmodified and still given by Eq. ( 11), whereas the the post-quench Hamiltonian now reads We choose to work with a non-negative κ, since κ < 0 can introduce post-quench instabilities similar to the one noted for the linear case.The κ > 0 case exhibits no such instabilities.Yet again, we see that this interaction preserves the photon number â † â, so that it proves convenient to write ĤF = ∞ n=0 ĤF,n where each ĤF,n can be cast in the form where Ω n ≡ ω m + 2 κ n and Σ n ≡ 2κ n.Within each such fixed photon-number manifold, we notice the appearance of a modified mechanical frequency, together with a squeezing operator for the mechanical mode whose argument is conditioned on the photon number.The evolution operator relative to the post-quench Hamiltonian can subsequently be expressed as (26) Our next task is to disentangle each exponential operator in the sum.By using the commutation relations between the operators involved in Eq. ( 26), which provide a twoexcitation realization of the su(1, 1) algebra [23], we find e − i ĤF,nu = e where with κ ≡ κ/ω m being a dimensionless quench parameter.We further have the complex quantity ξ n ≡ |ξ n |e iφn whose phase is φ n ≡ η n + π 2 and modulus Armed with this tool we can thus compute the characteristic function of the work distribution, which reads and comes in the form of a thermal average with respect to the cavity distribution-as in Eq. ( 15)-of algebraic functions.Each of the latter is the reciprocal of the square-root of a second degree polynomial in the mean number of phonons N m , whose coefficients are concisely related to each other.Indeed, we can split χ n,0 into its real and imaginary parts, which read Re(χ n,0 ) = cos(ω m u) cos(ω m u √ 1 + 4κ n) (32) We thus have χ n,1 = 2(χ n,0 − 1) and χ n,2 = 2 Re(χ n,0 ) − 1 .As before, since the Fourier transform of Eq. ( 30) cannot be directly evaluated, in order to compute the probability distribution of the work Eq. ( 2) we proceed by diagonalizing the post-quench Hamiltonian ĤF .First, we keep in mind that ĤI is the same as before.However, within any fixed photon number manifold, ĤF,n be diagonalized via a squeezing operation Ŝ(z) = exp(z * b2 /2 − z b † 2 /2) on the mechanical mode conditioned on the photon number n [24].Once again denoting the post-quench quantities with a prime, and expressing the states in the eigenbasis of ĤI , we find eigenstates ĤF,n where the squeezing parameter is given by ζ n ≡ 1 4 log 1+ 4(κ/ω m ) n , and the eigenvalue As sketched in Fig. 7, for the manifold corresponding to n photons, the quench results in a modification of the oscillation frequency which, is multiplied by a factor 1 + 4(κ/ω m ) n , a relative shift of the mechanical levels by ω m 1 + 4(κ/ω m ) n − 1 , and a squeezing of the state by a factor ζ n .Putting everything together, the probability distribution of the work is thus given by where being x the floor function of argument x, which yields the largest integer not greater than x.The probability distribution for the work done on the oscillator in the case of a quadratic interaction, as derived in this section, is illustrated for some representative cases in Figs. 8 and 9.
In order to characterize quantitatively the key features of the distribution of work, here we mention that the average work generated by a quench of the quadratic optomechanical Hamiltonian is different from zero and is then given by hence increasing with respect the occupation numbers of both the cavity and the mechanical mode, as made apparent by inspecting the different panels in Fig. 8.The variance of the distribution reads Finally, the most striking feature of the probability distribution in the case of a quadratic quench is that it is very asymmetrical, fact witnessed by its skewness (38) We note that, for N m 1, it acquires the values 5/ √ 3N c for N c 1 and 74/5 √ 5 for N c 1; both these values are independent of the strength of the quench.As for the linear case the dynamics brings the initial bipartite state of cavity and mechanical mode into a separable sate, given by ˆ entanglement is generated between the two modes.Proceeding in the same manner as before, we can show that the free energy can be cast in the form In this case, too, a suitable Wick-like rotation to imaginary u can be performed to obtain ∆F from χ(u).In practice, however, this calculation is frought with technical difficulties and it is far easier to compute ∆F from an explicit diagonalisation of the Hamiltonian, as was done above.The behavior of the irreversible work for this case has been shown in Fig. (10), and once again we can see how it drops lowering the temperature and increases by increasing the coupling strength.
As in the linear case, is easier to extract a physical meaning behind the various features of these plots by inspecting the respective coarse-grained distributions.First, we see that the positive-W tail still exhibits an approximately exponential decay.It is also apparent that the distribution is, in this case, significantly more skewed towards the right than in the linear case, which can be understood simply through the fact that the post-quench mechanical oscillator frequency is always larger ; even for the case when k = k, therefore, which at least for small κ/ω m has a large probability of occurring, the work done is positive.

IV. CONCLUSIONS AND OUTLOOK
The exploration of out-of-equilibrium features of small systems working in the quantum regime is attracting ever-increasing attention.Optomechanical systems, more so than other systems, offer the tantalizing perspective of naturally bridging the study of quantum thermodynamics with the macroscopic domain.We actually believe that this class of systems offers the possibility of a captivating analogy: Movable mirrors and cavity fields closely resemble pistons and working media in a pistonchamber engine; in turn, this embodies the archetypal example of a thermal machine.In this sense, such systems may serve as the paradigm for understanding a new class of machines, operating both in the quantum regime and far from equilibrium.However, an adequate description of optomechanical systems involves a fully quantum treatment, and a detailed analysis of the thermodynamical properties of them, carried out at a fundamental level and retaining the full nonlinearity of the interaction, has not been conducted thus far.In this work we discussed the generation of work induced by a non-equilibrium transformation in an isolated optomechanical system, quantitatively assessing how an instan-taneous quench of the light-matter coupling affects the thermodynamical response of the system.Our study was grounded through several analytic results, presenting expressions for both the characteristic function of the work distribution and the full statistics of the work generated for two different situations of much relevance for current and future optomechanical experiments.For a quench of linear coupling between light and the position of an oscillator, we found that no work is generated on average, whilst quenching a quadratically-coupled optomechanical interaction requires work to be performed on the system.
Besides being interesting in itself, and allowing for a full analytical treatment, the scenario we addressed comprises the fundamental ingredients necessary in order to gain knowledge about the microscopic origin of the work generated by quenching an optomechanical interaction, from a fully quantum perspective.An in-depth understanding of the thermodynamical response of such an isolated quantum system represents the cornerstone for future investigations.For instance, the implementation of protocols for extracting work out of such systems will require benchmarks based on the analysis that we have performed here, which will in turn be necessary to help uncover fundamental advantages or limitations for possible future thermal machines working in the quantum regime and that exploit the optomechanical interaction.
FIG. 2. Schematic diagram (not to scale) of the energy-level structure of the pre-quench, ĤI,n, and post-quench, ĤF,n, Hamiltonians for the n-photon manifold.Quenching the linear optomechanical interaction results both in an energy shift and a displacement of the machanical oscillator.Two possible transitions induced by the quench-having different values of ∆k = k − k-are shown as an example.

FIG. 3 .
FIG. 3. Logarithmic plot of the probability distribution of the stochastic work variable, W (in units of ωm) for different values of the average number of cavity photons Nc, average number of mechanical phonons Nm and coupling g.Panel (a) is for (Nc, Nm, g) = (0.001, 0.1, 0.2ωm), (b) is for (Nc, Nm, g) = (0.1, 1, 0.1ωm) while (c) for (Nc, Nm, g) = (0.1, 1, 0.8ωm).In the inset is shown the behavior against the time-like variable u (multiplied by ωm) of the real, Re(χ) (solid blue, left), and imaginary, Im(χ) (dashed red, right) parts of the characteristic function.

FIG. 7 .
FIG. 7. Schematic diagram (not to scale) of the energy-level structure of the pre-quench, ĤI,n, and post-quench, ĤF,n, Hamiltonians for the n-photon manifold.Quenching the quadratic optomechanical interaction results both in an energy shift and a squeezing of the frequency of the machanical oscillator.Two possible transitions induced by the quenchhaving different values of ∆k = k − k-are shown as an example.
FIG.9.Logarithmic plot of the probability distribution of work (in units of ωm) corresponding to the parameters (Nc, Nm, κ) = (0.19, 9, 0.7ωm).We also show the coarse grained version of the work distribution (solid magenta line).The coarse graining is realized by convolving the discrete distribution with a Gaussian function of standard deviation 0.9 ωm.