Out-of-equilibrium thermodynamics of quantum optomechanical systems

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_{t\gt 0}^{(1)}={{\omega }_{c}}/L\equiv g/{{x}_{{\rm 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 ${{\hat{H}}_{I}}$ the (initial) uncoupled Hamiltonian

$$\begin{eqnarray}&&{{\hat{H}}_{t=0}}=\hbar {{\omega }_{c}}{{\hat{a}}^{\dagger }}\hat{a}+\hbar {{\omega }_{m}}({{\hat{b}}^{\dagger }}\hat{b}+\frac{1}{2})\equiv {{\hat{H}}_{I}},\end{eqnarray} \tag{ 11 }$$

and by ${{\hat{H}}_{F}}$ the (final) interacting one

$$\begin{eqnarray}&&{{\hat{H}}_{t\gt 0}}={{\hat{H}}_{I}}+\hbar \;g\;{{\hat{a}}^{\dagger }}\hat{a}(\hat{b}+{{\hat{b}}^{\dagger }})\equiv {{\hat{H}}_{F}}.\end{eqnarray} \tag{ 12 }$$

It is straightforward to prove that

$$\begin{eqnarray}\begin{array}{rcl} {{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\hat{H}}}_{F}}u}} & = & {{{\rm e}}^{-{\rm i}{{\omega }_{c}}u\,{{{\hat{a}}}^{\dagger }}\hat{a}+{\rm i}\frac{{{g}^{2}}}{\omega _{m}^{2}}({{\omega }_{m}}u-{\rm sin} {{\omega }_{m}}u){{({{{\hat{a}}}^{\dagger }}\hat{a})}^{2}}}}{{{\rm e}}^{-\frac{g}{{{\omega }_{m}}}{{{\hat{a}}}^{\dagger }}\hat{a}(\eta {{{\hat{b}}}^{\dagger }}-{{\eta }^{*}}\hat{b})}}{{{\rm e}}^{-{\rm i}{{\omega }_{m}}u\,{{{\hat{b}}}^{\dagger }}\hat{b}}}, \\ \end{array}\end{eqnarray} \tag{ 13 }$$

where $\eta =(1-{{{\rm e}}^{-{\rm i}{{\omega }_{m}}u}})$ [17]. 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 Kerr-like term.

The characteristic function in equation (10) can then be explicitly worked out. The form of the interaction suggests taking the trace over the number states ${{\{|n\rangle }_{c}}\}$ for mode $\hat{a}$ and over the coherent states ${{\{|\alpha \rangle }_{m}}\}$ for $\hat{b}$ (we reserve Latin letters for Fock-state labels and Greek letters for coherent-state labels throughout), i.e.,

$$\begin{eqnarray}&&\chi (u)=\mathop{\sum }\limits_{n=0}^{\infty }{{\int }_{\mathbb{C}}}{{{\rm d}}^{2}}\alpha \;p_{n}^{(c)}\;{{\mathcal{P}}^{(m)}}(\alpha )\langle n,\alpha |{{{\rm e}}^{\frac{{\rm i}}{\hbar }{{{\hat{H}}}_{F}}u}}{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\hat{H}}}_{I}}u}}|n,\alpha \rangle ,\end{eqnarray} \tag{ 14 }$$

where ${{\mathcal{P}}^{(m)}}(\alpha )={\rm exp} (-|\alpha {{|}^{2}}/{{N}_{m}})/(\pi {{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,\alpha \rangle \equiv |n{{\rangle }_{c}}\;\otimes \;|\alpha {{\rangle }_{m}}$. It is possible to gather the following analytical expression for the characteristic function

$$\begin{eqnarray}&&\chi (u)=\mathop{\sum }\limits_{n=0}^{\infty }\frac{N_{c}^{n}{{{\rm e}}^{-\frac{{{g}^{2}}{{n}^{2}}}{\omega _{m}^{2}}[{\rm i}({{\omega }_{m}}u-{\rm sin} {{\omega }_{m}}u)+(1+2{{N}_{m}})(1-{\rm cos} {{\omega }_{m}}u)]}}}{{{(1+{{N}_{c}})}^{n+1}}},\end{eqnarray} \tag{ 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 equation (15) cannot be explicitly worked out, we evaluate the probability distribution of the work by calculating equation (2) directly. To do this, energy eigenvalues and eigenstates of ${{\hat{H}}_{I}}$ and ${{\hat{H}}_{F}}$ are required. As ${{\hat{H}}_{I}}$ is the free Hamiltonian of the uncoupled system, it satisfies the eigenvalue equation ${{\hat{H}}_{I}}|n,k\rangle ={{E}_{n,k}}|n,k\rangle$, where $|n,k\rangle =|n{{\rangle }_{c}}\;\otimes \;|k{{\rangle }_{m}}$, and ${{E}_{n,k}}=\hbar {{\omega }_{c}}n+\hbar {{\omega }_{m}}(k+\frac{1}{2})$. Owing to the fact that $[{{\hat{a}}^{\dagger }}\hat{a},{{\hat{H}}_{F}}]=0$, the post-quench Hamiltonian can be written as ${{\hat{H}}_{F}}=\oplus {{_{n=0}^{\infty }}^{{}}}{{\hat{H}}_{F,n}}$, where ${{\hat{H}}_{F,n}}=|n\rangle \langle n{{|}_{c}}\left[ \hbar \omega _{c}^{{}}n+\hbar {{\omega }_{m}}({{{\hat{b}}}^{\dagger }}\hat{b}+\frac{1}{2})+\hbar \;g\;n(\hat{b}+{{{\hat{b}}}^{\dagger }}) \right]$ refers to the Hamiltonian of the n-photon manifold. Each ${{\hat{H}}_{F,n}}$ can then be diagonalized using a displacement operator $\hat{D}(z)={\rm exp} (z{{\hat{b}}^{\dagger }}-{{z}^{*}}\hat{b})$ on the mechanical mode, whose amplitude we take conditioned to the photon number n [18]. Denoting the quantities referring to ${{\hat{H}}_{F,n}}$ with a prime we find the energy eigenstates, written in the energy eigenbasis of the initial Hamiltonian ${{\hat{H}}_{I}}$, $|n^{\prime} {{\rangle }_{c}}\;\otimes \;{{\hat{D}}^{\dagger }}(\frac{g\;n^{\prime} }{{{\omega }_{m}}})|k^{\prime} {{\rangle }_{m}}$, with eigenvalues ${{E}_{n^{\prime} ,k^{\prime} }}=\hbar {{\omega }_{c}}n^{\prime} +\hbar {{\omega }_{m}}(k^{\prime} +\frac{1}{2})-\hbar \frac{{{g}^{2}}}{{{\omega }_{m}}}n{{^{\prime} }^{2}}$. A pictorial view of pre- and post-quench eigenstates in the subspace at fixed number n of photons is sketched in figure 2. As stated by equation (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

$$\begin{eqnarray}\begin{array}{rcl} P(W) & = & \mathop{\sum }\limits_{n,n^{\prime} ,k,k^{\prime} }p_{n}^{(c)}p_{k}^{(m)}{{|}_{m}}{{\langle k^{\prime} |\hat{D}[(g/{{\omega }_{m}})n^{\prime} ]|k\rangle }_{m}}{{|}^{2}}\delta [W-({{E}_{n^{\prime} ,k^{\prime} }}-{{E}_{n,k}})]{{\delta }_{n,n^{\prime} }} \\ {} & = & \mathop{\sum }\limits_{n,k,k^{\prime} }p_{n}^{(c)}p_{k}^{(m)}\frac{k!}{k^{\prime} !}{{{\rm e}}^{-{{(g/{{\omega }_{m}})}^{2}}{{n}^{2}}}}{{[(g/{{\omega }_{m}})n]}^{2(k^{\prime} -k)}}{{\left\{ \mathcal{L}_{k}^{(k^{\prime} -k)}[{{(g/{{\omega }_{m}})}^{2}}{{n}^{2}}] \right\}}^{2}}\delta \{W-\hbar {{\omega }_{m}}[k^{\prime} -k-{{(g/{{\omega }_{m}})}^{2}}{{n}^{2}}]\}, \\ \end{array}\end{eqnarray} \tag{ 16 }$$

where $\mathcal{L}_{a}^{b}(x)$ are the generalized Laguerre polynomials coming from the evaluation of the overlap between pre- and post-quench mechanical oscillator eigenstates [19]. A comparison with equation (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.

Zoom In Zoom Out Reset image size

The probability distribution of the work, together with real and imaginary parts of the characteristic function, is shown in figure 3, for different values of N_c, N_m, and coupling strength. By differentiating the expression of characteristic function equation (15) and evaluating it in the origin, according to the prescription in equation (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 $\chi (u)$, shown in the inset of figure 3, which approaches u = 0 with zero derivative. Having access to the characteristic function also gives us information about the statistical moments of P(W); for instance, the variance of the distribution is given by

$$\begin{eqnarray}&&\langle {{W}^{2}}\rangle -{{\langle W\rangle }^{2}}={{\hbar }^{2}}\;{{g}^{2}}{{N}_{c}}(1+2{{N}_{c}})(1+2{{N}_{m}}).\end{eqnarray} \tag{ 17 }$$

As expected, this quantity increases both with respect to the intensity of the quench, as quantified by $g/{{\omega }_{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/{{\omega }_{m}}=0.2$, to the other two, both obtained for ${{N}_{c}}=0.1$ and ${{N}_{m}}=1$—thus varying the ratio ${{\omega }_{c}}/{{\omega }_{m}}$—but corresponding to $g/{{\omega }_{m}}=0.1$ and $g/{{\omega }_{m}}=0.8$ respectively, i.e., increasing both the temperature and the coupling strength.

Zoom In Zoom Out Reset image size

Let us first analyze P(W) as illustrated for a few representative cases in figure 3, where we consider small values of $g/{{\omega }_{m}}\lesssim 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 $\hbar {{\omega }_{m}}$. These peaks originate from different initially-populated Fock states of the mechanical subsystem. Indeed, the number of peaks with appreciable amplitude increases strongly with N_m. In figure 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 $\hbar {{\omega }_{m}}$ and are biased towards less positive values of W. This is directly caused by the Kerr-like term in ${{\hat{H}}_{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

$$\begin{eqnarray}\begin{array}{rcll} \gamma & = & \frac{\langle {{(W-\langle W\rangle )}^{3}}\rangle }{{{\langle {{(W-\langle W\rangle )}^{2}}\rangle }^{3/2}}}= & \frac{{{\omega }_{m}}/g}{(1+2{{N}_{m}})\sqrt{{{N}_{c}}(1+2{{N}_{c}})}}, \\ \end{array}\end{eqnarray} \tag{ 18 }$$

and is more apparent in the low-temperature regime; indeed, by simply looking at the distribution shown in figure 3(b), it is possible to appreciate the positive skewness of the distribution.

Shifting our attention from figures 3 to 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 figure 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\;\hbar {{\omega }_{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.

Zoom In Zoom Out Reset image size

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 $\int _{-\infty }^{W}{\rm d}w\;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 figure 4. Another way of getting a smooth, continuous distribution would be to take the classical limit for the dynamics of both the optical and mechanical system, which would then result in a stochastically frequency-shifted harmonic oscillator [20]. While this approach would not entail the enforcement of any smoothing procedure, its study would deserve a full-fledged analysis that goes beyond the scopes of this work.

We have shown, both analytically and numerically, that despite turning on a nonlinear interaction between the two modes, on average there is no net production of work. This feature contrasts both with the case of a quench in the frequency of the harmonic potential of a single oscillator [21], and of the linear interaction between two bosonic modes [10], where 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 non-zero average work.

Using equation (13) we can easily compute the evolution of the initial Gibbs state, as defined by $\hat{\varrho }(t)={{{\rm e}}^{-\frac{{\rm i}}{\hbar }t{{{\hat{H}}}_{F}}}}\hat{\varrho }_{\beta }^{(c)}\;\otimes \;\hat{\varrho }_{\beta }^{(m)}{{{\rm e}}^{\frac{{\rm i}}{\hbar }t{{{\hat{H}}}_{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\pi r/{{\omega }_{m}}$ ($r\in \mathbb{Z}$), the system goes back to the initially factorized state, while for $t=(2r+1)\pi /{{\omega }_{m}}$ ($r\in \mathbb{Z}$), one gets the maximally (classically) correlated state.

Equation (13) also allows us to compute the partition function of the system, via a suitable Wick rotation of the argument, i.e., $u\to -{\rm i}\hbar \beta$, 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 ${{\mathcal{Z}}_{I}}=\mathcal{Z}_{\beta }^{(c)}\mathcal{Z}_{\beta }^{(m)}\equiv {{[(1-{{{\rm e}}^{-\hbar \beta {{\omega }_{c}}}})(1-{{{\rm e}}^{-\hbar \beta {{\omega }_{m}}}})]}^{-1}}$, while for the coupled system we obtain

$$\begin{eqnarray}&&{{\mathcal{Z}}_{F}}={{(1-{{{\rm e}}^{-\hbar \beta {{\omega }_{m}}}})}^{-1}}\mathop{\sum }\limits_{n=0}^{\infty }{{{\rm e}}^{-\hbar \beta {{\omega }_{c}}n}}{{{\rm e}}^{\hbar \beta ({{g}^{2}}/{{\omega }_{m}}){{n}^{2}}}}.\end{eqnarray} \tag{ 19 }$$

The free energy difference is correspondingly given by

$$\begin{eqnarray*}\begin{array}{rcll} \Delta F & = & -\frac{1}{\beta }{\rm ln} \left[ \mathop{\sum }\limits_{n=0}^{\infty }\frac{N_{c}^{n}}{{{(1+{{N}_{c}})}^{n+1}}}{{{\rm e}}^{\hbar \beta ({{g}^{2}}/{{\omega }_{m}}){{n}^{2}}}} \right]= & -\frac{1}{\beta }{\rm ln} \left[ 1-{{{\rm e}}^{-\hbar \beta {{\omega }_{c}}}} \right]-\frac{1}{\beta }{\rm ln} \left[ \mathop{\sum }\limits_{n=0}^{\infty }{{{\rm e}}^{-\hbar \beta {{\omega }_{c}}n}}{{{\rm e}}^{\hbar \beta \frac{{{g}^{2}}{{n}^{2}}}{{{\omega }_{m}}}}} \right], \\ \end{array}\end{eqnarray*}$$

which, as can be verified, agrees with the Jarzynski equality $\Delta F=-\frac{1}{\beta }{\rm ln} \chi ({\rm i}\beta )$. 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/{{\omega }_{m}}$, and ${{\omega }_{c}}/{{\omega }_{m}}$, there exists $\bar{n}=\bar{n}(g,r)$ such that $\forall n\gt \bar{n}$, we have that ${{g}^{2}}n\gt r$. This causes the sum to diverge exponentially, such that $\Delta F$ is formally undefined. This divergent term can be traced back to the part of ${{\hat{H}}_{F}}$ that reads ${{\omega }_{c}}{{\hat{a}}^{\dagger }}\hat{a}-{{g}^{2}}/{{\omega }_{m}}{{({{\hat{a}}^{\dagger }}\hat{a})}^{2}}$. As is apparent, the spectrum of this Hamiltonian is not bounded from below. Occupation of levels with $n\geqslant \bar{n}$, which occurs naturally for any non-zeroβ, can thus be mapped into a negative temperature with respect to ${{\hat{H}}_{F}}$. To resolve this issue, we impose a cutoff on the number of terms in the series; when $g/{{\omega }_{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 $\Delta F$, upon quenching the system. Alternatively, one can also devise ad hoc mechanisms to stabilize the divergent character of free energy difference. For example, one may add to the Hamiltonian ${{\hat{H}}_{I}}$ a quadratic term in the photon number operator $\hbar \;\epsilon {{({{\hat{a}}^{\dagger }}\hat{a})}^{2}}$; this can be done by adding a Kerr-like medium into the cavity, along the lines of [22]. When $\epsilon \geqslant {{g}^{2}}/{{\omega }_{m}}$ this contribution will compensate the divergent one, ensuring the converge of $\Delta F$. Another possibility is quenching both the linear and the quadratic therms in the mechanical displacement i.e., retaining both the contributions k = 1, 2 in equation (9); for the mechanical oscillator this would imply—beside a displacement and an energy shift both dependent on n—also a renormalization of the frequency (dependent on the photon number as well, see section 3.3) which eventually allow for the convergence of series for suitable values of the parameters involved. However it must be stressed that in both cases the regularization comes at the expense of modifying the probability distribution of the work and henceforth, for the rest of this work, we will therefore restrict ourselves to the physical domain in which the series does converge, leaving a quantitative study of the modification entailed by the additions mentioned above to further studies.

An explicit calculation of $\Delta F$, as illustrated in figure 5, shows that the free energy difference is negative, in agreement with the statement of the second law $\Delta F\leqslant \langle W\rangle \equiv 0$. Moreover, the irreversible work reduces to ${{W}_{{\rm irr}}}=-\Delta 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/{{\omega }_{m}}$, the free energy difference grows in modulus.

Zoom In Zoom Out Reset image size

Finally we point out that, in the spirit of [23], the divergence of the free energy difference is related, via the Jarzynski equality equation (7), to the divergence of exponentiated average work ${{\int }_{\mathbb{R}}}{\rm d}WP(W){{{\rm e}}^{-\beta W}}$, which in turn imposes some constraint on the behavior of the tails of the work probability distribution. In our case a divergent $\Delta F$ would imply a sub-exponential decay of the tails of the distribution.

In the previous subsection we observed how $\langle W\rangle =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 $\hat{x}$, is an odd function in the mechanical field operators, such that

$$\begin{eqnarray}&&\langle W\rangle =-g\;{{N}_{c}}{\rm Tr}\left\{ (\hat{b}+{{{\hat{b}}}^{\dagger }})\hat{\varrho }_{\beta }^{(m)} \right\}=0.\end{eqnarray} \tag{ 20 }$$

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 $\mathcal{E}\;{{\omega }_{m}}\in \mathbb{R}$ to the mechanical mode $\hat{b}$ of the Hamiltonian (9) so that the initial and final Hamiltonians will now read ${{\hat{H}}_{I,F,\mathcal{E}}}={{\hat{H}}_{I,F}}+\hbar \;\mathcal{E}\;{{\omega }_{m}}(\hat{b}+{{\hat{b}}^{\dagger }})$. It can be shown that ${{\hat{H}}_{I,\mathcal{E}}}=\hat{D}(\mathcal{E}){{\hat{H}}_{I}}{{\hat{D}}^{\dagger }}(\mathcal{E})$ and ${{\hat{H}}_{F,\mathcal{E}}}=\hat{D}(\mathcal{E})({{\hat{H}}_{F}}+2\hbar \;g\;\mathcal{E}{{\hat{a}}^{\dagger }}\hat{a}){{\hat{D}}^{\dagger }}(\mathcal{E})$ with $\hat{D}(\mathcal{E})$ a local displacement of amplitude $\mathcal{E}$. Proceeding as before, the characteristic function of the work distribution can be computed as

$$\begin{eqnarray}\begin{array}{rcl} \chi (u) & = & \mathop{\sum }\limits_{n=0}^{\infty }\frac{N_{c}^{n}}{{{(1+{{N}_{c}})}^{n+1}}}{{{\rm e}}^{-{\rm i}{{(g/{{\omega }_{m}})}^{2}}{{n}^{2}}({{\omega }_{m}}u-{\rm sin} {{\omega }_{m}}u)}}{{{\rm e}}^{-{{(g/{{\omega }_{m}})}^{2}}{{n}^{2}}(1+2{{N}_{m}})(1-{\rm cos} {{\omega }_{m}}u)}}{{{\rm e}}^{2{\rm i}gn\mathcal{E}u}}, \\ \end{array}\end{eqnarray} \tag{ 21 }$$

which differs from equation (15) by a phase factor. This extra factor is actually responsible for positive derivative of the imaginary part ${\rm Im}[\chi (u,\mathcal{E})]$ at the origin and hence to a non-zero value of the average work. Indeed, applying equation (4), one finds that the average work done by quenching the optomechanical interaction is given by

$$\begin{eqnarray}&&\langle W\rangle =2\hbar \;g\;\mathcal{E}{{N}_{c}},\end{eqnarray} \tag{ 22 }$$

which depends linearly on the displacement $\mathcal{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

$$\begin{eqnarray}&&\Delta F=-\frac{1}{\beta }{\rm ln} \left[ \mathop{\sum }\limits_{n=0}^{\infty }\frac{N_{c}^{n}}{{{(1+{{N}_{c}})}^{n+1}}}{{{\rm e}}^{\hbar \beta \frac{{{g}^{2}}{{n}^{2}}}{{{\omega }_{m}}}-2\hbar \beta gn\mathcal{E}}} \right].\end{eqnarray} \tag{ 23 }$$

The behavior of the irreversible work W_irr is reported in figure 6, with respect to the inverse temperature and the magnitude of the displacement.

Zoom In Zoom Out Reset image size

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)}}=\kappa /x_{{\rm zpf}}^{2}$, in analogy with the linear case. The initial Hamiltonian H_I is unmodified and still given by equation (11), whereas the the post-quench Hamiltonian now reads

$$\begin{eqnarray}&&{{\hat{H}}_{F}}={{\hat{H}}_{I}}+\hbar \;\kappa {{\hat{a}}^{\dagger }}\hat{a}{{(\hat{b}+{{\hat{b}}^{\dagger }})}^{2}}.\end{eqnarray} \tag{ 24 }$$

We choose to work with a non-negative κ, since $\kappa \lt 0$ can introduce post-quench instabilities similar to the one noted for the linear case. The $\kappa \gt 0$ case exhibits no such instabilities. Yet again, we see that this interaction preserves the photon number ${{\hat{a}}^{\dagger }}\hat{a}$, so that it proves convenient to write ${{\hat{H}}_{F}}=\oplus {{_{n=0}^{\infty }}^{{}}}{{\hat{H}}_{F,n}}$ where each ${{\hat{H}}_{F,n}}$ can be cast in the form

$$\begin{eqnarray}&&{{\hat{H}}_{F,n}}=\left[ \hbar {{\omega }_{c}}n+\hbar {{\Omega }_{n}}\left( {{{\hat{b}}}^{\dagger }}\hat{b}+\frac{1}{2} \right)+\hbar {{\Sigma }_{n}}\left( {{{\hat{b}}}^{\dagger \,2}}+{{{\hat{b}}}^{2}} \right) \right]|n\rangle \langle n{{|}_{c}},\end{eqnarray} \tag{ 25 }$$

where ${{\Omega }_{n}}\equiv {{\omega }_{m}}+2\;\kappa n$ and ${{\Sigma }_{n}}\equiv 2\kappa 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

$$\begin{eqnarray}&&{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\hat{H}}}_{F}}u}}=\mathop{\sum }\limits_{n=0}^{\infty }{{{\rm e}}^{-{\rm i}u[{{\omega }_{c}}n+{{\Omega }_{n}}({{{\hat{b}}}^{\dagger }}\hat{b}+\frac{1}{2})+{{\Sigma }_{n}}({{{\hat{b}}}^{\dagger \,2}}+{{{\hat{b}}}^{2}})]}}|n\rangle \langle n{{|}_{c}}.\end{eqnarray} \tag{ 26 }$$

Our next task is to disentangle each exponential operator in the sum. By using the commutation relations between the operators involved in equation (26), which provide a two-excitation realization of the $\mathfrak{s}\mathfrak{u}(1,1)$ algebra [24], we find

$$\begin{eqnarray}&&{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\hat{H}}}_{F,n}}u}}={{{\rm e}}^{\frac{1}{2}[\xi _{n}^{*}{{{\hat{b}}}^{2}}-{{\xi }_{n}}{{{\hat{b}}}^{\dagger 2}}]}}{{{\rm e}}^{-{\rm i}{{\eta }_{n}}({{{\hat{b}}}^{\dagger }}\hat{b}+\frac{1}{2})}}{{{\rm e}}^{-{\rm i}{{\omega }_{c}}u\,n}}|n\rangle \langle n{{|}_{c}},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&{{\eta }_{n}}\equiv {\rm arctan} \left[ \frac{1+2\tilde{\kappa }n}{\sqrt{1+4\tilde{\kappa }n}}{\rm tan} \left( {{\omega }_{m}}u\sqrt{1+4\tilde{\kappa }n} \right) \right]\end{eqnarray} \tag{ 28 }$$

with $\tilde{\kappa }\equiv \kappa /{{\omega }_{m}}$ being a dimensionless quench parameter. We further have the complex quantity ${{\xi }_{n}}\equiv |{{\xi }_{n}}|{{{\rm e}}^{{\rm i}{{\phi }_{n}}}}$ whose phase is ${{\phi }_{n}}\equiv {{\eta }_{n}}+\frac{\pi }{2}$ and modulus

$$\begin{eqnarray}&&|{{\xi }_{n}}|\equiv {\rm arcsinh}\left[ \frac{2\tilde{\kappa }n}{\sqrt{1+4\tilde{\kappa }n}}{\rm sin} \left( {{\omega }_{m}}u\sqrt{1+4\tilde{\kappa }n} \right) \right].\end{eqnarray} \tag{ 29 }$$

Armed with this tool we can thus compute the characteristic function of the work distribution, which reads

$$\begin{eqnarray}&&\chi (u)=\mathop{\sum }\limits_{n=0}^{\infty }\frac{N_{c}^{n}}{{{(1+{{N}_{c}})}^{n+1}}}\frac{1}{\sqrt{\mathop{\sum }\limits_{j=0}^{2}{{\chi }_{n,j}}N_{m}^{j}}},\end{eqnarray} \tag{ 30 }$$

and comes in the form of a thermal average with respect to the cavity distribution—as in equation (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 ${{\chi }_{n,0}}$

into its real and imaginary parts, which read

$$\begin{eqnarray}\begin{array}{rcl} \operatorname{Re}({{\chi }_{n,0}}) & = & {\rm cos} ({{\omega }_{m}}u){\rm cos} ({{\omega }_{m}}u\sqrt{1+4\tilde{\kappa }n})+\frac{1+2\tilde{\kappa }n}{\sqrt{1+4\tilde{\kappa }n}}{\rm sin} ({{\omega }_{m}}u){\rm sin} \left( {{\omega }_{m}}u\sqrt{1+4\tilde{\kappa }n} \right), \\ \end{array}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}\begin{array}{rcl} {\rm Im}({{\chi }_{n,0}}) & = & {\rm sin} ({{\omega }_{m}}u){\rm cos} \left( {{\omega }_{m}}u\sqrt{1+4\tilde{\kappa }n} \right)-\frac{1+2\tilde{\kappa }n}{\sqrt{1+4\tilde{\kappa }n}}{\rm cos} ({{\omega }_{m}}u){\rm sin} \left( {{\omega }_{m}}u\sqrt{1+4\tilde{\kappa }n} \right). \\ \end{array}\end{eqnarray} \tag{ 32 }$$

We thus have ${{\chi }_{n,1}}=2({{\chi }_{n,0}}-1)\ {\rm and}\ {{\chi }_{n,2}}=2\left[ \operatorname{Re}({{\chi }_{n,0}})-1 \right]$. As before, since the Fourier transform of equation (30) cannot be directly evaluated, in order to compute the probability distribution of the work equation (2) we proceed by diagonalizing the post-quench Hamiltonian ${{\hat{H}}_{F}}$. First, we keep in mind that ${{\hat{H}}_{I}}$ is the same as before. However, within any fixed photon number manifold, ${{\hat{H}}_{F,n}}$ be diagonalized via a squeezing operation $\hat{S}(z)={\rm exp} ({{z}^{*}}{{\hat{b}}^{2}}/2-z\;{{\hat{b}}^{\dagger \,2}}/2)$ on the mechanical mode conditioned on the photon number n [25, 26]. Once again denoting the post-quench quantities with a prime, and expressing the states in the eigenbasis of ${{\hat{H}}_{I}}$, we find eigenstates ${{\hat{H}}_{F,n}}|n^{\prime} {{\rangle }_{c}}\;\otimes \;\hat{S}({{\zeta }_{{{n}^{\prime }}}})|k^{\prime} {{\rangle }_{m}}={{E}_{n^{\prime} ,k^{\prime} }}|n^{\prime} {{\rangle }_{c}}\;\otimes \;\hat{S}({{\zeta }_{{{n}^{\prime }}}})|k^{\prime} {{\rangle }_{m}}$, where the squeezing parameter is given by ${{\zeta }_{{{n}^{\prime }}}}\equiv \frac{1}{4}{\rm log} \left[ 1+4(\kappa /{{\omega }_{m}}){{n}^{\prime }} \right]$, and the eigenvalue

$$\begin{eqnarray}&&{{E}_{n^{\prime} ,k^{\prime} }}=\hbar {{\omega }_{c}}n^{\prime} +\hbar {{\omega }_{m}}\sqrt{1+4(\kappa /{{\omega }_{m}}){{n}^{\prime }}}(k^{\prime} +\frac{1}{2}).\end{eqnarray} \tag{ 33 }$$

As sketched in figure 7, for the manifold corresponding to ${{n}^{\prime }}$ photons, the quench results in a modification of the oscillation frequency which, is multiplied by a factor $\sqrt{1+4(\kappa /{{\omega }_{m}}){{n}^{\prime }}}$, a relative shift of the mechanical levels by $\hbar {{\omega }_{m}}\left[ \sqrt{1+4(\kappa /{{\omega }_{m}}){{n}^{\prime }}}-1 \right]$, and a squeezing of the state by a factor ${{\zeta }_{{{n}^{\prime }}}}$. Putting everything together, the probability distribution of the work is thus given by

$$\begin{eqnarray}\begin{array}{rcl} P(W) & = & \mathop{\sum }\limits_{n,n^{\prime} ,k,k^{\prime} }p_{n}^{(c)}p_{k}^{(m)}{{\left| \langle k^{\prime} |S({{\zeta }_{{{n}^{\prime }}}})|k\rangle \right|}^{2}}\;\delta \left[ W-({{E}_{n^{\prime} ,k^{\prime} }}-{{E}_{n,k}}) \right]{{\delta }_{n,n^{\prime} }} \\ {} & = & \mathop{\sum }\limits_{n,k,k^{\prime} }p_{n}^{(c)}p_{k}^{(m)}\frac{k!\;k^{\prime} !}{{{({\rm cosh} {{\zeta }_{n}})}^{2k+1}}}{{\left[ \mathcal{S}(k,k^{\prime} ,{{\zeta }_{n}}) \right]}^{2}}\delta \left\{ W-\hbar {{\omega }_{m}}\left[ \sqrt{1+\frac{4n\kappa }{{{\omega }_{m}}}}k^{\prime} -k \right] \right\}, \\ \end{array}\end{eqnarray} \tag{ 34 }$$

where $\mathcal{S}(k,k^{\prime} ,{{\zeta }_{n}})$ is given by

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{S}(k,k^{\prime} ,{{\zeta }_{n}}) & = & \mathop{\sum }\limits_{m=0}^{\lfloor \frac{k^{\prime} }{2}\rfloor }\mathop{\sum }\limits_{l=0}^{\lfloor \frac{k}{2}\rfloor }\frac{{{(-1)}^{3m+2l}}}{{{2}^{m+l}}m!\;l!}\frac{{{({\rm tanh} {{\zeta }_{n}})}^{m+l}}}{(k-2l)!}{{({\rm cosh} {{\zeta }_{n}})}^{2l}}{{\delta }_{k^{\prime} -2m,k-2l}}, \\ \end{array}\end{eqnarray} \tag{ 35 }$$

being $\lfloor x\rfloor$ the floor function of argument x, which yields the largest integer not greater than x.

Zoom In Zoom Out Reset image size

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 figures 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

$$\begin{eqnarray}&&\langle W\rangle =\hbar \kappa {{N}_{c}}(1+2{{N}_{m}}),\end{eqnarray} \tag{ 36 }$$

hence increasing with respect the occupation numbers of both the cavity and the mechanical mode, as made apparent by inspecting the different panels in figure 8. The variance of the distribution reads

$$\begin{eqnarray}&&\langle {{W}^{2}}\rangle -{{\langle W\rangle }^{2}}={{\hbar }^{2}}{{\kappa }^{2}}{{N}_{c}}(3+5{{N}_{c}}){{(1+2{{N}_{m}})}^{2}}.\end{eqnarray} \tag{ 37 }$$

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

$$\begin{eqnarray}&&\gamma =\frac{4+8{{N}_{c}}+(g/{{\omega }_{m}})(15+81{{N}_{c}}+74N_{c}^{2}){{(1+2{{N}_{m}})}^{2}}}{(g/{{\omega }_{m}})\sqrt{{{N}_{c}}}{{(3+5{{N}_{c}})}^{3/2}}{{(1+2{{N}_{m}})}^{2}}}.\end{eqnarray} \tag{ 38 }$$

We note that, for ${{N}_{m}}\gg 1$, it acquires the values $5/\sqrt{3{{N}_{c}}}$ for ${{N}_{c}}\ll 1$ and $74/5\sqrt{5}$ for ${{N}_{c}}\gg 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 $\hat{\varrho }(t)={{{\rm e}}^{-\frac{{\rm i}}{\hbar }t{{{\hat{H}}}_{F}}}}\hat{\varrho }_{\beta }^{(c)}\;\otimes \;\hat{\varrho }_{\beta }^{(m)}{{{\rm e}}^{\frac{{\rm i}}{\hbar }t{{{\hat{H}}}_{F}}}}={{\sum }_{n}}p_{n}^{(c)}|n\rangle \langle n{{|}_{c}}\;\otimes \;\int {{{\rm d}}^{2}}\alpha \;{{\mathcal{P}}^{(m)}}(\alpha )|{{{\rm e}}^{{\rm i}{{\eta }_{n}}}}\alpha ,{{\xi }_{n}}\rangle \langle {{{\rm e}}^{{\rm i}{{\eta }_{n}}}}\alpha ,{{\xi }_{n}}{{|}_{m}}$, where $|{{{\rm e}}^{{\rm i}{{\eta }_{n}}}}\alpha ,{{\xi }_{n}}{{\rangle }_{m}}=\hat{S}({{\xi }_{n}})\hat{D}({{{\rm e}}^{{\rm i}{{\eta }_{n}}}}\alpha )|0{{\rangle }_{m}}$ is a squeezed coherent state of the mechanical mode, and hence no 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

$$\begin{eqnarray}\begin{array}{rcl} \Delta F & = & -\frac{1}{\beta }{\rm ln} \left[ {\rm sinh} \left( \frac{\beta }{2} \right) \right]-\frac{1}{\beta }{\rm ln} \left[ \mathop{\sum }\limits_{n=0}^{\infty }\frac{N_{c}^{n}\;{\rm cosech}\left( \sqrt{1+\frac{4n\kappa }{{{\omega }_{m}}}}\frac{\beta }{2} \right)}{{{(1+{{N}_{c}})}^{n+1}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

In this case, too, a suitable Wick-like rotation to imaginary u can be performed to obtain $\Delta F$ from $\chi (u)$. In practice, however, this calculation is frought with technical difficulties and it is far easier to compute $\Delta F$ from an explicit diagonalization of the Hamiltonian, as was done above. The behavior of the irreversible work for this case has been shown in figure 10, and once again we can see how it drops lowering the temperature and increases by increasing the coupling strength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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}^{\prime }}=k$, therefore, which at least for small $\kappa /{{\omega }_{m}}$ has a large probability of occurring, the work done is positive.