\begin{document}

    
\title{Optomechanical cooling with coherent and squeezed light: the thermodynamic cost of opening the heat valve}

\author{Juliette Monsel}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}

\author{Nastaran Dashti}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}

\author{Sushanth Kini Manjeshwar}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}

\author{Jakob Eriksson}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}
\affiliation{University of Gothenburg, S-412 96 G\"oteborg, Sweden}
\author{Henric Ernbrink}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}
\author{Ebba Olsson}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}
\author{Emelie Torneus}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}

\author{Witlef Wieczorek}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}
\author{Janine Splettstoesser}
\affiliation{Department of Microtechnology and Nanoscience (MC2), Chalmers University of Technology, S-412 96 G\"oteborg, Sweden\looseness=-1}

\delta\hat{a}te{\today}                                                                                                                                                                                                                                                 
\begin{abstract}
Ground-state cooling of mechanical motion by coupling to a driven optical cavity has been demonstrated in various optomechanical systems. In our work, we provide a so far missing thermodynamic performance analysis of optomechanical sideband cooling in terms of a heat valve. As performance quantifiers we examine not only the lowest reachable effective temperature (phonon number), but also the evacuated-heat flow as an equivalent to the cooling power of a standard refrigerator, as well as appropriate thermodynamic efficiencies, which all can  be experimentally inferred from measurements of the cavity output light field. Importantly, in addition to the standard optomechanical setup fed by coherent light, we investigate two recent alternative setups for achieving ground-state cooling: replacing the coherent laser drive by squeezed light or using a cavity with a frequency-dependent (Fano) mirror. We study the dynamics of these setups within and beyond the weak-coupling limit and
give concrete examples based on parameters of existing experimental systems. By applying our thermodynamic framework, we gain detailed insights into these three different optomechanical cooling setups, allowing a comprehensive understanding of the thermodynamic mechanisms at play.
\text{e}nd{abstract}

\text{mec}aketitle	

\text{s}ection{Introduction}

Nano- and micromechanical resonators constitute an excellent platform for exploring thermodynamics on the nanoscale~\cite{Seifert2012Nov,lutz_information_2015,ciliberto_experiments_2017}. Importantly, stochastic and quantum fluctuations can play a major role determining nanomechanical motion~\cite{klaers_squeezed_2017,oconnell_quantum_2010}, enabling tests of stochastic~\cite{Seifert2012Nov} and quantum thermodynamics concepts~\cite{QuThermo} in experiments. In this respect, optomechanical systems~\cite{Aspelmeyer2014Dec} are a pertinent platform as they allow for precise control over the classical and quantum dynamics of nanomechanical motion by using electromagnetic fields. This capability motivated proposals of optomechanics-based work measurement schemes~\cite{Elouard2015May, Dong2015Sep, Monsel2018Nov} or heat engines~\cite{Zhang2014Apr, Zhang2014Aug, gelbwaser-klimovsky_work_2015, Dechant2015May, Mari2015Jul, Bathaee2016Aug, Dechant2017Sep, Bennett2020Oct} and already lead to the experimental quantification of non-equilibrium thermodynamic processes such as irreversible entropy production  \cite{brunelli_experimental_2018,rossi_experimental_2020}.

An ubiquitous thermodynamic process is cooling. In the context of optomechanics, cooling is employed to reduce the entropy of eigenmodes of the mechanical resonator, which is a necessary step for exerting quantum control over mechanical motion and, thus, a crucial requirement for the implementation of any optomechanics-based quantum technology. Furthermore, any optomechanics-based heat engine exploits cooling in its operation cycle. Cooling to the ground state of mechanical motion has been theoretically analyzed (see, e.g., Refs.~\cite{Marquardt2007Aug, wilson-rae_theory_2007, Genes2008Mar, Wilson-Rae2008Sep}) and experimentally realized by direct thermal cooling~\cite{oconnell_quantum_2010}, optomechanical sideband cooling~\cite{Chan2011,teufel_sideband_2011,Delic2020} or feedback cooling~\cite{rossi_measurement-based_2018,magrini_optimal_2020}. In the context of thermodynamics, optomechanical cooling has been examined in the light of entropy production \cite{brunelli_experimental_2018,rossi_experimental_2020} or as a quantum absorption refrigerator \cite{mitchison_realising_2016, Mitchison2018, Naseem2020Jun} fed by thermal light.
However, a detailed analysis of the desired output, such as the reduction in phonon number, related to the \textit{thermodynamic cost} of optomechanical cooling with non-thermal light sources in the steady-state has so far been missing.  

\begin{figure}[htb!]
    \includegraphics[width=\linewidth]{fig_setup.pdf}             
    \caption{\label{fig:principle}
        \textbf{(a)} Thermodynamic framework:  a mechanical resonator (in black) is in contact with a hot phonon bath (in red). To cool it down, it is coupled to a photon bath (in blue) whose temperature is much smaller than the relevant photon energies. The coupling can be regarded as opening a valve (in green), allowing heat to flow to the cold bath. The arrows indicate the chosen direction for positive heat flows, but they do not reflect the actual direction of the heat flows in the setup.

        \textbf{(b)} Standard setup for optomechanical cooling: a cavity with a moving-end mirror is driven by a laser. The mechanical resonator is in contact with a thermal phonon reservoir at temperature $T_\text{mec}$ and the cavity with a photon reservoir at temperature $T_\text{cav}$. The light blue dashed box represents the effective phonon bath that can be identified in the weak-coupling regime. 
        See Sec.~\ref{sec:II:standard setup} for definitions of the notations. 
    }
\text{e}nd{figure}

In our work, we present this missing analysis of the thermodynamic performance of different optomechanical sideband cooling schemes, which is of crucial relevance to optimize the performance of future devices.
We therefore provide a comprehensive theoretical framework for sideband cooling beyond previously performed approximations.
We find that this analysis not only yields information about actual performance benchmarks, but also provides insights into the thermodynamic mechanisms at play and the required resources for realizing optomechanical cooling in different experimental setups. 
First, we notice that optomechanical cooling does not correspond to the commonly employed refrigerator framework, in which the work provided to a system makes heat flow out of a cold bath\footnote{See also Appendix \ref{appendix:comp absorption refrigerator} for an overview of the differences between our framework and absorption refrigerators.}. 
The thermodynamic system to be cooled consists of a mechanical resonator (from which heat should be taken away) that is in contact with a hot phonon bath, which leads to Brownian motion of the resonator. On the other hand, there is a photon bath, which is determined by the noise properties of the light field and which appears as cold, since its temperature is much smaller than the energy scale of the photons. This photon bath couples to the mechanical degrees of freedom via a driven optomechanical cavity [see Fig.~\ref{fig:principle}]. We identify that the light field acts as a tunable valve that controls the flow of heat from the mechanical resonator to the cold bath, i.e., the noise of the light field. Therefore, unlike as in usual refrigerators, heat flows from hot to cold, but, importantly, there is a cost associated with maintaining the valve open (via the driving of the cavity). The heat-valve scenario that we suggest here thereby gives a very clear picture of this setting with two heat baths and an additional work input. It is particularly helpful for the understanding of the mechanisms at play when comparing sideband cooling in different optomechanical setups.

We furthermore find it insightful to look at optomechanical cooling as a special example of reservoir engineering \cite{Poyatos1996Dec}: a low temperature \textit{photon} bath is prepared with the help of a laser-irradiated cavity, such that it is transformed into an effective cold \textit{phonon} bath \cite{Marquardt2007Aug}. Engineered reservoirs have shown to be especially appealing to thermodynamics as they allow to explore non-thermal resources~\cite{Brunner2012May, Abah2014Apr, Niedenzu2018Jan, Ghosh2018, Sanchez2019Nov} and study the impact of genuine quantum features, such as squeezing~\cite{Correa2014Feb, Manzano2016May, Agarwalla2017Sep, Manzano2018Oct}, entanglement~\cite{Francica2017Mar} or quantum coherence~\cite{Scully2003Feb, Harrington2019May, Monsel2020Mar} on energy and entropy exchanges. Recently, it was demonstrated that a classically engineered reservoir can be used to beat the classical Carnot bound of heat machines~\cite{klaers_squeezed_2017}. 
Here, the situation is different: we start from a photon bath, in which thermal fluctuations are negligible at optical frequencies, in order to engineer a \textit{thermal} bath of a completely different type, namely an effective phonon bath. In this sense the \textit{valve} that is established is an important ingredient of the engineered bath.
In most cases, the cost to create the engineered bath is not taken into account in the efficiency. 

In our work, we provide and calculate appropriate cooling efficiencies of optomechanical sideband cooling, which account for the cost of the heat valve (i.e.~the bath engineering) and which yield an additional benchmark in optomechanics, besides the resonator's phonon occupation. We then apply our comprehensive thermodynamic analysis of optomechanical cooling to three different setups, which can be used to cool mechanical motion to the ground state. These are: (i) a standard optomechanical setup driven by a coherent laser input [Fig.~\ref{fig:principle}], (ii) a standard optomechanical setup driven by \textit{squeezed-light} [Fig.~\ref{fig:Fano_and_squeezing}(a)], and (iii) an optomechanical cavity that employs a frequency-dependent (Fano) mirror as one of the cavity boundaries and that is driven by a coherent laser [Fig.~\ref{fig:Fano_and_squeezing}(b)]. Importantly, we illustrate our analysis with concrete experimental parameters, which allow us to benchmark different optomechanical platforms against each other. Setup (i)~\cite{Marquardt2007Aug, wilson-rae_theory_2007, Genes2008Mar, Wilson-Rae2008Sep} is widely employed for ground-state cooling provided the optomechanical system is in the resolved-sideband regime~\cite{Chan2011,teufel_sideband_2011}. Setup (ii)~\cite{Asjad2016Nov,Clark2017Jan,Lau2020Mar} is motivated by shaping the quadrature noise of the cold photon bath. This setup can lead to a reduced backaction of the light field onto the mechanical resonator and, thus, to ground-state cooling in a regime, where a coherent drive fails, as we here illustrate. Finally, setup (iii)~\cite{Cernotik2019Jun} results in a Fano-like cavity line shape. This setup aims at suppressing the Stokes process and, thus, the light-induced heating of the mechanical motion originating from the optomechanical interaction itself. Importantly, our theoretical analysis, which employs a Langevin as well as a master equation approach, is valid both for weak and strong optomechanical coupling\footnote{We are referring here to the effective optomechanical coupling $g = g_0\abs{\alpha}$, enhanced by the light field in the cavity, of amplitude $\abs{\alpha}$, see Eq.~(\ref{cooperativity}). On the contrary, the single-photon coupling $g_0$ is assumed to be very small compared to the other frequencies involved.} It thereby extends the validity of previous work~\cite{Asjad2016Nov}.

This paper is structured as follows. In Sec.~\ref{sec:II} we describe the theoretical models we use for each setup, outline the thermodynamic picture that we develop here and present the specific experimental platforms we use as examples. In Sec.~\ref{sec:III} we derive the dynamics of the three setups and apply these results to a thermodynamic analysis of optomechanical cooling illustrated with realistic optomechanics experiments. We conclude in Sec.~\ref{sec:IV}. The appendix contains detailed information on the employed theoretical methods and approximations and thereby allows to straightforwardly follow all calculations performed in this work.

\text{s}ection{Setups and models} \label{sec:II}

\text{s}ubsection{Sideband cooling in cavity optomechanical setups}

In the following, we introduce the models for the three different setups that we consider for sideband cooling and start with the standard setup of driving an optomechanical cavity with coherent light.

\text{s}ubsubsection{Standard setup}\label{sec:II:standard setup}

A typical optomechanical setup consists of a Fabry-Pérot cavity of frequency $\omega_\text{cav}$ with a moving-end mirror of mechanical resonance frequency $\Omega_\text{mec}$ [see Fig.~\ref{fig:principle}(b)]. The cavity is driven by a laser of frequency $\omega_\text{las}$. It is described by the Hamiltonian~\cite{Aspelmeyer2014Dec}, in the frame rotating at $\omega_\text{las}$,
 \begin{align}\label{H_sys}
   \hat{H}=\,&\hbar \Omega_\text{mec} \hat{b}^\delta\hat{a}gger \hat{b} + \hbar \Delta_0 \hat{a}^\delta\hat{a}gger \hat{a} +\hbar (\text{e}psilon \hat{a}^\delta\hat{a}gger+ \text{e}psilon^{*} \hat{a})\nonumber\\
   &-\hbar g_0  \hat{a}^\delta\hat{a}gger \hat{a} (\hat{b} + \hat{b}^\delta\hat{a}gger) ,
\text{e}nd{align}
where $\Delta_0 = \omega_\text{cav} - \omega_\text{las}$ is the bare detuning between the cavity and the laser. We have denoted $\hat{a}$ the photon-annihilation operator of the cavity and $\hat{b}$ the phonon-annihilation operator. The single-photon optomechanical coupling strength is given by $g_0$ and $\text{e}psilon$ is the drive amplitude. 

The cavity is coupled to a photon and a phonon bath, where
$\kappa$ denotes the optical loss rate and $\gamma$ the mechanical damping rate. The evolution of this system can be described with a set of Langevin equations~\cite{Genes2008Mar}
\begin{subequations}\label{standard:Langevin}
    \begin{align}
        \dot{\hat q}&=\Omega_\text{mec} \hat p, \\
        \dot{\hat p}&=-\Omega_\text{mec} \hat q-\gamma \hat p+g_0 \text{s}qrt{2}\hat{a}^\delta\hat{a}gger \hat{a}  +\text{s}qrt{\gamma } \hat\xi, \\
        \dot{\hat a}&=-(\kappa+i\Delta_0 ) \hat{a} +ig_0 \text{s}qrt{2}\hat{a} \hat q -i \text{e}psilon +\text{s}qrt{2 \kappa} \hat{a}_\text{in},
    \text{e}nd{align}
\text{e}nd{subequations}
where we have defined the mechanical quadratures $\hat q=\left(\hat{b}+\hat{b}^\delta\hat{a}gger  \right)/\text{s}qrt{2}$ and $\hat p=\left(\hat{b}- \hat{b}^\delta\hat{a}gger  \right)/i \text{s}qrt{2}$. The thermal noise of the mechanics is given by the operator $\hat{\xi}$, while vacuum noise of the light field is associated with the operator $\hat{a}_\text{in}$. 

We assume in the following that the mechanical quality factor $Q_\text{mec} = \Omega_\text{mec}/\gamma$ is large and that the temperature $T_\text{mec}$ of the mechanical environment is high, that is $\hbar\Omega_\text{mec} \ll k_\text{B} T_\text{mec}$, see Sec.~\ref{sec:II:experimental setups} for concrete implementations in this regime. We can then make the white-noise approximation for the correlation function for the thermal noise of the mechanics
\begin{equation}
    \text{mec}ean{\hat \xi(t)\hat \xi(t')}=(2\bar{n}_\text{mec}+1) \delta (t-t'),\label{correlations_xi}
\text{e}nd{equation}
where $\bar{n}_\text{mec} = (\text{e}^{\hbar\Omega_\text{mec}/k_\text{B} T_\text{mec}} - 1)^{-1}$ is the average number of phonons in the mechanical environment.
This approximation is equivalent to the Born-Markov approximation made in the derivation of the master equation~\cite{Wilson-Rae2008Sep} (see Appendix~\ref{appendix:standard:approxs} for a comparison of the different approaches and level of approximations frequently used in optomechanics). In the examples studied in this article (presented in Sec.~\ref{sec:II:experimental setups}), this approximation yields satisfactory results, but we explain how to relax it in Appendix~\ref{appendix:standard:beyond white noise}.

In order to quantify the ratio between the optomechanical coupling and the coupling to the optical and mechanical baths, we use the quantum cooperativity,
\begin{equation}
    C=\frac{2g^2}{\kappa \gamma \bar{n}_\text{mec}}.\label{cooperativity}
\text{e}nd{equation}
Here, $g=g_0|\alpha|$ is the effective coupling with $|\alpha|$ being the amplitude of the coherent field in the cavity, which is proportional to the driving amplitude of the laser $\text{e}psilon$, see Eq.~\text{e}qref{standard:alpha}. 

\begin{figure}[bt]
    \includegraphics[width=\linewidth]{fig_ppe_cooling.pdf}             
    \caption{\label{fig:weak coupling cooling}
        Principle of resolved-sideband cooling in the weak-coupling limit in a Raman-scattering picture. The Stokes and anti-Stokes processes are represented by the magenta and cyan arrows, respectively. 
        \textbf{(a)} Cavity density of states. $\Delta$ is the effective detuning (defined in Sec.~\ref{sec:II:Thermodynamic picture}) between the laser frequency (orange line) and the cavity resonance frequency (dashed black line), which is shifted with respect to $\Delta_0$ by the optomechanical interaction.

        \textbf{(b)} Partial energy diagram of the optomechanical system around the state $\ket{p, n}$ with $p$ photons and $n$ phonons. 
    }
\text{e}nd{figure}

 
The laser driving the cavity is modeled as a classical drive, included in the Hamiltonian~\text{e}qref{H_sys}. Its drive amplitude $\text{e}psilon$ is related to the input laser power ${P}_\text{las}$ by \begin{equation}
    {P}_\text{las} = \hbar\omega_\text{las} \frac{\abs{\text{e}psilon}^2}{2\kappa}.\label{Pin}
\text{e}nd{equation}  
The cavity and the resonator are typically at the same temperature, i.e. $T_\text{mec} = T_\text{cav}$, however the cavity frequency is typically orders of magnitude larger than the mechanical frequency and $\hbar \omega_\text{cav} \gg k_\text{B} T_\text{cav}$ (see Sec.~\ref{sec:II:experimental setups}). Therefore the average number of photons in the optical environment is negligible,  $\bar{n}_\text{cav} \text{s}imeq 0$, and when driving with coherent light, there is only one non-zero correlation function of the cavity input noise, in contrast to the squeezed case, see Eqs.~(\ref{squeezing:correlations_a_in}). This correlation function is
\begin{equation}
     \text{mec}ean{\hat{a}_\text{in}(t)  \hat{a}_\text{in}^\delta\hat{a}gger (t')}= \delta (t-t'). \label{standard:correlations_a_in}
\text{e}nd{equation}

To sketch the principle of sideband cooling, we consider temporarily the resolved-sideband regime, that is $\kappa \ll \Omega_\text{mec}$, and a weak optomechanical coupling: $g \ll \kappa, \Omega_\text{mec}$. In this regime, as illustrated by Fig.~\ref{fig:weak coupling cooling}, three processes can occur: A photon from the laser can (1) be absorbed by the cavity without changing the state of the resonator (in orange), (2) create a photon in the cavity at lower frequency $\omega_\text{las} - \Omega_\text{mec}$ and a phonon (Stokes process, in magenta), or (3) in combination with a phonon, create a photon at higher frequency $\omega_\text{las} + \Omega_\text{mec}$ (anti-Stokes process, in cyan). The Stokes process heats up the resonator while the anti-Stokes process cools it down. By choosing a cavity-laser detuning close to the mechanical frequency, the rates of the first two processes can be reduced, allowing  effective cooling~\cite{Schliesser2008May}. Note that the results presented in this article focus on, but are not restricted to, the resolved-sideband regime and they are valid beyond the weak-coupling limit. 
 

\begin{figure}[bt]
    \includegraphics[width=\linewidth]{fig_Fano_squeezing.pdf}
    \caption{\label{fig:Fano_and_squeezing}
        \textbf{(a)} Squeezed-light setup: We add to the laser driving the cavity a squeezed vacuum field. Since the control of the environment of the cavity is not perfect, there is some residual vacuum noise entering the cavity. The squeezed noise is generated with a degenerate parametric amplifier (DPA) driven at $2\omega_\text{las}$ thanks to a second harmonic generator (SHG). See Sec.~\ref{sec:II:squeezed-light setup} for the definitions of the notations. 

        \textbf{(b)} Fano-mirror setup: The left-hand side mirror of the cavity is replaced by a Fano mirror. Unlike in the canonical setup [Fig.~\ref{fig:principle}(b)], we differentiate the noise coming through the left and right mirrors. See Sec.~\ref{sec:II:Fano-mirror setup} for the definitions of the notations. 

        \textbf{(c)} Cavity density of states (solid gray line) which is asymmetrical in the Fano-mirror setup. The dashed gray line represents the density of states of a standard cavity of equivalent linewidth. 

    }
\text{e}nd{figure}

\text{s}ubsubsection{Squeezed-light setup}\label{sec:II:squeezed-light setup}

To improve the cooling scheme presented above, one can use squeezed light instead of coherent light to go beyond the quantum backaction limit~\cite{Asjad2016Nov, Clark2017Jan}. In the model considered here, as depicted in Fig.~\ref{fig:Fano_and_squeezing}(a), this means replacing the vacuum noise input by squeezed noise, generated with a degenerate parametric amplifier (DPA)  (see Ref.~\cite{Gardiner2004}, chapter 10). This is an example of a cascaded open quantum system~\cite{Gardiner1993Apr}, where the output $\hat{c}_\text{out,s}$ of the DPA is used as an input for the cavity, namely $\hat{c}_\text{out,s}(t) = \hat{a}_\text{in,s}(t)$ (we neglect the delay due to propagation, which is irrelevant for our purposes). We consider that the DPA mode, with annihilation operator $\hat{c}$, is resonant with the cavity pump, i.e., at frequency $\omega_\text{las}$. The DPA pump is tuned with the parametric oscillator, i.e., at frequency $2\omega_\text{las}$.

The Langevin equation for the DPA reads~\cite{Asjad2016Nov}
\begin{equation}
    \dot{\hat c}=-\varkappa \hat{c} +\chi \hat{c}^\delta\hat{a}gger +\text{s}qrt{2 \varkappa}\hat c_\text{in},
\text{e}nd{equation}    
where $\hat{c}_\text{in} = \frac{1}{\text{s}qrt{\varkappa}}(\text{s}qrt{\varkappa_\text{s}}\hat{c}_\text{in,s} + \text{s}qrt{\varkappa_\text{v}}\hat{c}_\text{in,v})$ is the total input field, $\varkappa = \varkappa_\text{s} + \varkappa_\text{v}$ the corresponding total loss rate and $\chi =\abs{\chi}\text{e}^{-i2\theta}$ is the squeezing parameter with $\theta$ the squeezing angle [see Fig.~\ref{fig:Fano_and_squeezing}(a)]. The input vacuum noise $\hat{c}_\text{in,v}$ with associated coupling $\varkappa_\text{v}$ corresponds to the uncontrolled losses of the DPA. The desired output from the DPA is then finally obtained from $\hat{c}_\text{out,s}=\text{s}qrt{2\varkappa_\text{s}}\hat{c}-\hat{c}_\text{in,s}$. The Langevin equations for the optomechanical system are still Eqs.~\text{e}qref{standard:Langevin}, but now the total optical input noise reads $\hat{a}_\text{in} = \frac{1}{\text{s}qrt{\kappa}}(\text{s}qrt{\kappa_\text{s}}\hat{c}_\text{out,s} + \text{s}qrt{\kappa_\text{v}}\hat{a}_\text{in,v})$, where $\kappa_\text{s}$ is the coupling to squeezed vacuum noise and $\hat{a}_\text{in,v}$ corresponds to the uncontrolled losses of the cavity with associated coupling $\kappa_\text{v}$~\cite{Asjad2016Nov}. The total loss rate of the cavity is therefore $\kappa = \kappa_\text{s} + \kappa_\text{v}$. In the white-noise approximation, the correlation functions of the input noise read~\cite{Gardiner2004} (see Appendix~\ref{appendix:squeezing:noise})
\begin{subequations}\label{squeezing:correlations_a_in}
    \begin{align}
        \text{mec}ean{\hat{a}^{\delta\hat{a}gger}_\text{in}(t)\hat{a}_\text{in}(t')} &= \pi_\text{s} N_\text{s}\delta(t - t'),\\
        \text{mec}ean{\hat{a}_\text{in}(t)\hat{a}_\text{in}(t')} &= \pi_\text{s} M_\text{s} \delta(t - t').
    \text{e}nd{align}
\text{e}nd{subequations}  
The purity of the squeezing is $\pi_\text{s} = \varkappa_\text{s} \kappa_\text{s} / \varkappa \kappa$ with $\pi_\text{s} = 1$ for pure squeezing and $\pi_\text{s} = 0$ for vacuum noise only. Furthermore, we have  
\begin{subequations}\label{squeeezing:Ns and Ms}
    \begin{align}
    N_\text{s} &= \frac{4 r_\text{s}^2}{(1 - r_\text{s}^2)^2}, \label{squeezing:N_s}\\
    M_\text{s} &= \frac{2 r_\text{s}(1 + r_\text{s}^2)}{(1 - r_\text{s}^2)^2}\text{e}^{-i2\theta}. \label{squeezing:M_s}
\text{e}nd{align}
\text{e}nd{subequations}
with the squeezing ratio $r_\text{s} = \abs{\chi}/\varkappa$ and this ratio is such that $0 < r_\text{s} < 1$ so that the white-noise approximation holds. 

\text{s}ubsubsection{Fano-mirror setup}\label{sec:II:Fano-mirror setup}

An alternative strategy to improve the optomechanical cooling is to make the cavity density of states asymmetric in order to further suppress the detrimental Stokes process [see Fig.~\ref{fig:Fano_and_squeezing}(c)]. This can be achieved by replacing 
the left-hand side mirror of the cavity by a Fano mirror (e.g., a photonic crystal membrane)~\cite{Cernotik2019Jun}.
We then treat separately the cavity input noise $\hat{a}_\text{in,L}$ and $\hat{a}_\text{in,R}$ coming through the left-hand side and right-hand side mirrors, denoting $\kappa_\text{L}$ and $\kappa_\text{R}$ the corresponding loss rates. Both are vacuum noise and their correlation functions obey Eq.~\text{e}qref{standard:correlations_a_in}. The setup is depicted in Fig.~\ref{fig:Fano_and_squeezing}(b) and modeled by adding the internal mode of the Fano mirror to the optomechanical description. We consider an even Fano-mirror mode\footnote{This is the one studied in Ref.~\cite{Cernotik2019Jun}, but it could also be an odd mode, as treated in the Supplementary material of the same article.}, so the setup's Hamiltonian reads~\cite{Cernotik2019Jun}
\begin{equation}
    \hat{H}_\text{Fano} = \hat{H} + \hbar\omega_d \hat{d}^\delta\hat{a}gger \hat{d} + \hbar\lambda(\hat{a}^\delta\hat{a}gger \hat{d} + \hat{a} \hat{d}^\delta\hat{a}gger), \label{H_Fano}
\text{e}nd{equation}
where $\hat{H}$ is the Hamiltonian of the standard setup, given by Eq.~\text{e}qref{H_sys}. We have denoted $\hat{d}$ the annihilation operator for the considered mode of the Fano mirror, $\omega_d $ its frequency and $\lambda$ the coupling strength between this mode and the cavity mode.

The Langevin equations for such a system read~\cite{Cernotik2019Jun} (in the frame rotating at $\omega_\text{las}$)
\begin{subequations}\label{Fano:Langevin}
    \begin{align}
        \dot{\hat q}=\,&\Omega_\text{mec} \hat p, \\
        \dot{\hat p}=\,&-\Omega_\text{mec} \hat q-\gamma \hat p+g_0 \text{s}qrt{2}\hat{a}^\delta\hat{a}gger \hat{a}  +\text{s}qrt{\gamma } \hat\xi, \\
        \dot{\hat a}=\,&-(\kappa+i\Delta_0) \hat{a} +ig_0 \text{s}qrt{2}\hat{a} \hat q -i \text{e}psilon - \text{mec}athcal{G}\hat{d} \nonumber\\&+\text{s}qrt{2 \kappa_\text{R}}\hat{a}_\text{in,R} +\text{s}qrt{2 \kappa_\text{L}}\hat{a}_\text{in,L},\\
        \dot{\hat d}=\,&-(\gamma_d +i\Delta_d) \hat{d} - \text{mec}athcal{G}\hat{a}  +\text{s}qrt{2 \gamma_d }\hat{a}_\text{in,L},
    \text{e}nd{align}
\text{e}nd{subequations}
with $\gamma_d$ the coupling strength between the Fano-mirror mode and the left-hand side input, $\Delta_d = \omega_d-\omega_\text{las}$ the detuning with the laser drive and $\kappa = \kappa_\text{L} + \kappa_\text{R}$ the total loss rate of the cavity. We have also defined $\text{mec}athcal{G} = i\lambda + \text{s}qrt{\kappa_\text{L}\gamma_d }$. In the model considered here, the coupling between the Fano-mirror mode and the cavity mode is given by  $\lambda = \text{s}qrt{\kappa_0\gamma_d}$~\cite{Cernotik2019Jun}, where $\kappa_0$ is the loss rate that the left-hand side mirror would have if it were frequency independent. For simplicity we take $\kappa_0 = \kappa_\text{R}$.
 

For the Fano-mirror setup, the total loss rate of the cavity is very large (a lot larger than for the standard setup). Such a large loss rate would place a standard setup in the non-sideband resolved limit. However, this does not reflect the actual linewidth of the cavity. Instead, the presence of the Fano mirror leads to a different effective loss rate, corresponding to the \textit{effective} linewidth of the cavity, $\kappa_\text{eff}$ (see Sec.~\ref{sec:II:experimental setups:Fano}). We therefore use this effective linewidth in the cooperativity
 \begin{equation}
    C\rightarrow\frac{2g^2}{\kappa_\text{eff} \gamma \bar{n}_\text{mec}}\ ,\label{cooperativity_Fano}
\text{e}nd{equation}
 for this setup, in what follows.

\text{s}ubsection{Thermodynamic picture} \label{sec:II:Thermodynamic picture}

The thermodynamic system of interest is the mechanical resonator. The phonon number of this mechanical resonator fluctuates, due to its coupling to a mechanical environment, which constitutes a hot phonon bath at temperature $T_\text{mec}$. 
The goal of optomechanical (sideband) cooling is to reduce the amount of fluctuations in the resonator, thereby ``cooling it down". 
In order to do so, we have at our disposal a cold photon bath,   in the sense that it contains negligibly few thermal excitations.\footnote{Laser phase noise may be treated as a non-zero temperature bath, see \cite{rabl_phase-noise_2009}.}. However, heat can flow from the fluctuating mechanical resonator to this zero-temperature environment, only if a coupling mechanism between phonons and photons is established. This coupling is provided by the laser that is driving the cavity, thereby constituting a heat valve [see Fig.~\ref{fig:principle} for a sketch of this process]. Hence, the established heat flow from hot to cold comes with a cost, arising from coupling the system to the cold bath, i.e., from keeping the valve opened. 
Note that this \textit{cooling process} is very different from a typical thermodynamic \textit{refrigerator}, where one has the goal to make heat flow from a cold bath to a hot bath by providing work to the system. 

In order to characterize this cooling process, we analyze three thermodynamic quantities, which are complementary to each other and to a certain extent equivalent to those typically quantifying the performance of a thermodynamic cooling device. This is first of all the coldest temperature that can be reached, here represented by the lowest steady-state phonon number of the mechanical resonator. This quantity is experimentally accessible via calibrated optomechanical thermometry or sideband asymmetry determined from the cavity output light field. This is the quantity that is commonly considered when studying optomechanical cooling. Note, however, that the mechanical resonator is a microscopic system and does not constitute a bath itself. By its coupling to the hot mechanical phonon bath and the cold photon bath it is brought into a nonequilibrium state and the lowest reachable  steady-state phonon number can hence only be related to an \textit{effective} temperature. 

Importantly, in our work we complement this study by analyzing two further thermodynamic quantities. Typically, in a refrigerator the  cooling power is analyzed as a performance quantifier, namely the heat flow out of the bath to be cooled. As mentioned above, the mechanical resonator that is intended to be cooled is, however, \textit{not} a bath, and the total heat flow out of or into the resonator is always zero in the steady state. We therefore introduce instead the \textit{evacuated-heat flow} as an equivalent to the cooling power, namely the heat current, $J_\text{c}$, carried by phonons at frequency $\Omega_\text{mec}$ flowing out of the mechanical resonator and into the effective cold bath.
 We define $J_\text{c} = \hbar\Omega_\text{mec} I_\text{c}^\text{phonon}$ with the associated phonon flux, $I_\text{c}^\text{phonon}$. The latter quantity is experimentally accessible from the cavity output light field by measuring the amplitude and phase quadrature of the light field, see Appendix~\ref{appendix:standard:flows measurement}.

Finally, we analyze the efficiency of this cooling process.
As described above, maintaining the valve opened in a given position has a cost. The full cost is given by the laser power ${P}_\text{las}$ required to create the coherent field in the cavity mode. This motivates us to define a cooling efficiency
\begin{equation}
    \text{e}taL = \frac{J_\text{c}}{{P}_\text{las}}\ ,\label{etaL}
\text{e}nd{equation}
comparing the evacuated-heat flow (the desired output) with the full laser power (the exploited resource). We define an alternative cooling efficiency later in Eq.~(\ref{etaC}), which, at the cost of neglecting parts of the resource, gives a better picture of the conversion process from phonons to photons. 
This alternative efficiency [Eq.~(\ref{etaC})] can actually be more relevant for device realizations, when the laser light, which was not used in the conversion process, is partly re-used for further operations and, hence, is not lost.
Note that standard refrigerators or heat pumps are often characterized by ``coefficients of performance"; here we refrain from using this terminology and rather use the more generic notion ``(cooling) efficiency", since the desired output of the present process is not  the heat flow into or out of one of the baths, but rather the relevant part of the flux out of the microscopic system attached to them.

In the following, we want to consolidate the above described thermodynamic picture and therefore start with an analysis of the standard setup. As a first step, we linearize the Hamiltonian~\text{e}qref{H_sys} around the system's semiclassical steady state (see Appendix~\ref{appendix:standard:lin} for details)\footnote{This approximation is justified in most optomechanical setups as the optomechanical coupling $g_0$ is typically very small: $g_0 \ll \gamma(2\bar{n}_\text{mec} + 1), \kappa$~\cite{Aspelmeyer2014Dec, Bowen2015Nov}, see \Tref{tab:Params}.}, in which the cavity mode is in the coherent state of amplitude $\alpha$ and the resonator's rest position is $\bar{q}$~\cite{Genes2008Mar}, with
\begin{subequations}\label{standard:average values}
    \begin{align}
        \alpha&= \dfrac{-i\text{e}psilon}{\kappa+i \Delta}, \label{standard:alpha}\\
        \bar{q} &=\text{s}qrt{2}\dfrac{g_0}{\Omega_\text{mec}} \abs{\alpha}^2.
    \text{e}nd{align}
\text{e}nd{subequations} 
We define the effective detuning $\Delta=\Delta_0-g_0 \text{s}qrt{2}\bar{q}$ and choose the phase in $\text{e}psilon$ so that $\alpha=|\alpha|$ is real and positive. We split the operators into their average and fluctuation terms:  $\hat{a}=\alpha + \delta\hat{a}$, $\hat{q}=\bar{q} + \delta\hat{q}$, $\hat{p}=\delta\hat{p}$. The linearized Hamiltonian reads
\begin{align}\label{H_Lin}
    \hat{H}_\text{lin}=\,&\hbar {\Delta} \delta\hat{a}^\delta\hat{a}gger \delta\hat{a} +\hbar \Omega_\text{mec} \delta\hat{b}^\delta\hat{a}gger \delta\hat{b}\nonumber\\& -\hbar g (\delta\hat{a}+\delta\hat{a}^\delta\hat{a}gger  ) (\delta\hat{b} + \delta\hat{b}^\delta\hat{a}gger) ,
\text{e}nd{align}
where we have denoted $\delta\hat{b} = (\delta\hat{q} + i\delta\hat{p})\text{s}qrt{2}$. Indeed, we see from the expression of the effective coupling $g = g_0 \abs{\alpha}$ that the light field in the cavity plays a key role in the cooling process, acting like a valve that allows phonons to flow towards the cavity.

For the sake of obtaining a simple and more usual thermodynamics picture, we---in this section---furthermore resort to the weak-coupling regime ($g \ll \Omega_\text{mec}, \kappa$).  Tracing out the cavity degrees of freedom, we then find the  master equation~\cite{Wilson-Rae2008Sep} describing the dynamics of the mechanical resonator,
\begin{align}\label{weak-coupling:master_eq}
    \dv{\hat{\rho}_\text{mec}}{t}=\,&-i [\Omega_\text{mec}^\text{eff}\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}, \hat{\rho}_\text{mec}]+ (\gamma \bar{n}_\text{mec} + A_+) D[\delta\hat{b}^\delta\hat{a}gger]\hat{\rho}_\text{mec} \nonumber\\ & + (\gamma (\bar{n}_\text{mec} + 1) + A_-)D[\delta\hat{b}]\hat{\rho}_\text{mec}\ .
\text{e}nd{align}
Here, $\hat{\rho}_\text{mec}$ is the density operator of the resonator, $\Omega_\text{mec}^\text{eff}$ the effective mechanical frequency (see Appendix~\ref{appendix:standard:solution Langevin}) and $D[\hat{O}]\hat{\rho} = \hat{O}\hat{\rho} \hat{O}^\delta\hat{a}gger - \frac{1}{2} \{\hat{O}^\delta\hat{a}gger \hat{O}, \hat{\rho}\}$. We have defined the rates
\begin{equation}
    A_\pm = \frac{2g^2\kappa}{\kappa^2 + (\Omega_\text{mec} \pm {\Delta})^2}\label{Stokes_rates_perturbative},
\text{e}nd{equation}
which correspond to the Stokes and anti-Stokes processes, respectively.  This is the situation depicted in Fig.~\ref{fig:weak coupling cooling}. 
Equivalently, Eq.~\text{e}qref{weak-coupling:master_eq} can be rewritten in the form
\begin{align}
    \!\!\dv{\hat{\rho}_\text{mec}}{t}=\,&-i [\Omega_\text{mec}^\text{eff}\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}, \hat{\rho}_\text{mec}]\nonumber\\
    &+ \Gamma_\text{opt}\left( \bar{n}_\text{opt} D[\delta\hat{b}^\delta\hat{a}gger] + (\bar{n}_\text{opt} + 1)D[\delta\hat{b}]\right)\hat{\rho}_\text{mec}\quad\quad\,\,\,\nonumber\\
    &+ \gamma\left( \bar{n}_\text{mec} D[\delta\hat{b}^\delta\hat{a}gger] + (\bar{n}_\text{mec} + 1)D[\delta\hat{b}]\right)\hat{\rho}_\text{mec}.
\text{e}nd{align}
This equation shows that the driven cavity acts on the resonator like an effective thermal \textit{phonon} bath of associated damping rate $\Gamma_\text{opt} \text{e}quiv A_- - A_+$ and containing a finite number of phonons $\bar{n}_\text{opt} = A_+/\Gamma_\text{opt}$. Hence, the cavity and its environment play the role of an effective cold bath (dashed box in Fig.~\ref{fig:principle}(b)), of finite effective temperature $T_\text{opt}$ such that $\bar{n}_\text{opt} = (\text{e}^{\hbar\Omega_\text{mec}/k_\text{B} T_\text{opt}} - 1)^{-1}$. As a result, the average phonon number associated to the fluctuations of the resonator in the steady state $\bar{n}_\text{fin} = \text{mec}ean{\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}}_\text{ss}$ can be interpreted as the steady state of a system coupled to two thermal baths~\cite{Aspelmeyer2014Dec},
\begin{equation}
    \bar{n}_\text{fin} = \frac{\gamma\bar{n}_\text{mec} + \Gamma_\text{opt}\bar{n}_\text{opt}}{\gamma + \Gamma_\text{opt}}.\label{weak-coupling:neff}
\text{e}nd{equation}
We can see from Eq.~\text{e}qref{Stokes_rates_perturbative} that $\Gamma_\text{opt}$ is proportional to $g^2 = g_0^2 |\alpha|^2$. This confirms that the light field in the cavity is like a heat valve allowing more or less heat to flow between the resonator and the cold phonon bath. Note that the effective cold bath is an \textit{engineered bath} and therefore it comes with a cost. This further motivates the above described thermodynamics picture and the resulting definition for the cooling efficiency as given in Eq.~\text{e}qref{etaL}.

In the weak-coupling regime, the evacuated-heat flow is explicitly given by
\begin{align}\label{weak-coupling:Jc}
    J_\text{c} &= \hbar\Omega_\text{mec}\Gamma_\text{opt}(\bar{n}_\text{fin} - \bar{n}_\text{opt}) \nonumber\\
    &= \hbar\Omega_\text{mec}\frac{\gamma\Gamma_\text{opt}}{\gamma + \Gamma_\text{opt}}(\bar{n}_\text{mec} - \bar{n}_\text{opt}).
\text{e}nd{align}
Using Eq.~\text{e}qref{standard:average values}, the full efficiency, Eq.~\text{e}qref{etaL}, becomes
\begin{equation}
    \text{e}taL = \frac{\Omega_\text{mec} 2\kappa}{\omega_\text{las}\alpha^2(\kappa^2 + \Delta^2)}\frac{\gamma\Gamma_\text{opt}}{\gamma + \Gamma_\text{opt}}(\bar{n}_\text{mec} - \bar{n}_\text{opt})\ .
\text{e}nd{equation}
This weak-coupling thermodynamic picture also gives us further insights into the improvements provided by alternative setups. Squeezing the input noise, as described in Sec.~\ref{sec:II:squeezed-light setup}, allows to reduce the position quadrature noise of the intra-cavity field, which couples to the mechanics. As the cavity transfer function converts part of the input momentum noise into position noise and conversely, it is necessary to rotate the input squeezed state by the adequate angle $\theta$ (given in Sec.~\ref{sec:III:neff:squeezing}) so as to compensate for this effect. This reduction of the intra-cavity position noise can be interpreted as lowering the temperature of the effective cold phonon bath, i.e., $\bar{n}_\text{opt}$ is smaller while $\Gamma_\text{opt}$ remains the same as in the standard setup. On the other hand, with the Fano-mirror setup from Sec.~\ref{sec:II:Fano-mirror setup}, the imbalance between rates of the Stokes and anti-Stokes processes is increased, which results in a change in $\Gamma_\text{opt}$ and a decrease in $\bar{n}_\text{opt}$. Therefore, the effective cold bath gets colder and its coupling to the resonator is modified.

\text{s}ubsection{Examples for specific realizations}\label{sec:II:experimental setups}

We now turn to inspect specific optomechanical realizations in order to illustrate our thermodynamic analysis of optomechanical cooling with concrete examples and also to benchmark these realizations against each other. With the evacuated-heat flow and the cooling efficiency, our analysis yields complementary insights beyond the achievable phonon occupation as a benchmark only. To this end, we select four representative optomechanical systems, whose parameters are summarized in \Tref{tab:Params} in the appendix. From here on, we refer to them as system (1) to system (4).

System (1) employs mechanical resonators in the MHz regime and is placed in the sideband-resolved regime. Such systems can be realized, for example, with free-space Fabry-Pérot-based optomechanics~\cite{thompson_strong_2008,groblacher_demonstration_2009}. Importantly, ground-state cooling has already been demonstrated with such a system in the membrane-in-the-middle configuration~\cite{underwood_measurement_2015}. An advantage of the Fabry-Pérot system is its versatility as it can be directly modified to realize a Fano-mirror setup~\cite{Cernotik2019Jun,fitzgerald_cavity_2021}. System (2) is an integrated optomechanics device made of an optomechanical crystal with mechanical motion in the GHz regime that is placed in the sideband-resolved regime. This system is a prominent implementation for realizing optomechanics-based quantum networks with demonstrations of ground-state cooling~\cite{Chan2011} and non-classical state generation~\cite{riedinger_non-classical_2016,riedinger_remote_2018}. Our thermodynamic analysis of this system yields complementary insights beyond its use in quantum networks. System (3) is an optically levitated nanoparticle coupled to an optical cavity. Recently, ground-state cooling of the center-of-mass motion of such a levitated particle has been demonstrated~\cite{Delic2020}. Levitated optomechanics constitutes a pertinent platform for exploring thermodynamic heat engines or stochastic and quantum thermodynamics phenomena. Multiple experiments are already performed along these lines~\cite{gieseler_levitated_2018}. System (4) is based on the membrane-in-the-middle configuration~\cite{thompson_strong_2008}, but in the non-sideband-resolved regime. Recently, feedback-based cooling to the ground state~\cite{Rossi2018} has been demonstrated in such a system. We here show that one can actually achieve ground-state cooling when employing a squeezed-light drive instead of a coherent one, even in the absence of feedback, as evidenced in Refs.~\cite{Asjad2016Nov, Clark2017Jan, Lau2020Mar}. 

We analyze the considered four systems in the context of three different setups: (i) the standard setup, (ii) the squeezed-light setup and (iii) the Fano-mirror setup. The corresponding choice of parameters is shortly discussed in the following.

\text{s}ubsubsection{Standard setup}

For the standard setup, we assume driving of the optomechanical system with coherent laser light and use the parameters of the four systems as given in \Tref{tab:Params}. In this way, we basically follow the experimental implementations of Refs.\cite{Eichenfield2009Nov, Chan2011,Delic2020,Delic2020_1,Rossi2018}.

\text{s}ubsubsection{Squeezed-light setup}

In principle, the four systems can also be driven with a squeezed-light source. However, this has so far only been done and used for demonstrating ground-state cooling in microwave optomechanics~\cite{Clark2017Jan}. Squeezed states of light are customarily generated using non-linear optics processes. For example, squeezing levels of 12.7\,dB~\cite{eberle_quantum_2010} and 13\,dB~\cite{schonbeck_13db_2018} at wavelengths of 1064\,nm and 1550\,nm, respectively, have been generated. These wavelengths are also typically used in optomechanics experiments, see \Tref{tab:Params}. Squeezing in these experiments \cite{eberle_quantum_2010,schonbeck_13db_2018} is customarily observed at MHz sideband frequencies, which matches optomechanics implementations with mechanical motion in the MHz regime. Squeezing at GHz sideband frequencies, which is required when employing mechanical motion at GHz frequencies, has been generated at 5\,dB levels~\cite{ast_high-bandwidth_2013,kashiwazaki_continuous-wave_2020,tasker_silicon_2021}. In our analysis, we assume a squeezed-light source with realistic squeezing levels\footnote{The amount of squeezing in dB is calculated as $-10\log_{10}(2\text{mec}ean{\Delta P^\theta_\text{in,s}})$ with the variance of shot noise being $\frac{1}{2}$ and $\text{mec}ean{\Delta P^\theta_\text{in,s}}$ the variance of the squeezed quadrature.} of 0.87\,dB, 0.59\,dB, 2.7\,dB and 15.4\,dB for the four systems, respectively. 

\text{s}ubsubsection{Fano-mirror setup}\label{sec:II:experimental setups:Fano}

Free-space optomechanical devices can be directly adapted to realize a Fano-mirror setup, see, e.g., Refs.~\cite{Cernotik2019Jun,fitzgerald_cavity_2021}. To be able to compare the Fano-mirror setup with the standard setup, we use the fact that the effective linewidth of the cavity with the Fano mirror is $\kappa_\text{eff} = \gamma_d/\zeta_0^2$~\cite{Cernotik2019Jun}, where $\zeta_0$ is the polarizability of the left mirror at frequency $\omega_d$. Since we have chosen $\kappa_\text{R} = \kappa_0 $, $\zeta_0$ is also the (frequency-independent) polarizability of the right mirror. Denoting $\Gamma$ the free spectral range of the cavity that is  given in \tref{tab:Params}, we have $\kappa_\text{L} = 2\Gamma$ and $\kappa_0 = \Gamma/2\zeta_0^2$.
Moreover, choosing to use the optimal value $\gamma_d^* = 4\Omega_\text{mec}\zeta_0$~\cite{Cernotik2019Jun}, we can determine $\zeta_0$ for each system  such that $\kappa_\text{eff}$ is equal to the system's cavity loss rate. 

Optomechanical cavities with Fano mirrors can be realized by using suspended photonic crystal slabs, see for example Refs.~\cite{bui_high-reflectivity_2012,norte_mechanical_2016,chen_high-finesse_2017,gartner_integrated_2018,kini_manjeshwar_suspended_2020}. The parameters of the photonic crystal slab~\cite{fan_analysis_2002} allow engineering of a desired reflectivity at frequency $\omega_d$ and, thus, a corresponding $\kappa_\text{eff}$. Fine tuning of the optical resonance of the photonic crystal mirror can be achieved in post-fabrication, for instance by using thermal tuning~\cite{Peter2005Aug},  strain tuning~\cite{Wong2004Feb} or by depositing~\cite{Yang2007Oct} or etching~\cite{bernard_precision_2016} atomic layers of the suspended photonic crystal. Therefore it is reasonable to assume that $\gamma_d$ close to the optimal value can be achieved in an experimental setup.

\text{s}ection{Thermodynamic performance}\label{sec:III}

In \text{s}ref{sec:II} we have described the cooling process in the three setups and have outlined how to set up a consistent thermodynamic picture. In this section, we analyze the thermodynamic performance of four example realizations and compare their performance for the operation in the three different setups. The parameters used for the four example realizations are given in \tref{tab:Params}. The analysis of the thermodynamic performance includes the lowest reachable phonon number, the evacuated-heat flow as well as the cooling efficiency. In order to obtain all these quantities, we first need to analyze the dynamics of the systems in the three setups.

\text{s}ubsection{Dynamics}

The first step to the study of the systems' dynamics in each of the three setups is to linearize the Langevin equations, as we have previously done in Sec.~\ref{sec:II:standard setup} for the standard setup. Here we present a short overview for all three setups. We then derive the evolution of the second-order moments, and in particular of the photon and phonon numbers.

\text{s}ubsubsection{Standard setup}\label{sec:III:dynamics:standard}

Applying the linearization introduced in Sec.~\ref{sec:II:Thermodynamic picture} to the Langevin equations for the standard setup, Eq.~\text{e}qref{standard:Langevin}, we obtain
\begin{subequations}\label{standard:Langevin lin}
    \begin{align}
        \delta \dot{\hat{q}}&=\Omega_\text{mec}  \delta\hat{p}, \\
        \delta\dot{\hat{p}}&=-\Omega_\text{mec}  \delta\hat{q}-\gamma  \delta\hat{p} +g\text{s}qrt{2}  (\delta\hat{a}^\delta\hat{a}gger + \delta\hat{a})  +\text{s}qrt{\gamma } \hat{\xi}, \\
        \delta \dot{\hat{a}}&=-(\kappa+i\Delta )  \delta\hat{a} +ig\text{s}qrt{2}  \delta\hat{q} +\text{s}qrt{2 \kappa} \hat{a}_\text{in}.
    \text{e}nd{align}
\text{e}nd{subequations}
From these equations and the noise correlation functions, given by Eqs.~\text{e}qref{correlations_xi} and~\text{e}qref{standard:correlations_a_in}, we derive the evolution of the photon and phonon numbers associated to the fluctuations\footnote{Note that the environment contribution in Eq.~\text{e}qref{standard:evolution phonon number} has an extra term in $\text{mec}ean{\delta\hat{b}^2 + \delta\hat{b}^{\delta\hat{a}gger 2}}$ compared to the one in Ref.~\cite{Wilson-Rae2008Sep} because Ref.~\cite{Wilson-Rae2008Sep} used a different noise model. Both models give very similar results in the high mechanical quality factor regime considered here. See Appendix~\ref{appendix:standard:approxs} for a discussion of the differences between the two approaches.} 
\begin{subequations}\label{standard:evolution photon and phonon numbers}
    \begin{align}
    \!\dv{\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}}}{t} =\,
    & i g \left(\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger (\delta\hat{b}^\delta\hat{a}gger +  \delta\hat{b})} - \text{mec}ean{\delta\hat{a} (\delta\hat{b}^\delta\hat{a}gger+ \delta\hat{b})}\right)\nonumber\\
    & - 2\kappa\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}}, \label{standard:evolution photon number}\\ 
    \!\dv{ \text{mec}ean{\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}}}{t} =\,
    & ig\left(\text{mec}ean{( \delta\hat{a} +  \delta\hat{a}^\delta\hat{a}gger)\delta\hat{b}^\delta\hat{a}gger} - \text{mec}ean{(\delta\hat{a}^\delta\hat{a}gger+\delta\hat{a})\delta\hat{b} }\right)\nonumber\\
    &\!\!\!+ \gamma \left(\bar{n}_\text{mec} - \text{mec}ean{\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}}  + \frac{\text{mec}ean{\delta\hat{b}^2 + \delta\hat{b}^{\delta\hat{a}gger 2}}}{2}\right).\!\label{standard:evolution phonon number}
    \text{e}nd{align}
\text{e}nd{subequations}
Details about the derivation and the evolution equations of the other second-order moments can be found in Appendix~\ref{appendix:standard:Lyapunov}.
For each of the above equations, two contributions stand out: the optomechanical one, proportional to $g$, and the one from the optical and mechanical environment, proportional to $\kappa$ and $\gamma$, respectively. When one neglects the non-resonant processes $\delta\hat{a}\delta\hat{b}$ and $\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger$ (rotating wave approximation, valid for small coupling $g$, see Appendix~\ref{appendix:standard:approxs}), the optomechanical contribution corresponds to the scattering picture from Fig.~\ref{fig:weak coupling cooling}(b), with heating from the Stokes process and cooling from the anti-Stokes process.

\text{s}ubsubsection{Squeezed-light setup}

In the case of the squeezed-light setup, described in Sec.~\ref{sec:II:squeezed-light setup}, the DPA and cavity are only coupled through the noise. Therefore, the average fields in the cavity and resonator are still given by Eqs.~\text{e}qref{standard:average values} and the linearized Langevin equations are the same as Eqs.~\text{e}qref{standard:Langevin lin}. In contrast, the cavity input noise correlation functions are given by Eqs.~\text{e}qref{squeezing:correlations_a_in}. As a consequence, the photon number evolution becomes
\begin{align}
    \dv{\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}}}{t} =\,
    & i g \left(\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger (\delta\hat{b}^\delta\hat{a}gger +  \delta\hat{b})} - \text{mec}ean{\delta\hat{a} (\delta\hat{b}^\delta\hat{a}gger+ \delta\hat{b})}\right)\nonumber\\
    & +2\kappa(\pi_\text{s} N_\text{s} - \text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}} )\label{squeezing:evolution photon number},
\text{e}nd{align}
while the evolution equation of the phonon number is unchanged (Eq.~\text{e}qref{standard:evolution phonon number}).
Further changes to the evolution of the second-order moments are detailed in Appendix~\ref{appendix:squeezing:Lyapunov}.

\text{s}ubsubsection{Fano-mirror setup}

Like in Sec.~\ref{sec:II:Thermodynamic picture}, we linearize the Langevin equations~\text{e}qref{Fano:Langevin} by splitting operators into their average terms and fluctuation terms. But in this case, we have one extra operator, $\hat{d}=\delta + \delta\hat{d}$, which modifies the semiclassical steady state into (see Appendix~\ref{appendix:Fano:lin} for details)
\begin{subequations}\label{Fano:average values}
    \begin{align}
    \alpha&=\frac{-i\text{e}psilon}{\kappa+i \Delta - \text{mec}athcal{G}^2/(\gamma_d  + i\Delta_d )}\label{Fano:alpha} \\
    \bar{q} &=\text{s}qrt{2}\dfrac{g_0}{\Omega_\text{mec}} \abs{\alpha}^2,\\
    \delta&=\frac{-\text{mec}athcal{G}}{\gamma_d  + i\Delta_d } \alpha.
\text{e}nd{align}
\text{e}nd{subequations}
The linearized Langevin equations read
\begin{subequations}\label{Fano:Langevin lin}
    \begin{align}
        \delta \dot{\hat{q}}=\,&\Omega_\text{mec}  \delta\hat{p}, \\
        \delta\dot{\hat{p}}=\,&-\Omega_\text{mec}  \delta\hat{q}-\gamma  \delta\hat{p} +g\text{s}qrt{2}  (\delta\hat{a}^\delta\hat{a}gger + \delta\hat{a})  +\text{s}qrt{\gamma } \hat{\xi}, \\
        \delta \dot{\hat{a}}=\,&-(\kappa+i\Delta )  \delta\hat{a} +ig\text{s}qrt{2}  \delta\hat{q} -\text{mec}athcal{G}\delta\hat{d} \nonumber\\ &+\text{s}qrt{2 \kappa_\text{L}} \hat{a}_\text{in,L}+\text{s}qrt{2 \kappa_0} \hat{a}_\text{in,R},\\
        \delta \dot{\hat{d}}=\,&-(\gamma_d +i\Delta_d  )  \delta\hat{d} -\text{mec}athcal{G}\delta\hat{a} +\text{s}qrt{2 \gamma_d } \hat{a}_\text{in,L}.
    \text{e}nd{align}
\text{e}nd{subequations}
The free parameter $\Delta_d$ of the Fano-mirror mode was found to be constrained by the cavity mode and the optical loss rates by  $\Delta_d = \Delta - 2\text{s}qrt{\kappa_0\kappa_\text{L}}$~\cite{Cernotik2019Jun}.

From these equations and the noise correlation functions, given by Eqs.~\text{e}qref{correlations_xi} and \text{e}qref{standard:correlations_a_in}, we derive the evolution of the photon number (see Appendix~\ref{appendix:Fano:Lyapunov}):
\begin{align}
    \dv{\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}}}{t} =\,
    & i g \left(\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger (\delta\hat{b}^\delta\hat{a}gger +  \delta\hat{b})} - \text{mec}ean{\delta\hat{a} (\delta\hat{b}^\delta\hat{a}gger+ \delta\hat{b})}\right) \nonumber\\
    & -2 \text{R}e{\text{mec}athcal{G}\text{mec}ean{\delta\hat{d} \delta\hat{a}^\delta\hat{a}gger}}  -2\kappa\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}} \label{Fano:evolution photon number},
\text{e}nd{align}
The evolution equation of the phonon number is unchanged (Eq.~\text{e}qref{standard:evolution phonon number}).

\text{s}ubsection{Steady-state phonon number}

The first performance indicator we consider is the steady-state phonon number $\bar{n}_\text{fin} = \text{mec}ean{\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}}_\text{ss}$, which is the standard quantity that is used to characterize sideband cooling and experimentally accessible via calibrated noise thermometry or sideband asymmetry. In a thermodynamic picture, this quantity would correspond to the lowest temperature that can be reached by a refrigerator.   
Indeed, the average phonon number in the fluctuations can be related to $T_\text{eff}$, the \textit{effective} temperature of the resonator, by
\begin{equation}
   \bar{n}_\text{fin} = \left(\text{e}^{\hbar\Omega_\text{mec}/k_\text{B} T_\text{eff}} - 1\right)^{-1}.
\text{e}nd{equation}
 Since the position of the resonator has been shifted (see Sec.~\ref{sec:II:Thermodynamic picture}), we actually have the mean stationary phonon number $\text{mec}ean{\hat{b}^\delta\hat{a}gger \hat{b}}_\text{ss} = (g/\Omega_\text{mec})^2|\alpha|^2 + \bar{n}_\text{fin}$. 
 
 In the following, we describe how to obtain the steady-state phonon number, $\bar{n}_\text{fin}$, in the three setups and discuss some of its general properties. We then discuss differences between the different setups and example realizations in the end of this section.
 
 

\text{s}ubsubsection{Standard setup} 

For the standard setup, $\bar{n}_\text{fin}$ can be found analytically, by solving the equations for the second-order moments for the steady state. However, the obtained expression is rather complex and we therefore present it in Appendix~\ref{appendix:standard:Lyapunov}, see also Refs.~\cite{Genes2008Mar,Wilson-Rae2008Sep}. Instead, in order to illustrate our results, we here present plots, see Fig.~\ref{fig:neff}, using the realistic experimental parameters from Tab.~\ref{tab:Params}. The only free parameter in the standard setup is the detuning $\Delta$. In the resolved-sideband regime, $\kappa \ll \Omega_\text{mec}$, the lowest phonon number is obtained for $\Delta \text{s}imeq \Omega_\text{mec}$, as illustrated in Fig.~\ref{fig:weak coupling cooling}. Note however, that the optimal detuning is exactly given by the mechanical frequency only in the weak-coupling approximation. All the plots for the standard setup are done at $\Delta = \Omega_\text{mec}$, except for the experimental system (4) (see Tab.~\ref{tab:Params}), where we use $\Delta = \kappa$. System (4) is not in the resolved-sideband regime and the cavity linewidth is larger than the mechanical frequency. As a consequence, if $\Delta = \Omega_\text{mec}$, the Stokes and anti-Stokes rates are almost equal, not leading to effective cooling. The best cooling is obtained when the difference between the two rates is the largest, that is on the slope of the cavity density of states, for a detuning close to $\kappa$. 

\text{s}ubsubsection{Squeezed-light setup} \label{sec:III:neff:squeezing}

As in the standard case, we obtain an analytical expression of the steady-state phonon number (see Appendix~\ref{appendix:squeezing:Lyapunov}) for the squeezed-light setup.
Again, the lowest phonon number is obtained for $\Delta \text{s}imeq \Omega_\text{mec}$. However, we have three additional parameters: the squeezing purity $\pi_\text{s}$, the squeezing ratio $r_\text{s}$ and the squeezing angle $\theta$ (see Sec.~\ref{sec:II:squeezed-light setup}). As expected, the best results with respect to the lowest reachable phonon number are obtained for pure squeezing, i.e. $\pi_\text{s} = 1$. Besides, there exists an optimal value of the squeezing angle $\theta^*$ that minimizes the phonon number.
Interestingly, this angle does not depend on the other squeezing parameters (see Eq.~\text{e}qref{squeezing:optimal angle} in Appendix~\ref{appendix:squeezing:optimization} for the full expression) and, in the limit of weak coupling and neglecting the mechanical damping rate, it is given by~\cite{Asjad2016Nov}
\begin{equation}
    \tan(2\theta^*) = \frac{2  {\Delta} {\kappa}}{{\Delta}^{2} - {\Omega_\text{mec}}^{2} - {\kappa}^{2}}.
\text{e}nd{equation}
Therefore it is possible to determine this optimal angle from the optomechanical setup's parameters. There exists also an optimal value $r^*_\text{s}$ of the squeezing ratio which depends on the squeezing angle (see Appendix~\ref{appendix:squeezing:optimization} for details).
All the plots for this setup are done using the parameters from Tab.~\ref{tab:Params}, with $\Delta = \Omega_\text{mec}$ (except for system (4) for which $\Delta = \kappa$), $\pi_\text{s} = 1$ and the  values of $\theta$ and $r_\text{s}$, which are optimal for reducing the steady-state phonon number.

Going to the limit of weak coupling and neglecting the mechanical damping rate, we derive the phonon number in the effective cold bath (see Sec.~\ref{sec:II:Thermodynamic picture}) 
\begin{eqnarray}
    \bar{n}_\text{opt}^\text{squ} & = & \bar{n}_\text{opt}\left( 1 + \frac{2\pi_\text{s} N_\text{s}(\Delta^2 + \Omega_\text{mec}^2 + \kappa^2)}{\kappa^2 + (\Delta - \Omega_\text{mec})^2} \right.\\
&&   \left. -2 \pi_\text{s} \abs{M_\text{s}}\text{s}qrt{\frac{\kappa^2 + (\Delta + \Omega_\text{mec})^2}{\kappa^2 + (\Delta - \Omega_\text{mec})^2}}\cos(2\theta - 2\theta^*) \right),\nonumber
\text{e}nd{eqnarray} 
where $\bar{n}_\text{opt}$ is the phonon number in the effective cold bath in the absence of squeezing, in agreement with the results from~\cite{Asjad2016Nov}. This explicitly shows that, with the appropriate squeezing parameters, we can get a \textit{colder} effective cold phonon bath than with the standard setup, thus allowing to reach a lower $\bar{n}_\text{fin}$. In this limit, the optimal squeezing ratio $r_\text{s}^*$ is such that 
\begin{equation}
    \frac{\abs{M_\text{s}(r_\text{s}^*)}}{N_\text{s}(r_\text{s}^*)} = \text{s}qrt{\frac{\kappa^2 + (\Delta + \Omega_\text{mec})^2}{\kappa^2 + (\Delta - \Omega_\text{mec})^2}}
\text{e}nd{equation}
and then $\bar{n}_\text{opt}^\text{squ} = \bar{n}_\text{opt}(1 - \pi_\text{s})$~\cite{Asjad2016Nov}. Therefore, in this ideal case and with perfect squeezing, the lowest attainable temperature of the effective cold phonon bath is zero.

\text{s}ubsubsection{Fano-mirror setup} 

In this setup, the cavity density of states is centered around $\omega_d$ \cite{Cernotik2019Jun} [see Fig.~\ref{fig:Fano_and_squeezing}(c)], therefore the lowest phonon number can be obtained for $\Delta_d \text{s}imeq \Omega_\text{mec}$~\cite{Cernotik2019Jun}. As explained in Sec.~\ref{sec:II:experimental setups:Fano}, Fano mirrors with the desired loss rate $\gamma_d$ can be fabricated by engineering the photonic crystal lattice parameters. It is therefore realistic to use the loss rate $\gamma_d^*$, for which the phonon number gets minimized, see Tab.~\ref{tab:Params} for the numerical values, and $\Delta_d = \Omega_\text{mec}$ in the plots. Once again, to obtain the best cooling, we use a different detuning for system (4): $\Delta_d = 2.5\Omega_\text{mec}$ (see Appendix \ref{appendix:Fano:non-sideband-resolved}).
We solved the equations for the second-order moments in order to obtain $\bar{n}_\text{fin}$  
numerically (see Appendix~\ref{appendix:Fano:Lyapunov}).

\text{s}ubsubsection{Discussion}

\begin{figure}[tb!]
    \includegraphics[width=\linewidth]{fig_neff.pdf}
    \caption{\label{fig:neff}
        Steady-state phonon number as a function of cooperativity for the three setups (indicated by different lines in each plot). The dotted and dashed black lines respectively indicate the average phonon number in the hot mechanical bath and in the effective cold bath (weak-coupling limit) for the standard setup. The vertical solid lines indicate the cooperativities at which the standard and squeezed-light setups (in blue) and Fano-mirror setup (in green) become unstable (see Appendices \ref{appendix:standard:lin} and \ref{appendix:Fano:lin}).The number in each sub-figure corresponds to the system number as given in Tab.~\ref{tab:Params}.
    }
\text{e}nd{figure}

We have plotted $\bar{n}_\text{fin}$ as a function of cooperativity for each setup and for each system from Tab.~\ref{tab:Params} in Fig.~\ref{fig:neff}.  For all four example systems in all three cooling setups, we notice  that at low cooperativity, the thermal noise prevails over the optomechanical coupling and $\bar{n}_\text{fin} \text{s}imeq \bar{n}_\text{mec}$. Then, when the cooperativity increases, the optomechanical coupling allows to reach lower phonon numbers, in most cases down to the ground state of the mechanical resonator (namely $\bar{n}_\text{fin}<1$). The minimum is  reached, when the optomechanical coupling becomes so large that $\bar{n}_\text{fin}$ is governed by the occupation of the effective cold phonon bath, namely $\bar{n}_\text{fin}\approx \bar{n}_\text{opt}$ (see Eq.~\text{e}qref{weak-coupling:neff} obtained in the weak-coupling limit). Finally, when further increasing the cooperativity, the coupling becomes so strong that the optomechanical back-action increases the fluctuations of the resonator and the phonon number increases as well\footnote{In the strong coupling limit, lower phonon numbers can be obtained in an instantaneous instead of the steady state, see Ref.~\cite{Liu2014May}.}. This is the case when non-resonant processes of the type $\delta\hat{a}\delta\hat{b}$ and $\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger$ start to play a role. This behavior due to back-action is not captured by the weak-coupling approximation (Eq.~\text{e}qref{weak-coupling:neff}, see also Appendix~\ref{appendix:standard:approxs}). Eventually, the phonon number diverges\footnote{
The system is no longer stable, see Appendices \ref{appendix:standard:lin} and \ref{appendix:standard:Lyapunov}.}. Interestingly, Fig.~\ref{fig:quadratures} shows that the divergence of $\bar{n}_\text{fin} = (\text{mec}ean{\delta\hat{q}^2} + \text{mec}ean{\delta\hat{p}^2}  -1)/2$ fully comes from the position fluctuations while the momentum fluctuations keep decreasing with cooperativity.

There are, however, also differences between the different systems and between the different setups. 
Comparing the different setups, we find that both the driving with squeezed light as well as the use of a Fano mirror allow to reach lower phonon numbers than the standard setup, as expected. Interestingly, the steady-state phonon number is however only little decreased with respect to the standard setup. For systems (1) to (3), the difference is hardly visible on the log-scale used in Fig.~\ref{fig:neff}; while the absolute difference in the phonon number is indeed small, 
the relative difference with respect to the standard setup is quite significant, as visible in Tab.~\ref{tab:rel diff x}.

\begin{table}[tb]
    \text{s}etlength{\tabcolsep}{10pt}
    \begin{tabular}{c|cccc}
        \hline
         System&  (1) &  (2) &  (3) &  (4) \\ \hline\hline
        $x^\text{squeezed}$ & 14.4\% & 22.5\% & 54.8\% &99.4\% \\\hline
        $x^\text{Fano}$ &  5.0\%& 22.7\%& 56.5\% & 66.4\%\\ \hline
    \text{e}nd{tabular}
    \caption{\label{tab:rel diff x}
       Relative difference between the minimum phonon number $\bar{n}_\text{fin}$ reached with the standard and the squeezed-light/Fano-mirror setup, $x^\text{squeezed/Fano}=\left(\bar{n}_\text{fin}^\text{standard}-\bar{n}_\text{fin}^\text{squeezed/Fano}\right)/\bar{n}_\text{fin}^\text{standard}$.
    }
\text{e}nd{table}

\begin{figure}[b!]
    \includegraphics[width=\linewidth]{fig_quadratures.pdf}
    \caption{\label{fig:quadratures}
        Fluctuations of the mechanical position (solid lines) and momentum (dashed lines) quadratures as a function of cooperativity for the three setups: standard setup in blue (dark gray), squeezed-light setup in red and orange (medium gray) and Fano-mirror setup in green (light gray). The number in each sub-figure corresponds to the system number as given in Tab.~\ref{tab:Params}. For systems (1) to (3) the curves for the standard and squeezed-light setups are almost identical.
    }
\text{e}nd{figure}

Furthermore, the dips in the curves for the Fano-mirror setup are always shifted to larger cooperativity values, which means that a larger optomechanical coupling is required to cool the resonator. This can be reached by increasing the power of the laser driving. The shift can be understood by looking at the weak-coupling picture. For the systems we consider, the coupling $\Gamma_\text{opt}$ to the effective cold phonon bath, at a given $C$, is weaker in the Fano-mirror setup than in the standard setup (see Appendix~\ref{appendix:Fano:weak-coupling}).

For all systems in the resolved-sideband regime (the first three in Tab.~\ref{tab:Params}), all setups enable ground-state cooling. While the overall behavior seems very similar for systems (1) to (3), the initial phonon number has very different orders of magnitude in the different systems, meaning that also the overall reduction in order to reach the ground state is not the same. In contrast system (4), which is not in the resolved-sideband regime, does not allow for ground state cooling when using a coherent drive\footnote{This would require feedback cooling (see for instance~\cite{Genes2008Mar}), as was done in the experiment using this setup~\cite{Rossi2018}.}. However, driving with squeezed light  reduces $\bar{n}_\text{fin}$ below one. Numerical values for the lowest reachable phonon numbers in all setups and systems are given in \Tref{tab:Params}.

\text{s}ubsection{Evacuated-heat flow}

Traditionally when studying refrigerators, one looks at the cooling power, namely the heat that is extracted from the object to be cooled. Here, heat is extracted from the mechanical resonator by phonons flowing into the engineered, cold environment, thereby compensating the heat flow from the hot mechanical bath. The evacuated-heat flow, constituting an equivalent to the cooling power here, hence  corresponds to the phonon flow into this cold environment times the frequency of the mechanical resonator, namely the energy that each of these phonons carries away, $J_\text{c} = \hbar\Omega_\text{mec} I_\text{c}^\text{phonon}$. This quantity can be inferred from measurements of the cavity output light field
, see Appendix~\ref{appendix:standard:flows measurement} for details.

The flux of phonons from the resonator to the cavity, $I_\text{c}^\text{phonon}$, can be identified in the evolution of the phonon number (Eq.~\text{e}qref{standard:evolution phonon number}) as
\begin{equation}
   I_\text{c}^\text{phonon} = ig\text{mec}ean{(\delta\hat{a} + \delta\hat{a}^\delta\hat{a}gger)(\delta\hat{b} - \delta\hat{b}^\delta\hat{a}gger)}.\label{cooling_power}
\text{e}nd{equation}
Note that this expression is very similar to the contribution $\text{mec}u_b$ of the entropy production rate analyzed in Ref.~\cite{brunelli_experimental_2018} with $I_\text{c}^\text{phonon} = -\text{mec}u_b(\bar{n}_\text{mec} +\frac{1}{2})$ in the steady state. The flux $-I_\text{c}^\text{phonon}$ is then equal to the flow of phonons going into the mechanical environment, which is an irreversible process that contributes to the entropy production rate and, thus, is reflected in $\text{mec}u_b$.

\begin{figure}[htb!]
    \includegraphics[width=\linewidth]{fig_cooling_power.pdf}
    \caption{\label{fig:cooling_power}
        Evacuated-heat flow as a function of cooperativity for the three setups (indicated by different lines in each plot). The insets show the evacuated-heat flow at low cooperativity with a linear scale. The number in each sub-figure corresponds to the system number in the parameter Tab.~\ref{tab:Params}. 
    }
\text{e}nd{figure}

We obtain the steady-state evacuated-heat flow in the three setups and for the four different systems and show  $J_\text{c}$ in Fig.~\ref{fig:cooling_power} as a function of the cooperativity, up to those values, where $\bar{n}_\text{fin}$ diverges. 
For small cooperativity values, the evacuated-heat flow is linear (visible in the insets) but vanishingly small, namely when the optomechanical coupling is small. The result for the evacuated-heat flow from the weak-coupling approximation (Eq.~\text{e}qref{weak-coupling:Jc}), $J_\text{c}\approx\hbar\Omega_\text{mec}\Gamma_\text{opt}\bar{n}_\text{mec}$,  clearly shows that the phonon flow for small $C$ is constrained by the small coupling to the effective cold phonon bath.
It then rapidly increases at cooperativity values at which $\Gamma_\text{opt}\approx\gamma$. This is the same parameter range for which the steady-state phonon number $\bar{n}_\text{fin}$ starts to differ from the value imposed by the hot environment, $\bar{n}_\text{mec}$. The evacuated-heat flow then saturates at a constant value, $J_\text{c}\approx\hbar\Omega_\text{mec}\gamma\bar{n}_\text{mec}$,  namely when the coupling to the effective optical bath is so large that the flow of phonons into the cold bath is only constrained by
the inflow of phonons from the hot mechanical bath. This is the maximum amount of (phononic) heat that can flow through the system, as long as the effective optical phonon bath is colder than the mechanical bath. Note that this saturation happens at a much smaller cooperativity than the one at which the minimum $\bar{n}_\text{fin}$ is reached. Indeed, the bottleneck of the heat flow through the system is reached when $\Gamma_\text{opt}$ exceeds $\gamma$; however, equilibration of the resonator at some $\bar{n}_\text{fin}$ is governed by the effective cold bath only when the exchange rate $\Gamma_\text{opt}\bar{n}_\text{opt}$ exceeds $\gamma\bar{n}_\text{mec}$.
Importantly, the evacuated-heat flow is \textit{not} affected by the backaction effect leading to an increase of $\bar{n}_\text{fin}$ at large cooperativity. This means that the opposite heat flow from the non-resonant terms $\delta\hat{a}\delta\hat{b}$ and $\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger$ exactly counterbalance each other, see  Appendix~\ref{appendix:standard:heat flows} and the plot therein. The strong hybridization between cavity photons and phonons in the limit of large cooperativity can hence be understood as leading to a direct steady-state heat flow between mechanical bath and light field, with less cooling effect on the local mechanical sub-system.  Note that this crucial difference between the behaviors of the lowest reachable phonon number and the evacuated heat flow clearly show the complementarity of these two performance quantifiers.

The described qualitative behavior is globally the same in all setups and for all studied systems. 
The threshold, which is needed to significantly cool down the resonator, is however higher for the Fano mirror, since, as discussed previously, higher cooperativities are needed to establish a $\Gamma_\text{opt}$ that is comparable to the one of the standard and the squeezed-light setup. 
The standard setup and the setup fed by squeezed light only differ by small numerical values. This means that the evacuated-heat flow is not influenced by the noise properties of the input light. Indeed, in the weak-coupling regime that governs the physics of the evacuated-heat flow, squeezing only enters via $\bar{n}_\text{opt}$, see Eq.~(\ref{fig:weak coupling cooling}), which is negligibly small both in the standard and the squeezing setup. 

\text{s}ubsection{Efficiencies}

The squeezed-light and Fano-mirror setups allow to reach lower phonon numbers than the standard setup, without altering the saturation value of the evacuated-heat flow. In this section, we aim to understand if these cooling schemes are more  or less efficient, namely whether they require extra resources to evacuate the same amount of heat and reach resonator states close to the ground state. We therefore look at the cooling efficiency $\text{e}taL$ defined as the evacuated-heat flow divided by the laser power (Eq.~\text{e}qref{etaL}) as well as a complementary cooling efficiency $\text{e}taC$, characterizing the efficiency of the phonon-photon conversion process.   

\text{s}ubsubsection{Full cooling efficiency}

In order to obtain the full cooling efficiency for all setups and systems, we insert the expression for the laser power, Eq.~\text{e}qref{Pin}, and for the evacuated-heat flow, Eq.~\text{e}qref{cooling_power}, into the full efficiency, $\text{e}taL = J_\text{c}/{P}_\text{las}$. The drive amplitude $\text{e}psilon$ of the laser, entering the expression for the laser power is determined from the steady-state amplitude of the field in the cavity,  $|\alpha| = g/g_0$, using Eq.~\text{e}qref{standard:alpha} in the standard and squeezed-light setups, and Eq.~\text{e}qref{Fano:alpha} for the Fano-mirror setup.
The resulting full cooling efficiency is plotted in Fig.~\ref{fig:etaL} as a function of the cooperativity. 

For all setups and systems, the cooling efficiency  is constant until some cooperativity threshold, namely the one at which the evacuated-heat flow displays a crossover to a plateau. The reason for this is that both the evacuated-heat flow as well as the laser power are linear in $C$ for small cooperativities. The small magnitude of the full cooling efficiency derives from the low phonon energy compared to the energy of the laser photons and the weak optomechanical single-photon coupling $g_0$. At cooperativities at which the evacuated-heat flow sets in, the full cooling efficiency decreases rapidly. We chose to cut the curves at the cooperativities, in which the phonon number $\bar{n}_\text{fin}$ diverges.

\begin{figure}[b!]
    \includegraphics[width=\linewidth]{fig_etaL.pdf}
    \caption{\label{fig:etaL}
        Full cooling efficiency $\text{e}taL$ as a function of cooperativity for the three setups (indicated by different lines in each plot). The number in each sub-figure corresponds to the system number in parameter Tab.~\ref{tab:Params}. Note that the full cooling efficiency of the squeezed-light setup differs from the plotted one by a factor $R$ as given in Eq.~\text{e}qref{eq:etaprime}.
    }
\text{e}nd{figure}

The full cooling efficiencies for the three setups have similar features, but differ by their magnitude. The largest efficiency is the one of the standard setup, since it requires the least laser-power input. For the squeezed-light setup, only the reflected fraction $R$ of the input laser power is actually used for driving the cavity, while the remaining transmitted part feeds the DPA, generating the squeezing, see Fig.~\ref{fig:Fano_and_squeezing}(a). Therefore, the actual efficiency, $\text{e}taL'$, is smaller than the plotted one, $\text{e}taL$, by a factor $R$,
\begin{equation}\label{eq:etaprime}
    \text{e}taL'=\frac{J_\text{c}}{ {P}_\text{las}/R}=R\text{e}taL
\text{e}nd{equation}
The factor $R$ by which the full cooling efficiency of the squeezed-light setup is reduced with respect to the one of the standard setup hence depends on the cost of the squeezing generation.  We refrain from giving a specific value for $R$, because quantifying this cost requires to fix the details of the experimental setup used to generate the squeezed vacuum. With Eq.~\text{e}qref{eq:etaprime}, we provide an instruction on how to include the cost of squeezing in a concrete experimental realization.

The full cooling efficiency of the Fano-mirror setup is orders of magnitude smaller than the one of the standard (and squeezed-light) setup. Indeed, as visible in Eqs.~\text{e}qref{Fano:average values}, for the Fano-mirror setup, the laser field not only creates the steady-state field in the cavity but also the one in the mirror mode. So more laser power is required to reach the same value of $\alpha$, and therefore the same cooperativity. Thus, $\text{e}taL$ illustrates how much extra resources are required to cool down the resonator using a Fano mirror. 

Overall, reaching lower phonon numbers, either with squeezed light (as one has to account for the factor $R<1$ taking into account the cost to create squeezed light) or a Fano mirror, is hence more costly and reduces the efficiency of the cooling process.

Comparing the efficiencies of the example realizations of different systems, we note that the full cooling efficiency of the non-resolved sideband system (4) is similar to systems (1) and (3), even though the  phonon number $\bar{n}_\text{fin}$ of system (4) is larger than the one for the standard and Fano-mirror setups. This motivates us to define in the following section another cooling efficiency that is more sensitive to whether the system is in the resolved-sideband regime.

\begin{figure}[tb!]
    \includegraphics[width=\linewidth]{fig_neff_etaL.pdf}
    \caption{\label{fig:loop_etaL}
        Full cooling efficiency $\text{e}taL$ as a function of the effective phonon number for the three setups (indicated by different lines in each plot). The number in each sub-figure corresponds to the system number in parameter Tab.~\ref{tab:Params}. Note that the full cooling efficiency of the squeezed-light setup differs from the plotted one by a factor $R$ as given in Eq.~\text{e}qref{eq:etaprime}.
    }
\text{e}nd{figure}

As is commonly the case in thermodynamics, there is a trade-off between optimizing the efficiency and the desired output. While in a standard refrigerator, this trade-off is typically analyzed between the efficiency and the cooling power, we here identify the lowest reachable phonon number as the relevant quantity which should be optimized, possibly together with the efficiency of the heat evacuation process. We illustrate this trade-off in Fig.~\ref{fig:loop_etaL}. Indeed, $\text{e}taL$ starts to decrease at a cooperativity orders of magnitude lower than the one at which $\bar{n}_\text{fin}$ reaches its minimum value. 

\text{s}ubsubsection{Conversion efficiency}

The full cooling efficiency, $\text{e}taL$, discussed above compares the evacuated-heat flow to the \textit{full} laser power input needed to set up the cooling mechanism. While this yields a rather complete benefit/cost ratio, it gives less insight about the cooling efficiency of the optomechanical conversion process from phononic to photonic heat. Note however, that the light leaking out of the cavity could in principle partly be reused; thus, this part would further be available as a resource. In this subsection, we therefore introduce an alternative cooling efficiency that we call \textit{conversion efficiency}.

When looking at the microscopic process, it is of interest to consider the photons involved in the optomechanical interaction as the resource. We therefore start by identifying the flux of photons interacting with the mechanics 
\begin{equation}
    I^\text{photon} = ig\left(\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger (\delta\hat{b}^\delta\hat{a}gger +  \delta\hat{b})} - \text{mec}ean{\delta\hat{a} (\delta\hat{b}^\delta\hat{a}gger+ \delta\hat{b})}\right). \label{photon flow}
\text{e}nd{equation}
We have obtained this quantity by analyzing the evolution of the photon number, Eq.~\text{e}qref{standard:evolution photon number}. Note that this expression is very similar to the contribution $\text{mec}u_a$ of the entropy production rate analyzed in Ref.~\cite{brunelli_experimental_2018} with $I^\text{photon} = \text{mec}u_a/2$ in the steady state. Again, this is not surprising as $I^\text{photon}$ describes the flux of photons that ultimately leak out of the cavity to the optical environment, which is an irreversible process that contributes to the entropy production rate. Further, also $I^\text{photon}$ can be inferred from measurements of the output cavity light field, for details see Appendix~\ref{appendix:standard:flows measurement}.

With this, the conversion efficiency can be defined as the evacuated-heat flow (phonons flowing out of the resonator into the cold bath) with respect to the flow of photons that actually interact with the mechanics, 
\begin{align}
    \text{e}taC & = \frac{I_\text{c}^\text{phonon}}{I^{\text{photon}}}\ .\label{etaC-0}
\text{e}nd{align}
Using Eqs.~\text{e}qref{cooling_power} and \text{e}qref{photon flow}, we obtain
\begin{align}
    \text{e}taC 
    &=\frac{\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}-  \delta\hat{a}\delta\hat{b}^\delta\hat{a}gger } - \text{mec}ean{\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger -  \delta\hat{a}\delta\hat{b}}}{ \text{mec}ean{\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}-  \delta\hat{a}\delta\hat{b}^\delta\hat{a}gger } + \text{mec}ean{\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger -  \delta\hat{a}\delta\hat{b}}}\ .\label{etaC}
\text{e}nd{align}
In sideband cooling schemes, the detuning of the laser drive is close to resonance with the beam-splitter terms, $\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}$ and $\delta\hat{a}\delta\hat{b}^\delta\hat{a}gger$, while the two-mode-squeezing terms, $\delta\hat{a}\delta\hat{b}$ and $\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger$, are non-resonant. The latter are detrimental to the cooling, i.e. $\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger -  \delta\hat{a}\delta\hat{b}} > 0$ (a plot of the beam-splitter and two-mode-squeezing contributions to the heat flow is included in Appendix \ref{appendix:standard:heat flows}), leading to $\text{e}taC < 1$. However when the beam-splitter terms are perfectly resonant, that is for $\Delta=\Omega_\text{mec}$ in the resolved-sideband regime the two-mode-squeezing contributions become negligible in the weak-coupling limit and the conversion efficiency goes to 1.

\begin{figure}[tb!]
    \includegraphics[width=\linewidth]{fig_etaC.pdf}
    \caption{\label{fig:etaC}
        Conversion efficiency $\text{e}taC$ as well as $\text{e}taC'$ as a function of cooperativity for the three setups (indicated by different lines in each plot). The insets show the same plots but in log-log scale. The number in each sub-figure corresponds to the system number in parameter Tab.~\ref{tab:Params}. The insets show sections of the main plots as log-log plots.
    }
\text{e}nd{figure}

Note that in the Fano-mirror setup, there is an extra term in the evolution of the photon number (Eq.~\text{e}qref{Fano:evolution photon number}) due to the added mirror mode. This extra term does neither interact with the mechanics nor is it leaking out of the cavity. If interested in this extra cost of the interaction between the cavity and the Fano-mirror mode as well, one can take this term into account as one of the resources in the efficiency, which then gives
\begin{equation}\label{etaC_Fano}
    \!\!\!\text{e}taC' = \frac{g\Im( \text{mec}ean{\delta\hat{b}^\delta\hat{a}gger (\delta\hat{a}^\delta\hat{a}gger+\delta\hat{a})})}{ g\Im( \text{mec}ean{\delta\hat{a} (\delta\hat{b}^\delta\hat{a}gger+\delta\hat{b})}) - \text{R}e(\text{mec}athcal{G}\text{mec}ean{\delta\hat{d} \delta\hat{a}^\delta\hat{a}gger}) }.\!
\text{e}nd{equation}

Figure~\ref{fig:etaC} shows the conversion efficiency, $\text{e}taC$, as a function of cooperativity for all setups and systems. 
One first of all observes that the conversion efficiency is by orders of magnitude larger than the full cooling efficiency and can even reach values up to 1. This means that the largest reduction of resources provided by the laser comes from the re-usable laser power, which is accounted for in $\text{e}taC$ but not in $\text{e}taL$. 
At the same time, it shows that it is the conversion efficiency, $\text{e}taC$, which provides the deepest insights into the fundamental physics of our cooling process.

Interestingly, while the full efficiency $\text{e}taL$ rapidly decreases at those cooperativity values at which the evacuated-heat flow starts to be considerably large, the conversion efficiency, $\text{e}taC$ stays close to maximal up to those cooperativity values, at which the phonon number, $\bar{n}_\text{fin}$ approaches its minimum.

The conversion efficiency is overall the largest for the Fano-mirror setup, while it is very similar for the squeezed-light and standard setups. This similarity between standard and squeezed-light setup reflects that the effect of squeezing is mainly a reduction of the temperature of the cold effective phonon bath. Instead, the presence of the Fano mirror improves the conversion process. 
 However, the modified conversion efficiency of the Fano-mirror setup, $\text{e}taC'$,  which takes into account the cost of the interaction between the cavity and Fano-mirror modes (light green dotted line), is reduced by up to several orders of magnitude. That can be seen best in the insets of the plot panels of Fig.~\ref{fig:etaC}.  

\begin{figure}[tb!]
 \includegraphics[width=\linewidth]{fig_neff_etaC.pdf}
 \caption{\label{fig:loop_etaC}
     Conversion efficiency $\text{e}taC$ as a function of the effective phonon number for the three setups (indicated by different lines in each plot). The number in each sub-figure corresponds to the system number in parameter Tab.~\ref{tab:Params}. 
 }
\text{e}nd{figure}

 
This is different only for system (4), which is not in the sideband-resolved regime. For system (4), as expected, the conversion efficiency is smaller than for the other systems. Surprisingly, the Fano-mirror setup is significantly more efficient than the other two setups, even when looking at $\text{e}taC'$, which means that the asymmetry in the density of states of the cavity created by the Fano mirror can provide the same evacuated-heat flow out of the non-sideband resolved system, at the cost of much less resources than what is required for the standard setup using coherent or squeezed light. 

The relation between the conversion efficiency and the smallest reachable steady-state phonon number, $\bar{n}_\text{fin}$, is plotted in Fig.~\ref{fig:loop_etaC}. While also here a trade-off can be observed, importantly, the conversion efficiency is shown to be large down to phonon numbers that are very close to the lowest reachable values.

\text{s}ection{Conclusion}\label{sec:IV}

In this work, we have provided a detailed thermodynamic analysis of optomechanical sideband cooling for different setups of cavity optomechanics---a standard setup employing a coherent laser drive or a squeezed laser drive, as well as a setup with a Fano mirror as one of the cavity mirrors. We have demonstrated how optomechanical cooling is realized by powering a heat valve that connects the mechanical resonator, heated by its coupling to a mechanical bath, to an engineered cold phonon bath provided by the optical driving. We have analyzed the lowest steady-state phonon number that can be reached for the mechanical system, together with the evacuated-heat flow and the efficiencies of the cooling process. The latter are additional, complementary performance quantifiers that we have identified in our heat-valve setting and that turn out to be particularly insightful in the analysis put forward here. 

This comprehensive theoretical analysis, which goes beyond the validity regime of previous studies of sideband cooling \cite{Asjad2016Nov}, employs realistic numbers from a broad range of possible experimental implementations \cite{thompson_strong_2008,groblacher_demonstration_2009,Eichenfield2009Nov, Chan2011,Delic2020,Delic2020_1,Rossi2018}. Our analysis has two main benefits. First, it provides an improved understanding of the underlying processes and the required mechanisms and resources for different optomechanical cooling setups. Second, the understanding of the \textit{thermodynamics} underlying sideband cooling, as discussed in this work, will enable further development of thermodynamic machines exploiting the heat flow between hot phonon bath and cold bath, engineered by the laser-light induced phonon-photon coupling. 

\text{s}ection*{Acknowledgments}

We thank Alexia Auff\`eves, Andrew Jordan and Nikolai Kiesel for helpful discussions. We acknowledge funding from the Knut and Alice Wallenberg foundation through two Wallenberg Academy fellowships (J.S. and W.W.), from the Vetenskapsr\r{a}det, Swedish VR, under project numbers 2018-05061 (J.S. and J.M.) and 2019-04946 (W.W.) and the QuantERA project C’MON-QSENS! This work was performed in the context of the Excellence Initiative Nano at Chalmers.

\appendix

\text{s}ection{Specific parameters used in the figures}

All the parameters used in the figures are gathered in Tab.~\ref{tab:Params}. The columns correspond to each of the studied systems. The uppermost section of the table contains the actual experimental parameters, the second one gives some insightful parameter ratios, the third one corresponds to the assumed parameters we used for the squeezed-light and Fano-mirror setups and the last one indicates the minimum achievable phonon number $\bar{n}_\text{fin}$ in the mechanical resonator for each of the setups.

\begin{table*}[htb!]\text{s}etlength{\tabcolsep}{3pt}\renewcommand{1.6}{1.6}
    \centering
    \begin{tabularx}{\linewidth}{>{\footnotesize}c c |X X X X }
        \hline
        \normalsize Quantity of interest&Symbol& System (1)  & System (2)  & System (3) & System (4)  \\     
        &&  &  \cite{Eichenfield2009Nov, Chan2011} &  \cite{Delic2020,Delic2020_1} &  \cite{Rossi2018} \\
        \hline\hline
          &  &\footnotesize sideband-resolved, MHz &\footnotesize sideband-resolved, GHz &\footnotesize sideband-resolved, levitated &\footnotesize non-sideband-resolved \\
        \hline\hline 
        Mechanical frequency &     $\Omega_\text{mec}/2\pi$ & $1\cdot10^6$\,Hz     &  $ 3.7\cdot 10^{9}$\, Hz  & $ 3.05\cdot 10^{5}$\, Hz & $ 1.14\cdot 10^{6}$\, Hz\\
        \hline
        Mechanical damping rate &     $\gamma/2\pi$   & 0.1\,Hz     &  $3.5\cdot 10^{4}$\,Hz   & $1.6\cdot 10^{-4}$ \,Hz & $1.1\cdot 10^{-3}\,$ Hz\\
        \hline
        Mechanical environment temperature&     $T_\text{mec}$  & 4\,K       &   $20$\,K     & $300$\,K      & $10$\,K\\
        \hline
        Average phonon number at $T_\text{mec}$ &     $\bar{n}_\text{mec}$  & $8.3\cdot 10^4$      &   $7.1\cdot10^{2}$ &   $1.3\cdot10^{8}$ &   $1.1\cdot10^{6}$\\
        \hline
        Laser wavelength&     $\lambda_\text{las}$  & 1550\,nm       &   $1537$\,nm &   $1064$\,nm &   $795$\,nm\\
        \hline
        Cavity loss rate&     $\kappa/2\pi$    & $2\cdot10^5$\,Hz     &   $5\cdot10^{8}$\,Hz  &   $1.93\cdot10^{5}$\,Hz  &   $1.59\cdot10^{7}$\,Hz\\
        \hline
        Cavity length&     $L_\text{cav}$    & $3.75\cdot10^{-3}$\,m    &   $3\cdot10^{-6}\,$m  &  $1.07\cdot10^{-2}$\,m  &   $1.6\cdot10^{-3}$\,m\\
        \hline
        Cavity free spectral range &     $\Gamma/2\pi$    &  $4.0\cdot10^{10}$\,Hz &  $1.44\cdot10^{13}$\,Hz    &  $1.40\cdot10^{10}$\,Hz  &   $9.38\cdot10^{10}$\,Hz\\
        \hline
        Single-photon coupling strength &     $g_0/2\pi$  & 10\,Hz      &   $9.1\cdot10^{5}$\,Hz &   $0.3$\,Hz &   $129$\,Hz\\
        \hline
        \hline
        Sideband resolution &     $\Omega_\text{mec}/\kappa$ & $5$      &   $7.4$ &   $1.58$ &   $0.07$\\
        \hline
        Granularity parameter&     $g_0/\kappa$ &   $5.0\cdot10^{-5}$    & $ 1.8\cdot 10^{-3}$ & $ 1.5\cdot 10^{-6}$ & $ 8.1\cdot 10^{-6}$\\
        \hline
        Mechanical quality factor &     $Q_\text{mec}=\Omega_\text{mec}/\gamma$ & $1.0 \cdot 10^7$      & $ 1\cdot 10^{5}$ & $ 1.9\cdot 10^{9}$ & $ 1.03\cdot 10^{9}$\\
        \hline
        Single-photon quantum cooperativity &     $C_0 = 2g_0^2/\kappa\gamma\bar{n}_\text{mec}$   &  $ 6.0\cdot 10^{-8}$  & $ 6.6\cdot 10^{-5}$ & $ 2.2\cdot 10^{-11}$ & $ 8.6\cdot 10^{-7}$\\
        \hline
        Number of coherent oscillations &     $Q_\text{mec}/\bar{n}_\text{mec}$   & $120.0$  & $ 140.8$ & $ 14.8$ & $ 936.7$\\
        \hline\hline
        Effective detuning  &    $\Delta=\Delta_0-g_0 \text{s}qrt{2}\bar{q}$ & $\Omega_\text{mec}$      & $\Omega_\text{mec}$  & $\Omega_\text{mec}$  & $ \kappa$\\
        \hline
        Squeezing angle  &    $\theta$ & $0.835$ rad     & $0.819$ rad  & $0.939$ rad  & $1.11$ rad\\
        \hline
        Squeezing ratio  &    $r_\text{s}$ & $0.050$     & $0.034$  & $0.154$  & $0.711$\\
        \hline
        Squeezing level\footnote{This is the amount of squeezing in dB below shot noise, with the variance of shot noise being $1/2$ and $\text{mec}ean{\Delta P^\theta_\text{in,s}}$ the variance of the squeezed quadrature.}  & $\!\!\!-10\log_{10}(2\text{mec}ean{\Delta P^\theta_\text{in,s}})$ & 0.87\,dB & 0.59\,dB & 2.7\,dB & 15.4\,dB\\
        \hline
        Fano mode loss rate&     $\gamma_d/2\pi$    & $8.0\cdot10^7$\,Hz     &   $4.38\cdot10^{11}$\,Hz  &   $7.71\cdot10^{6}$\,Hz  &   $1.31\cdot10^{6}$\,Hz\\
        \hline
        Fano setup, left-side loss rate &     $\kappa_\text{L}/2\pi$    & $8.0\cdot10^{10}$\,Hz     &   $2.88\cdot10^{13}$\,Hz  &   $2.80\cdot10^{10}$\,Hz  &   $1.88\cdot10^{11}$\,Hz\\
        \hline
        Fano setup, right-side loss rate &     $\kappa_\text{R}/2\pi$    & $5.0\cdot10^7$\,Hz     &   $8.22\cdot10^{9}$\,Hz  &   $1.75\cdot10^{8}$\,Hz  &   $5.70\cdot10^{11}$\,Hz\\
        \hline\hline
        Minimum phonon number, standard setup & $\text{mec}in(\bar{n}_\text{fin}^\text{standard})$ & 0.077& 0.021& 0.21 & 6.7 \\
        \hline
        Minimum phonon number, squeezed-light setup & $\text{mec}in(\bar{n}_\text{fin}^\text{suqeezing})$& 0.066& 0.016& 0.096& 0.039 \\
        \hline
        Minimum phonon number, Fano-mirror setup &$\text{mec}in(\bar{n}_\text{fin}^\text{Fano})$& 0.073& 0.016& 0.091& 2.4 \\
        \hline
    \text{e}nd{tabularx}
    \caption{\label{tab:Params}
        Parameters of four different optomechanical systems that we consider for our thermodynamic analysis of optomechanical cooling. The last section in the table indicates the minimum achievable steady-state phonon number $\bar{n}_\text{fin}$ for each of the setups and the section above corresponds to the assumed parameters we used in the figures. 
    }
\text{e}nd{table*}

\text{s}ection{Standard setup} \label{appendix:standard}

\text{s}ubsection{Steady state and linearization} \label{appendix:standard:lin}

The semi-classical steady state of the optomechanical system is given by Eqs.~\text{e}qref{standard:average values}, which are actually a non-linear system of equations because the effective detuning $\Delta$ depends on the steady-state mechanical position $\bar{q}$~\cite{Genes2008Mar}. As a consequence, the average photon number in the cavity, $|\alpha|^2$,  is a solution of the cubic equation
\begin{equation}
    |\alpha|^2\left(\kappa^2 + \left(\Delta_0 - \frac{2g_0^2}{\Omega_\text{mec}}|\alpha|^2\right)^2\right) = \abs{\text{e}psilon}^2.
\text{e}nd{equation}
Therefore there are two different regimes depending on the optomechanical parameters and the laser power and detuning: one where the above equation has a single real root and one where it has three. The systems (see Tab.~\ref{tab:Params})  with the respective detunings considered in this article are in the first regime, therefore there is a unique steady state $(\alpha, \bar{q})$. On the other hand, in the second regime there are multiple possible steady states, see for instance~\cite{Bowen2015Nov}, Chap.~2, for more details.

We study the evolution of the optomechanical system around this steady state by splitting the Heisenberg operators in the Langevin equations~\text{e}qref{standard:Langevin} into a steady-state value and a fluctuation operator: $\hat{a} = \alpha + \delta\hat{a}$, $\hat{q} = \bar{q} + \delta\hat{q}$ and $\hat{p} = \bar{p} + \delta\hat{p}$, with $\bar{p} = 0$. Then, we use Eqs.~\text{e}qref{standard:average values} and neglect the second order terms (in $\delta\hat{a}\delta\hat{q}$ and $\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}$) to obtain the linearized Langevin equations~\text{e}qref{standard:Langevin lin}. In these equations, the non-dissipative part of the evolution is governed by the quadratic Hamiltonian $\hat{H}_\text{lin}$, given in Eq.~\text{e}qref{H_Lin}.

\begin{figure}[htb]
    \includegraphics[width=\linewidth]{fig_Routh-Hurwitz.pdf}
    \caption{\label{fig:appendix stability}
         Stability condition $S_\text{RH}$ as a function of the cooperativities considered in this article for each of the systems from Tab.~\ref{tab:Params}. $S_\text{RH}$ has been normalized by its value at zero coupling, $S_{\text{RH}, 0} = \Omega_\text{mec}(\kappa^2 + \Delta^2)$.
    }
\text{e}nd{figure}

The stability of the optomechanical system can be determined by applying the Routh-Hurwitz criterion, from which emerges the stability condition $S_\text{RH} > 0$ \cite{Wilson-Rae2008Sep, Genes2008Mar}, with
\begin{equation}
    S_\text{RH} = \Omega_\text{mec}(\kappa^2 + \Delta^2) - 4g^2\Delta. \label{standard:stability criterion}
\text{e}nd{equation}
As can be seen in Fig.~\ref{fig:appendix stability}, the four systems studied in this article fulfill this stability condition for all the considered values of the coupling $g$.

 

\text{s}ubsection{Solving the dynamics}\label{appendix:standard:Langevin}

\text{s}ubsubsection{Solution of the Langevin equations} \label{appendix:standard:solution Langevin}

The system of equations~\text{e}qref{standard:Langevin lin} can be solved in the frequency domain~\cite{Genes2008Mar}. We have used the following convention for the Fourier transform: $\hat{f}[\omega] = \int_{-\infty}^{+\infty} \dd t \text{e}^{i\omega t}\hat{f}(t)$. In particular, the expression of the mechanical position is given by
\begin{align} \label{standard:solution Langevin}
    \chi_\text{mec}^\text{eff}[\omega]^{-1}\delta\hat{q}[\omega] = & \frac{\text{s}qrt{2\kappa} 2g\left((\kappa - i\omega)\hat{X}_\text{in}[\omega] + \Delta\hat{P}_\text{in}[\omega] \right)}{(\kappa - i\omega)^2 + \Delta^2}\nonumber \\
    &+ \text{s}qrt{\gamma}\hat{\xi}[\omega], 
\text{e}nd{align}
with the input noise quadratures
 $\hat{X}_\text{in} = (\hat{a}_\text{in} + \hat{a}_\text{in}^\delta\hat{a}gger)/\text{s}qrt{2}$ and $\hat{P}_\text{in} = i(\hat{a}_\text{in}^\delta\hat{a}gger - \hat{a}_\text{in}) /\text{s}qrt{2}$.
The mechanical susceptibility is given by
\begin{equation}\label{standard:susceptibility}
    \chi_\text{mec}^\text{eff}[\omega] = \Omega_\text{mec}\left(\Omega_\text{mec}^2 - \omega^2 - i\omega\gamma -\frac{4g^2\Delta\Omega_\text{mec}}{(\kappa -i\omega)^2 + \Delta^2}\right)^{-1}\!\!\!.
\text{e}nd{equation}
We can put it in the usual form for a harmonic oscillator~\cite{Aspelmeyer2014Dec}, that is $\chi_\text{mec}^\text{eff}[\omega] = \Omega_\text{mec}\left(\Omega_\text{mec}^\text{eff}[\omega]^2 - \omega^2 - i\omega\gamma^\text{eff}[\omega]\right)^{-1}$, identifying the effective mechanical frequency $\Omega_\text{mec}^\text{eff}[\omega] = \Omega_\text{mec} + \delta\Omega_\text{mec}[\omega]$ and damping rate $\gamma^\text{eff}[\omega] = \gamma + \Gamma_\text{opt}[\omega]$. We have denoted
\begin{align}
    \delta\Omega_\text{mec}[\omega] &= - \frac{2g^2\Delta(\kappa^2 - \omega^2 + \Delta^2)}{(\kappa^2 + (\omega - \Delta)^2)(\kappa^2 + (\omega + \Delta)^2)},
    \\
    \Gamma_\text{opt}[\omega] &= \frac{8g^2\Delta\Omega_\text{mec}\kappa}{(\kappa^2 + (\omega - \Delta)^2)(\kappa^2 + (\omega + \Delta)^2)}.\label{standard:Gopt}
\text{e}nd{align}
$\delta\Omega_\text{mec}[\omega]$ is the frequency shift due to the optical spring effect~\cite{Aspelmeyer2014Dec}.
In the weak-coupling regime, we can neglect the contributions to the evolution of the non-resonant frequencies, therefore the effective frequency from the master equation~\text{e}qref{weak-coupling:master_eq} is equal to $\Omega_\text{mec}^\text{eff}[\Omega_\text{mec}]$ and the cold bath damping rate is given by $\Gamma_\text{opt}[\Omega_\text{mec}]$.

\text{s}ubsubsection{Lyapunov equation}\label{appendix:standard:Lyapunov}

From the Langevin equations~\text{e}qref{standard:Langevin lin} and using the correlation functions of the noise, Eqs.~\text{e}qref{correlations_xi} and \text{e}qref{standard:correlations_a_in}, we obtain the evolution of the second-order moments of the quadratures:
{\allowdisplaybreaks
\begin{widetext}
    \begin{subequations}\label{standard:eqs 2nd order moments}
        \begin{align}
            \dv{\text{mec}ean{\delta\hat{X}_a^2}}{t} =\,& 2 {\Delta} \text{mec}ean{\delta\hat{X}_a\delta\hat{P}_a} - 2 \kappa \text{mec}ean{\delta\hat{X}_a^2} + \kappa, \\
            \dv{\text{mec}ean{\delta\hat{X}_a\delta\hat{P}_a}}{t} =\,& {\Delta} (\text{mec}ean{\delta\hat{P}_a^2} - \text{mec}ean{\delta\hat{X}_a^2}) + 2 g \text{mec}ean{\delta\hat{X}_a\delta\hat{q}} - 2 \kappa \text{mec}ean{\delta\hat{X}_a\delta\hat{P}_a},  \\
            \dv{\text{mec}ean{\delta\hat{X}_a\delta\hat{q}}}{t} =\,& \Omega_\text{mec} \text{mec}ean{\delta\hat{X}_a\delta\hat{p}} + {\Delta} \text{mec}ean{\delta\hat{P}_a\delta\hat{q}} - \kappa\text{mec}ean{\delta\hat{X}_a\delta\hat{q}},  \\
            \dv{\text{mec}ean{\delta\hat{X}_a\delta\hat{p}}}{t} =\,& -\Omega_\text{mec} \text{mec}ean{\delta\hat{X}_a\delta\hat{q}} + {\Delta} \text{mec}ean{\delta\hat{P}_a\delta\hat{p}} + 2 g \text{mec}ean{\delta\hat{X}_a^2}  - (\gamma + \kappa) \text{mec}ean{\delta\hat{X}_a\delta\hat{p}}, \\
            \dv{\text{mec}ean{\delta\hat{P}_a^2}}{t} =\,& -2 {\Delta} \text{mec}ean{\delta\hat{X}_a\delta\hat{P}_a} + 4 g \text{mec}ean{\delta\hat{P}_a\delta\hat{q}} - 2 \kappa \text{mec}ean{\delta\hat{P}_a^2} + \kappa, \\
            \dv{\text{mec}ean{\delta\hat{P}_a\delta\hat{q}}}{t} =\,& -{\Delta} \text{mec}ean{\delta\hat{X}_a\delta\hat{q}} + \Omega_\text{mec} \text{mec}ean{\delta\hat{P}_a\delta\hat{p}} - \kappa\text{mec}ean{\delta\hat{P}_a\delta\hat{q}}  + 2 g \text{mec}ean{\delta\hat{q}^2}, \\
            \dv{\text{mec}ean{\delta\hat{P}_a\delta\hat{p}}}{t} =\,& -{\Delta} \text{mec}ean{\delta\hat{X}_a\delta\hat{p}} - \Omega_\text{mec} \text{mec}ean{\delta\hat{P}_a\delta\hat{q}}  - (\gamma + \kappa)\text{mec}ean{\delta\hat{P}_a\delta\hat{p}} + 2 g( \text{mec}ean{\delta\hat{X}_a\delta\hat{P}_a} +  \text{mec}ean{\delta\hat{q}\delta\hat{p}}), \\
            \dv{\text{mec}ean{\delta\hat{q}^2}}{t} =\,& 2  \Omega_\text{mec} \text{mec}ean{\delta\hat{q}\delta\hat{p}}, \\
            \dv{\text{mec}ean{\delta\hat{q}\delta\hat{p}}}{t} =\,& \Omega_\text{mec} (\text{mec}ean{\delta\hat{p}^2} - \text{mec}ean{\delta\hat{q}^2}) + 2 g \text{mec}ean{\delta\hat{X}_a\delta\hat{q}} - \gamma\text{mec}ean{\delta\hat{q}\delta\hat{p}}, \\
            \dv{\text{mec}ean{\delta\hat{p}^2}}{t} =\,& -2  \Omega_\text{mec} \text{mec}ean{\delta\hat{q}\delta\hat{p}} + 4 g \text{mec}ean{\delta\hat{X}_a\delta\hat{p}}   - 2 \gamma \text{mec}ean{\delta\hat{p}^2} + {\gamma}(2 \bar{n}_\text{mec} + 1),
        \text{e}nd{align}
    \text{e}nd{subequations}
\text{e}nd{widetext}}
\noindent where we have defined the optical quadratures $\delta\hat{X}_a = (\delta\hat{a} + \delta\hat{a}^\delta\hat{a}gger)/\text{s}qrt{2}$ and $\delta\hat{P}_a = (\delta\hat{a} - \delta\hat{a}^\delta\hat{a}gger)/i\text{s}qrt{2}$.
From these equations, we get the evolution of the photon and phonon numbers, $\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}} = \frac{1}{2}(\text{mec}ean{\delta\hat{X}_a^2} + \text{mec}ean{\delta\hat{P}_a^2} - 1)$ and $\text{mec}ean{\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}} = \frac{1}{2}(\text{mec}ean{\delta\hat{q}^2} + \text{mec}ean{\delta\hat{p}^2} - 1)$:
\begin{subequations}\label{2nd_order_moments_meq}
    \begin{align}
        \dv{\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}}}{t} &= 2 g \text{mec}ean{\delta\hat{P}_a\delta\hat{q}} - 2\kappa\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{a}},\\ 
        \dv{ \text{mec}ean{\delta\hat{b}^\delta\hat{a}gger \delta\hat{b}}}{t} &= 2 g \text{mec}ean{\delta\hat{X}_a\delta\hat{p}} + \gamma \left(\bar{n}_\text{mec} + \frac{1}{2}- \text{mec}ean{\delta\hat{p}^2} \right)\nonumber\\
        &=-I^\text{phonon}_\text{c} -I^\text{phonon}_\text{h}.
    \text{e}nd{align}
\text{e}nd{subequations}
We have identified the phonon flows $I^\text{phonon}_\text{c} = - 2 g \text{mec}ean{\delta\hat{X}_a\delta\hat{p}}$ and $I^\text{phonon}_\text{h} = - \gamma \left(\bar{n}_\text{mec} + \frac{1}{2}- \text{mec}ean{\delta\hat{p}^2} \right)$ directed towards the cold and hot bath, respectively. Similarly, the flow of photons coming from the cavity and interacting with the mechanics reads $I^\text{photon} = 2 g \text{mec}ean{\delta\hat{P}_a\delta\hat{q}}$.
Rewriting Eqs.~\text{e}qref{2nd_order_moments_meq} fully in terms of the annihilation and creation operators, we obtain Eqs.~\text{e}qref{standard:evolution photon and phonon numbers}.

The linearized Hamiltonian $\hat{H}_\text{lin}$ (Eq.~\text{e}qref{H_Lin}) is Gaussian. Therefore, the system of equations~\text{e}qref{standard:eqs 2nd order moments} can be rewritten in the form of a Lyapunov equation for $V$, the covariance matrix of the quadratures,
    \begin{align}\label{standard:Lyapunov equation}
        \dv{V}{t} = A V+V A^T+ B.
    \text{e}nd{align} 
The elements of covariance matrix $V$ are defined as 
\begin{align}\label{cov matrix elements}
    V_{ij}= \dfrac{1}{2} \text{mec}ean{\{ Y_i, Y_j\}}  - \text{mec}ean{Y_i} \text{mec}ean{Y_j},
\text{e}nd{align}
where $\vec{Y}$ is the vector of operators $\left( \delta\hat{X}_a, \delta\hat{P}_a, \delta\hat{q}, \delta\hat{p}  \right)$.
We have defined the matrices\\
\begin{subequations}
    \begin{align}
    A&= \begin{bmatrix}
        - \kappa  & \Delta & 0 & 0  \\
        -\Delta&  - \kappa  & 2g & 0  \\ 
        0 & 0 & 0 & \Omega_\text{mec} \\
        2g & 0 & -\Omega_\text{mec}&  -\gamma
    \text{e}nd{bmatrix},\label{standard:Lyapunov:A}\\

    B&= \begin{bmatrix}
        \kappa & 0 & 0 & 0  \\ 
        0& \kappa  & 0 & 0  \\ 
        0& 0 & 0 & 0 \\
        0 & 0 & 0& \gamma \left(2\bar{n}+1\right) 
    \text{e}nd{bmatrix}\label{standard:Lyapunov:B}.
\text{e}nd{align}
\text{e}nd{subequations}
We then determine the steady-state $\bar{V}$ by solving $0 = A \bar{V}+\bar{V} A^T+ B$ and obtain the steady-state phonon number
\begin{equation}
    \bar{n}_\text{fin} = \frac{1}{2}(\bar{V}_{33} + \bar{V}_{44} - 1). \label{neff_Lyapunov}
\text{e}nd{equation}
Furthermore, the steady-state efficiencies (Eqs.~\text{e}qref{etaC} and \text{e}qref{etaL}) can be expressed as
\begin{align}
    \text{e}taL &= -\frac{4\kappa\Omega_\text{mec} g\bar{V}_{14}}{\omega_\text{las} \abs{\text{e}psilon}^2},\label{etaL_Lyapunov}\\
    \text{e}taC &=  -\frac{\bar{V}_{14}}{\bar{V}_{23}}\label{etaC_Lyapunov}.
\text{e}nd{align}
We therefore obtain the following analytical expressions
\begin{align}
\bar{n}_\text{fin} &= \frac{g^2 c_1 + \gamma\bar{n}_\text{mec} c_2}{c_3  S_\text{RH}},
\\
\text{e}taL &= \frac{(4\Omega_\text{mec} g_{0}^{2} {\gamma} {\kappa}^{2})(c_4 + \bar{n}_\text{mec} c_5)}{{\left({\Delta}^{2} + {\kappa}^{2}\right)} \omega_\text{las} c_3},\\
 \text{e}taC &= \frac{\gamma  S_\text{RH} (c_4 + \bar{n}_\text{mec} c_5)}{\Omega_\text{mec}(c_6 + \gamma\bar{n}_\text{mec} c_7)},
\text{e}nd{align}
with  $S_\text{RH}$ defined by Eq.~\text{e}qref{standard:stability criterion} and
{\allowdisplaybreaks
\begin{widetext}
    \begin{subequations}\label{standard:full analytical expr}
    \begin{align}
        c_1 =\,
         &\Big[ {\Omega_\text{mec}}  ({\Delta}^{2} + {\kappa}^{2})^2 + {\left({\Delta}^{2} - 4 \, {\Delta} {\Omega_\text{mec}} + 2 \, {\Omega_\text{mec}}^{2} + {\kappa}^{2}\right)} S_\text{RH} \Big] \kappa^2  \nonumber\\
        &+ \gamma \Big[{\Delta}^{2}({\Delta}^{3} + {\Omega_\text{mec}}^{3}) - 2 \, {\Delta} {\Omega_\text{mec}}^{2} g^{2} + (2 {\Delta}^{3} + 2 {\Delta}^{2} {\Omega_\text{mec}} + {\Omega_\text{mec}}^{3}) {\kappa}^{2}+ {\left({\Delta} + 2 \, {\Omega_\text{mec}}\right)} {\kappa}^{4}\nonumber\\
        &\hspace{0.8cm} + {\left({\Delta}^{2} - 3 \, {\Delta} {\Omega_\text{mec}} + {\kappa}^{2}\right)} S_\text{RH} \Big]\kappa\nonumber\\
        &+ \gamma^2\Big[ 2  {\Delta}^{2} \Omega_\text{mec} g^{2} + {\left(2  {\Delta} + \Omega_\text{mec}\right)} {\kappa}^{4} +  {\Delta}^{2}{\left(2  {\Delta} + \Omega_\text{mec}\right)} {\kappa}^{2}\Big] + \gamma^3 (\Delta^2 + \kappa^2)\kappa\Delta,\\
        c_2 =\,
        & \Big[ \Omega_\text{mec} {\kappa}^{6} + {\left(3 {\Delta}^{2} \Omega_\text{mec} + 2 \Omega_\text{mec}^{3} - 2 {\Delta} g^{2}\right)} {\kappa}^{4} + {\left((3 {\Delta}^{4}  + \Omega_\text{mec}^{4})\Omega_\text{mec} - 2(2 {\Delta}^{2}  +3  \Omega_\text{mec}^{2}) \Delta g^{2}\right)} {\kappa}^{2} \nonumber\\
        &\hspace{0.2cm}+ \Delta{\left(({\Delta}^{2} -  \Omega_\text{mec}^{2})^2\Delta\Omega_\text{mec} - 2({\Delta}^{4}- 5 {\Delta}^{2} \Omega_\text{mec}^{2}   +2  \Omega_\text{mec}^{4} ) g^{2} - 8 {\Delta} \Omega_\text{mec} g^{4}\right)}\Big] \kappa\nonumber\\
        & + \gamma\Big[ 2 \Omega_\text{mec} {\kappa}^{6} + 2 {\left(2 {\Delta}^{2} \Omega_\text{mec} + \Omega_\text{mec}^{3} - 2 {\Delta} g^{2}\right)} {\kappa}^{4}+ 2\Delta {\left(\Delta\Omega_\text{mec} ({\Delta}^{2} + \Omega_\text{mec}^{2}) -( 2 {\Delta}^{2} + 3 \Omega_\text{mec}^{2} )g^{2}\right)} {\kappa}^{2} \nonumber\\
        &\hspace{0.75cm} + 2 g^2 {\Delta}^{2} \Omega_\text{mec} \left({\Delta} \Omega_\text{mec} - 2  g^{2} \right)\Big]  +\gamma^2 \Big[\Omega_\text{mec} {\kappa}^{4} + 2 {\left({\Delta} \Omega_\text{mec} -  g^{2}\right)}\Delta {\kappa}^{2} + {\left({\Delta} \Omega_\text{mec} - 2  g^{2}\right){\Delta}^{3}}\Big]\kappa,\\
        c_3 =\,
        &8 {\Delta} \Omega_\text{mec} g^{2} {\kappa}^{2}   + \gamma \left[{\kappa}^{4} + 2 {\left({\Delta}^{2} + \Omega_\text{mec}^{2}\right)} {\kappa}^{2} + {({\Delta}^{2} - \Omega_\text{mec}^{2})^2  + 8 {\Delta} \Omega_\text{mec} g^{2}} \right]\kappa \nonumber\\
        & + 2{\gamma}^{2}  {\left[{\Delta} \Omega_\text{mec} g^{2} + ({\Delta}^{2} + \Omega_\text{mec}^{2}) {\kappa}^{2} + {\kappa}^{4}\right]} + {\gamma}^{3}  {\left({\Delta}^{2}  + {\kappa}^{2}\right)} \kappa\\
        c_4 =\,
        & - \left({\kappa}^{2} + {\left({\Delta} - \Omega_\text{mec}\right)}^2 \right){\kappa} + {\left( {\Delta} \Omega_\text{mec} + {\kappa}^{2}- {\Delta}^{2}  \right)} {\gamma} \\
        c_5 =\,
        & 2  {\Delta} \Omega_\text{mec} {\left({\gamma} + 2  {\kappa}\right)}\\
        c_6 =\,
        & 4{\left({\Delta}^{2} + {\kappa}^{2}\right)} g^{2} {\kappa}^{2} + \gamma \Big[  ({\kappa}^{2} + \Delta^2)^2 +2g^2(4\kappa^2 + \Omega_\text{mec}^2) + (\Omega_\text{mec} - 2\Delta)S_\text{RH} \Big]\kappa \nonumber \\
        & + \gamma^2 \Big[ 2{\kappa}^{4} + {\left({\Delta}^{2} + ({\Delta} - \Omega_\text{mec})^{2} + 4g^{2}\right)} {\kappa}^{2} + {\Delta}^{2} \Omega_\text{mec}(\Omega_\text{mec}-{\Delta})+ 2 {\Delta} {\left(2{\Delta} -  \Omega_\text{mec}\right)} g^{2} \Big]  + \gamma^3 {\left({\Delta}^{2} + {\kappa}^{2}\right)} {\kappa} \\
        c_7 =\,
        &  \Big[ ({\kappa}^{2} + \Delta^2)^2 + \Omega_\text{mec} S_\text{RH} \Big]2\kappa 
         + \gamma \Big[ (\Omega_\text{mec}^2 + 4\kappa^2)(\kappa^2 + \Delta^2) + \Omega_\text{mec} S_\text{RH}\Big] + \gamma^2 {\left({\Delta}^{2} + {\kappa}^{2}\right)} 2 {\kappa}\\\nonumber
    \text{e}nd{align}
    \text{e}nd{subequations}
\text{e}nd{widetext}}

When $g$ reaches the critical value 
\begin{equation}
    g_\textbf{crit} = \text{s}qrt{\frac{\Omega_\text{mec}}{4\Delta}(\kappa^2 + \Delta^2)},
\text{e}nd{equation}
the matrix $A$ is not invertible and the Lyapunov equation \text{e}qref{standard:Lyapunov equation} does not have a steady state. In particular, $\bar{n}_\text{fin}$ diverges, as visible in Fig.~\ref{fig:neff}.  This critical value of the coupling $g$ corresponds to the point at which the stability condition $S_\text{RH} > 0$ (see Appendix \ref{appendix:standard:lin}) is no longer fulfilled. 

\text{s}ubsubsection{Beyond the white-noise approximation}\label{appendix:standard:beyond white noise}

   According to~\cite{Giovannetti2001Jan}, the white-noise approximation for the mechanical noise (Eq.~\text{e}qref{correlations_xi}) gives inconsistent results, even at high temperature because the commutation relation $[\hat q, \hat p] = i$ is not preserved and it yields a spurious term in the phase fluctuation spectrum. Depending on the quantities of interest this does not necessarily impact the results, like in our case (see Appendix \ref{appendix:standard:approxs}). Nevertheless, to restore the correct commutation relation, the thermal noise operator should fulfill the following relations
   \begin{align}
     \!\! \!\text{mec}ean{\hat\xi(t)\hat\xi(t')} &= \int\!\!\dd \omega \frac{\omega\text{e}^{- i \omega(t-t')}}{2\pi} \frac{\coth(\frac{\hbar\omega}{2k_\text{B} T_\text{mec}})+1}{\Omega_\text{mec}},\!\!\\
      \!\!\! [\hat\xi(t), \hat\xi(t')] &= 2i\frac{1}{\Omega_\text{mec}}\dv{t}\delta(t-t').
   \text{e}nd{align}
   To compute the steady-state phonon number without making the white-noise approximation, we use the same method as in~\cite{Genes2008Mar}. 
   The effective phonon number is given by
       $\frac{1}{2}(\text{mec}ean{\delta\hat{q}^2} + \text{mec}ean{\delta\hat{p}^2} - 1)$,
   therefore we can compute it from the position and momentum fluctuations:
   \begin{align}
       \text{mec}ean{\delta\hat{q}^2} &= \int_{-\infty}^{+\infty}\frac{\dd \omega}{2\pi} S_q[\omega],\\
       \text{mec}ean{\delta\hat{p}^2} &= \int_{-\infty}^{+\infty}\frac{\dd \omega}{2\pi}\frac{\omega^2}{\Omega_\text{mec}^2} S_q[\omega].
   \text{e}nd{align}
   The position fluctuation spectrum,
   \begin{align}
       S_q[\omega] &=\int_{-\infty}^{+\infty} \dd t \text{e}^{i\omega t}\text{mec}ean{\delta\hat{q}(t)\delta\hat{q}(0)}\nonumber\\
       &= \int_{-\infty}^{+\infty}\frac{\dd\omega'}{2\pi}\text{mec}ean{\delta\hat{q}[\omega]\delta\hat{q}[\omega']}, \label{fluctuation spectrum}
   \text{e}nd{align}
   is obtained from the solution of the Langevin equations~\text{e}qref{standard:solution Langevin}.
   In this case, the noise correlation functions are
   \begin{equation}
       \text{mec}ean{\hat{\xi}[\omega]\hat{\xi}[\omega']} = 2\pi\frac{\omega}{\Omega_\text{mec}}\left[ 1 + \coth(\frac{\hbar\omega}{2k_\text{B} T_\text{mec}})\right] \delta(\omega + \omega'),
   \text{e}nd{equation}

   and
   \begin{subequations}\label{correlation cavity input noise}
       \begin{align}
       \hspace*{-0.3cm}\text{mec}ean{\hat{X}_\text{in}[\omega]\hat{X}_\text{in}[\omega']} &=  \text{mec}ean{\hat{P}_\text{in}[\omega]\hat{P}_\text{in}[\omega']} = \pi\delta(\omega + \omega')\!\!\\
       \hspace*{-0.3cm}\text{mec}ean{\hat{X}_\text{in}[\omega]\hat{P}_\text{in}[\omega']} &= -\text{mec}ean{\hat{P}_\text{in}[\omega]\hat{X}_\text{in}[\omega']}  =  i\pi\delta(\omega + \omega')\!\!\!
   \text{e}nd{align}
   \text{e}nd{subequations}
   So, we get 
   \begin{equation}
       S_q[\omega] = \abs{\chi_\text{mec}^\text{eff}[\omega] }^2(S_\text{th}[\omega] + S_\text{rp}[\omega])\label{Sq}
   \text{e}nd{equation}
   where
   \begin{align}
       S_\text{th}[\omega] &= \gamma\frac{\omega}{\Omega_\text{mec}}\coth(\frac{\hbar\omega}{2k_\text{B} T_\text{mec}}), \label{Sth}\\
       S_\text{rp}[\omega] &= \frac{4g^2\kappa}{\kappa^2 + (\omega - \Delta)^2},\label{Srp}
   \text{e}nd{align}
   are the thermal and radiation pressure noise spectra. To avoid lengthy computation we did not calculate the analytical expression for the phonon number, but integrated it numerically instead. We confirmed that it gives almost identical results as the white-noise approximation for the parameters considered in this article [see Fig.~\ref{fig:appendix approx comparison}].\\
   
   For the other quantities, namely $J_\text{c}$, $\text{e}taL$ and $\text{e}taC$, we would also need to compute $\text{mec}ean{\delta\hat{X}_a\delta\hat{p}}$ and $\text{mec}ean{\delta\hat{P}_a\delta\hat{q}}$. This can be done by inserting the solution of the Langevin equations~\text{e}qref{standard:Langevin lin} in the frequency domain into 
   \begin{align}
       \text{mec}ean{\delta\hat{X}_a\delta\hat{p}} &= \int_{-\infty}^{+\infty}\frac{\dd \omega}{2\pi}\int_{-\infty}^{+\infty}\frac{\dd \omega'}{2\pi}\text{mec}ean{\delta\hat{X}_a[\omega]\delta\hat{p}[\omega']},\\
       \text{mec}ean{\delta\hat{P}_a\delta\hat{q}} &= \int_{-\infty}^{+\infty}\frac{\dd \omega}{2\pi}\int_{-\infty}^{+\infty}\frac{\dd \omega'}{2\pi}\text{mec}ean{\delta\hat{P}_a[\omega]\delta\hat{q}[\omega']}.
   \text{e}nd{align}      
   Note that this approach could be used to treat a case where the input noise in the cavity is not vacuum noise by replacing Eqs.~\text{e}qref{correlation cavity input noise} by the appropriate correlation functions. 

\text{s}ubsection{Contributions to the evacuated-heat flow}\label{appendix:standard:heat flows}

The evacuated-heat flow  $J_\text{c}$, Eq.~\text{e}qref{cooling_power}, can be split into two contributions: $J_\text{c} = J_\text{c}^\text{BS} + J_\text{c}^\text{TMS}$, where $J_\text{c}^\text{BS} = i g\text{mec}ean{\delta\hat{a}^\delta\hat{a}gger \delta\hat{b} - \delta\hat{a} \delta\hat{b}^\delta\hat{a}gger}$ comes from the beam-splitter interaction term $-\hbar g (\delta\hat{a}^\delta\hat{a}gger \delta\hat{b} + \delta\hat{a} \delta\hat{b}^\delta\hat{a}gger)$ in the Hamiltonian $\hat{H}_\text{lin}$ (Eq.~\text{e}qref{H_Lin}) and $J_\text{c}^\text{TMS} = i g\text{mec}ean{\delta\hat{a} \delta\hat{b} - \delta\hat{a}^\delta\hat{a}gger \delta\hat{b}^\delta\hat{a}gger}$ from the two-mode-squeezing term $-\hbar g (\delta\hat{a} \delta\hat{b} + \delta\hat{a}^\delta\hat{a}gger \delta\hat{b}^\delta\hat{a}gger)$. In the rotating wave approximation, valid in the weak coupling regime (see next section~\ref{appendix:standard:approxs}) for $\Delta$ close to $\Omega_\text{mec}$, only the beam-splitter contribution is taken into account. Figure~\ref{fig: appendix Jc contributions} indeed shows that in the resolved-sideband regime (systems (1) to (3)), $J_\text{c}^\text{TMS}$ is negligible at low $C$, then both contributions diverge while $J_\text{c}$ remains constant. In the non-resolved-sideband regime (system (4)), $J_\text{c}$ has the same behavior but the two contributions are always of the same order of magnitude.

\begin{figure}[htb]
    \includegraphics[width=\linewidth]{fig_heat_flows_standard.pdf}
    \caption{\label{fig: appendix Jc contributions}
        Beam splitter (dashed dark blue line) and two-mode squeezing (dotted orange line) contributions to the evacuated-heat flow $J_\text{c}$ (solid cyan line) as a function of cooperativity for each of the systems from Tab.~\ref{tab:Params}. We have actually plotted $-J_\text{c}^\text{TMS}$ so as to see all the curves on the same plot in log-log scale.
    }
\text{e}nd{figure}

\text{s}ubsection{Measuring the photon and phonon flows}\label{appendix:standard:flows measurement}

In an experiment, the power spectra of the quadratures $\hat{X}_\text{out}$ and $\hat{P}_\text{out}$ of the light field leaking out of the cavity can be measured using a double homodyne detection scheme. They are related to the fluctuation spectra of the quadratures of the intra-cavity field by
\begin{align}
    S_{{X}_\text{out}}[\omega] &= S_{{X}_\text{in}}[\omega] + 4\kappa |\alpha|^2\delta(\omega)  + 2\kappa S_{{X}_a}[\omega],\\
    S_{{P}_\text{out}}[\omega] &= S_{{P}_\text{in}}[\omega] + 2\kappa S_{{P}_a}[\omega],
\text{e}nd{align}
where we have used the input-output relation $\hat{a}_\text{out} = \hat{a}_\text{in} - \text{s}qrt{2\kappa}\hat{a}$. The power spectra are defined like in Eq.~\text{e}qref{fluctuation spectrum}.

Besides, using the equations for the second-order moment  \text{e}qref{standard:eqs 2nd order moments} in the steady-state, the phonon and photon flows, $I_\text{c}^\text{phonon}$ from Eq.~\text{e}qref{cooling_power} and $I^\text{photon}$ from Eq.~\text{e}qref{photon flow}, can be expressed as functions of $\text{mec}ean{\delta\hat{X}_a^2}$ and $\text{mec}ean{\delta\hat{P}_a^2}$:
\begin{align}
   \!\! I_\text{c}^\text{phonon} 
    &= \frac{2 \kappa }{\Delta\Omega_\text{mec}ega}\left[ \kappa^2\text{mec}ean{\delta\hat{X}_a^2} - \Delta^2\text{mec}ean{\delta\hat{P}_a^2} +\frac{\Delta^2 - \kappa^2}{2}\right],\!\\
   \!\! I^\text{photon} 
    &=  \kappa \left[\text{mec}ean{\delta\hat{X}_a^2} + \text{mec}ean{\delta\hat{P}_a^2} - 1\right].
\text{e}nd{align}
For a given operator $\hat{A}$, $\text{mec}ean{\hat{A}^2} = \int_{-\infty}^{+\infty} \frac{\dd{\omega}}{2\pi} S_A[\omega]$,
so we obtain
\begin{align}
    I_\text{c}^\text{phonon} 
    =\,& \frac{ 2\kappa }{\Delta\Omega_\text{mec}ega}\int_{-\infty}^{+\infty} \frac{\dd{\omega}}{2\pi} \left(\kappa^2 S_{{X}_a}[\omega]-\Delta^2 S_{{P}_a}[\omega]\right) \nonumber\\
    & + \frac{ \kappa (\Delta^2-\kappa^2)}{\Delta\Omega_\text{mec}ega},\\
    I^\text{photon} 
    =\,&\kappa\int_{-\infty}^{+\infty} \frac{\dd{\omega}}{2\pi} \left( S_{{X}_a}[\omega] +  S_{{P}_a}[\omega]\right) - \kappa.
\text{e}nd{align}
The cavity loss rate $\kappa$ and the effective detuning $\Delta$ can be measured using optomechanically induced transparency, therefore 
$I_\text{c}^\text{phonon}$ and $I^\text{photon}$ are accessible experimentally. Finally, since the laser input power ${P}_\text{las}$ can also be determined, we  have access to the evacuated-heat flow $J_\text{c}$ and both efficiencies $\text{e}taL$ and $\text{e}taC$.

\text{s}ubsection{Comparison of the different approximations}\label{appendix:standard:approxs}

Different approaches involving various approximations are used in the literature when studying optomechanical systems. 
Here we provide a quick overview of their differences and how this impacts the steady-state phonon number in the case of the setups considered in this article.

The simplest case is to consider the weak-coupling regime (see Ref.~\cite{Marquardt2007Aug, Aspelmeyer2014Dec} and Sec.~\ref{sec:II:Thermodynamic picture}) where $g \ll \kappa, \Omega_\text{mec}$. This allows to neglect the mechanical backaction in the fluctuation spectrum of the radiation pressure force and the steady-state phonon number in the resonator is given by Eq.~\text{e}qref{weak-coupling:neff}.

When going beyond the weak coupling limit, two commonly used approaches are master equations (e.g.~\cite{Wilson-Rae2008Sep}) or Langevin equations (e.g.~\cite{Genes2008Mar, Bowen2015Nov}). 
Modeling of the Brownian motion of a mechanical resonator is a complex problem, as highlighted in~\cite{Giovannetti2001Jan, Barchielli2015Aug} for instance, and
these two approaches usually do it differently which leads to different values for $\bar{n}_\text{fin}$. Nevertheless, these differences vanish in the limit of  high temperature $T_\text{mec}$ and high quality factor $Q_\text{mec}$, which is relevant for the setups studied in this article (see Tab.~\ref{tab:Params}). A detailed derivation of both the master equation and the Langevin equations valid beyond the high temperature regime for the mechanical bath is provided  in Ref.~\cite{Barchielli2015Aug}. Typically, in the master equation approach, the Brownian motion is modeled by the dissipative term $\gamma (1+\bar{n}_\text{mec} )  D[\hat{b}]\hat{\rho}+ \gamma \bar{n}_\text{mec}  D[\hat{b}^\delta\hat{a}gger]\hat{\rho}$ in the evolution of the density operator $\hat{\rho}$ , where $D[\hat{O}]\hat{\rho} = \hat{O}\hat{\rho} \hat{O}^\delta\hat{a}gger - \frac{1}{2} \{\hat{O}^\delta\hat{a}gger \hat{O}, \hat{\rho}\}$. This term corresponds to symmetrical noise affecting identically the $\hat{q}$ and $\hat{p}$ quadratures. Conversely, in the Langevin equations, the thermal noise typically only acts on $\hat{p}$ (like in Eqs.~\text{e}qref{standard:Langevin}). However, symmetrical noise can be obtained by making a rotating wave approximation (RWA) during the derivation of the Langevin equations (see Chap.~1 in Ref.~\cite{Bowen2015Nov}).
 

\begin{figure}[htb]
    \includegraphics[width=\linewidth]{fig_neff_comparison.pdf}
    \caption{\label{fig:appendix approx comparison}
        Comparison of the resonator's steady-state phonon number $\bar{n}_\text{fin}$ obtain with different methods for each of the systems from Tab.~\ref{tab:Params}.
    }
\text{e}nd{figure}

Figure~\ref{fig:appendix approx comparison} shows the resonator's steady-state phonon number $\bar{n}_\text{fin}$ computed with different methods: the weak-coupling limit using Eq.~\text{e}qref{weak-coupling:neff} (dashed green line), the master equation approach from Ref.~\cite{Wilson-Rae2008Sep} (dotted red line), the same master equation approach but making an extra RWA (yellow dash-dotted line) and a Langevin approach, going beyond the white-noise approximation for the mechanical noise like in Appendix~\ref{appendix:standard:beyond white noise} and Ref.~\cite{Genes2008Mar} (solid blue line). 
We can see that the master equation and Langevin approaches give almost identical results for the experimental systems considered in this article. Note that the approach followed in the main text to obtain, among others, the plots of $\bar{n}_\text{fin}$ in Fig.~\ref{fig:neff} uses an intermediate level of approximation: Langevin equations with noise on the resonator's momentum only (more accurate than the symmetrical noise in the master equation) but using the white-noise approximation for the correlations of the mechanical noise. 

The green and yellow curves show us at which cooperativity the weak-coupling approximation and RWA respectively stop working. The RWA in question here consists in neglecting the non-resonant terms $\delta\hat{a}\delta\hat{b}$ and $\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}^\delta\hat{a}gger$, and therefore leads to the Hamiltonian
\begin{equation}
    \hat{H}_\text{RWA}=\hbar{\Delta} \delta\hat{a}^\delta\hat{a}gger \delta\hat{a} +\hbar \Omega_\text{mec} \delta\hat{b}^\delta\hat{a}gger \delta\hat{b}\nonumber-\hbar g (\delta\hat{a}^\delta\hat{a}gger\delta\hat{b}+\delta\hat{a}\delta\hat{b}^\delta\hat{a}gger),
\text{e}nd{equation}
which  corresponds to an intuitive scattering picture since the total number of particles is conserved.

\text{s}ection{Squeezed-light setup}\label{appendix:squeezing}

\text{s}ubsection{Squeezed noise}\label{appendix:squeezing:noise}

    
The correlation functions for the squeezed noise generated by the setup described in Sec.~\ref{sec:II:squeezed-light setup} read~\cite{Asjad2016Nov, Gardiner2004}
\begin{subequations}\label{correl_squeezing}
    \begin{align}
        \text{mec}ean{\hat{a}^{\delta\hat{a}gger}_\text{in}(t)\hat{a}_\text{in}(t')} &= \pi_\text{s}  n_\text{s}(t - t'),\\
        \text{mec}ean{\hat{a}_\text{in}(t)\hat{a}_\text{in}(t')} &= \pi_\text{s}  m_\text{s}(t - t').
    \text{e}nd{align}
\text{e}nd{subequations}
We have defined
\begin{subequations}\label{squeeezing:ns and ms}
    \begin{align}
    n_\text{s}(\tau) &= \frac{r_+^2 - r_-^2}{4}\left(\frac{\text{e}^{-r_-\abs{\tau}}}{2r_-} - \frac{\text{e}^{-r_+\abs{\tau}}}{2r_+}\right), \label{squeezing:ms}\\
    m_\text{s}(\tau) &= \frac{r_+^2 - r_-^2}{4}\left(\frac{\text{e}^{-r_-\abs{\tau}}}{2r_-} + \frac{\text{e}^{-r_+\abs{\tau}}}{2r_+}\right)\text{e}^{-i2\theta},\label{squeezing:ns}
     \text{e}nd{align} 
\text{e}nd{subequations}
where we have denoted $r_\pm = \varkappa \pm \abs{\chi}$ the decay rates of the fluctuations of the squeezed and anti-squeezed quadratures. We have assumed that we have $\varkappa > \abs{\chi}$. Then, we make the white-noise approximation (\cite{Gardiner2004}, Chap.~10). Namely, we take the limit $r_+, r_- \rightarrow \infty$ in Eqs.~\text{e}qref{squeeezing:ns and ms} to obtain the expressions of $N_\text{s}$ and $M_\text{s}$ (Eqs.~\text{e}qref{squeeezing:Ns and Ms}) used in the noise autocorrelation functions (Eqs.~\text{e}qref{squeezing:correlations_a_in}). Rewriting $r_\pm = \varkappa (1 \pm r_\text{s})$, where $r_\text{s}$ is the squeezing ratio defined in Sec.~\ref{sec:II:squeezed-light setup}, we see that the white-noise approximation can be made for any value of $r_\text{s}$ in $[0, 1[$ provided that $\varkappa$ is sufficiently large.\\    

    

\text{s}ubsection{Solving the dynamics}\label{appendix:squeezing:Langevin}

The only differences between the squeezed-light and the standard setup are the correlation functions of the optical noise, given by Eqs.~\text{e}qref{squeezing:correlations_a_in}. As a consequence, the results from Appendix~\ref{appendix:standard:solution Langevin} still hold while there are some small modifications to the evolution of the second-order moments. In particular, the effective mechanical susceptibility is unchanged (Eq.~\text{e}qref{standard:susceptibility}). This shows us that in the weak-coupling limit, the resonator's coupling strength to the effective cold bath, $\Gamma_\text{opt}$, is the same while $\bar{n}_\text{opt}$ can be made smaller than in the standard setup.

\text{s}ubsubsection{Lyapunov equation}\label{appendix:squeezing:Lyapunov}

Like in Section~\ref{appendix:standard:Lyapunov}, we derive the Lyapunov equation $\dv{V}{t} = A V+V A^T+ B^\text{squ}$ for the covariance matrix $V$ of the quadratures. $A$ is unchanged (Eq.~\text{e}qref{standard:Lyapunov:A}) while $B^\text{squ} = B + \delta B$, where $B$ is the matrix for the usual vacuum noise, given by Eq.~\text{e}qref{standard:Lyapunov:B} and $\delta B$ the additional term coming from the squeezing:
\begin{equation}
    \!\!\!\frac{\delta B}{2\kappa \pi_\text{s}} = \begin{bmatrix}
        \!N_\text{s} + \abs{M_\text{s}}\cos(2\theta)\!\! & -\abs{M_\text{s}}\text{s}in(2\theta) &\! 0 &\! 0  \\ 
        -\abs{M_\text{s}}\text{s}in(2\theta)& \!\!N_\text{s} - \abs{M_\text{s}}\cos(2\theta)  &\! 0 &\! 0  \\ 
        0& 0 &\! 0 &\! 0 \\
        0 & 0 &\! 0& \!0
    \text{e}nd{bmatrix}\!\!.\!\!\! 

\text{e}nd{equation}
\\

As previously, we determine the steady state from the Lyapunov equation and compute $\bar{n}_\text{fin}$ from Eq.~\text{e}qref{neff_Lyapunov} and the efficiencies $\text{e}taC$ and $\text{e}taL$ from Eqs.~\text{e}qref{etaC_Lyapunov} and Eq.~\text{e}qref{etaL_Lyapunov}. We obtain the following analytical expressions
\begin{align}
    \!\!\bar{n}_\text{fin} &= \frac{g^2 (c_1 + \delta c_1) + \gamma\bar{n}_\text{mec} c_2}{c_3  S_\text{RH}},\\
    \!\!\text{e}taL &= \frac{(4\Omega_\text{mec} g_{0}^{2} {\gamma} {\kappa}^{2})(c_4 + \delta c_4 + \bar{n}_\text{mec} c_5)}{{\left({\Delta}^{2} + {\kappa}^{2}\right)} \omega_\text{las} c_3},\\
    \!\! \text{e}taC &= \frac{\gamma  S_\text{RH} (c_4 + \delta c_4 + \bar{n}_\text{mec} c_5)}{\Omega_\text{mec}(c_6 + \delta c_6 + \gamma\bar{n}_\text{mec} c_7)},\!\!
\text{e}nd{align}
where $c_1,\,c_2, \,c_3, \,c_4, \,c_5, \,c_6$ and $c_7$ are given by Eqs.~\text{e}qref{standard:full analytical expr} and  $S_\text{RH}$ is defined by Eq.~\text{e}qref{standard:stability criterion}. We have defined the corrections coming from the squeezing
\begin{widetext}
    \begin{subequations}
    \begin{align}
       \delta c_1 =\,
        &  2\kappa \pi_\text{s}  N_\text{s} \Big[ \left( \Omega_\text{mec} ({\kappa}^{2} + \Delta^2)^2 + (\Delta^2 + 2\Omega_\text{mec}^2)S_\text{RH} \right)\kappa \nonumber \\
       &\hspace{1.2cm} +\gamma \left( \Omega_\text{mec}\kappa^2(2\kappa^2 + 2\Delta^2 + \Omega_\text{mec}^2)  + (\Delta^2 + \kappa^2)S_\text{RH} + \Delta\Omega_\text{mec}^2(\Delta\Omega_\text{mec} - 2g^2)\right)+\gamma^2 {\left({\Delta}^{2} + {\kappa}^{2}\right)} \Omega_\text{mec} {\kappa}\Big]\nonumber\\
        &
        -2\kappa \pi_\text{s} \abs{M_\text{s}} \left(c^*\cos(2\theta) + s^*\text{s}in(2\theta)\right),\\
        \delta c_4 =\,
        & 2 \pi_\text{s} \abs{M_\text{s}} \Big[{\left({\Delta}^{2} - \Omega_\text{mec}^{2} - {\kappa}^{2}- {\gamma} {\kappa} \right)} {\cos(2\theta)} + \Delta{\left( {\gamma} + 2  {\kappa}\right)} {\text{s}in(2\theta)}\Big]\kappa   - 2 \pi_\text{s} N_\text{s} \Big[ ({\Delta}^{2} +  \Omega_\text{mec}^{2}  +{\kappa}^{2})\kappa +  \gamma{\left({\Delta}^{2} + {\kappa}^{2}\right)} \Big] \nonumber\\
        =\,&\delta c_6,
    \text{e}nd{align}
    \text{e}nd{subequations}
    with 
    \begin{subequations}\label{squeezing:optimal angle}
        \begin{align}
            s^* = 
            & \Big( 2  {\left[\Omega_\text{mec}({\Delta}^{2} +  {\kappa}^{2})+ S_\text{RH}\right]} {\kappa} + \gamma \left[ 5\Omega_\text{mec}(4  {\kappa}^{2} + \Omega_\text{mec}^{2}) + S_\text{RH}\right]  +2 \gamma^2 \Omega_\text{mec} {\kappa} \Big){\Delta} {\kappa},\\
            c^* =\,
            &  \Big[ \Omega_\text{mec}({\Delta}^{4} - {\kappa}^{4}) + {\left({\Delta}^{2}  - 2 \, {\Omega_\text{mec}}^{2}  - {\kappa}^{2}\right)} S_\text{RH} \Big]\kappa  + \gamma\Big[ - 3 \Omega_\text{mec} {\kappa}^{4} + {\left({\Delta}^{2} \Omega_\text{mec} - \Omega_\text{mec}^{3} + 4  {\Delta} g^{2}\right)} {\kappa}^{2} + 2  {\Delta} \Omega_\text{mec}^{2} g^{2}\Big] \nonumber\\
            &+ \gamma^2{\left({\Delta} + {\kappa}\right)} {\left({\Delta} - {\kappa}\right)} \Omega_\text{mec} {\kappa}.
        \text{e}nd{align}
    \text{e}nd{subequations}
\text{e}nd{widetext}
 We can see that the stability of the optomechanical system is still given by the condition $S_\text{RH} > 0$.

\text{s}ubsubsection{Parameter optimization}\label{appendix:squeezing:optimization}

The optimal squeezing angle $\theta^*$ is obtained by finding the minimum of $\bar{n}_\text{fin}(\theta)$. We find 
$\cos(2\theta^*) = c^*/\text{s}qrt{c^{*2} + s^{*2}}$ and $\text{s}in(2\theta^*) = s^*/\text{s}qrt{c^{*2} + s^{*2}}$, where $s^*$ and $c^*$ are defined in Eqs.~\text{e}qref{squeezing:optimal angle}. Therefore the optimal squeezing angle depends only on the optomechanical parameters and not on the other squeezing parameters.

We also minimize $\bar{n}_\text{fin}(r_\text{s})$ to find the optimal squeezing ratio $r_\text{s}^*$. We find that $r_\text{s}^*$ is a root of the fourth degree equation
\begin{equation}
    a^* r_\text{s}^4 - 4b^* r_\text{s}^3+ 6a^*r_\text{s}^2 -4 b^* r_\text{s} +  a^* = 0,
\text{e}nd{equation} 
with
\begin{subequations}
    \begin{align}
       a^* =&\, c^*\cos(2\theta) + s^*\text{s}in(2\theta)\\
       b^*=
       &\Big[\Omega_\text{mec} ({\kappa}^{2} + \Delta^2)^2+ S_\text{RH}(\Delta^2 + \kappa^2 + 2\Omega_\text{mec}^2) \Big]\kappa \nonumber\\
       &+\gamma \Big[ \Omega_\text{mec}(\kappa^2 + \Delta^2) (2\kappa^2 + \Delta^2) +(\Omega_\text{mec}^2 + \kappa^2)S_\text{RH}\nonumber\\
       &\hspace{0.75cm}   - 2\Delta^3g^2\Big]+\gamma^2{\left({\Delta}^{2} + {\kappa}^{2}\right)} \Omega_\text{mec} {\kappa}
    \text{e}nd{align}
\text{e}nd{subequations}
The only root that fulfills $0 \le r_\text{s} < 1$ is 
\begin{align}
    r_\text{s}^* =\,&  \frac{b^*}{a^*} + \text{s}qrt{\left(\frac{b^*}{a^*}\right)^{2} - 1}\nonumber\\
    &-\text{s}qrt{2} \text{s}qrt{\left(\frac{b^*}{a^*}\right)^{2} +\frac{b^*}{a^*} \text{s}qrt{\left(\frac{b^*}{a^*}\right)^{2} - 1}  - 1},
\text{e}nd{align}
which depends on both the optomechanical parameters and the squeezing angle.\\

\text{s}ubsubsection{Beyond the white-noise approximation}\label{appendix:squeezing:beyond}

Like for the standard setup in Appendix~\ref{appendix:standard:beyond white noise}, we can go beyond the white-noise approximation, but this time for both noises. The correlation functions for the optical noise, in the frequency domain, read
\begin{subequations}
    \begin{align}
        \text{mec}ean{(\hat{a}^{\delta\hat{a}gger}_\text{in})[\omega]\hat{a}_\text{in}[\omega']} &= \pi_\text{s} n_\text{s}[\omega] 2\pi\delta(\omega + \omega'),\\
        \text{mec}ean{\hat{a}_\text{in}[\omega]\hat{a}_\text{in}[\omega']} &= \pi_\text{s} m_\text{s}[\omega] 2\pi\delta(\omega + \omega'),
    \text{e}nd{align}
\text{e}nd{subequations}
with
\begin{subequations}
    \begin{align}
        n_\text{s}[\omega] &= \frac{r_+^2 - r_-^2}{4} \left(\frac{1}{r_-^2 + \omega^2} - \frac{1}{r_+^2 + \omega^2} \right),\\
        m_\text{s}[\omega] &= \frac{r_+^2 - r_-^2}{4} \left(\frac{1}{r_-^2 + \omega^2} + \frac{1}{r_+^2 + \omega^2} \right)\text{e}^{-i2\theta}.
    \text{e}nd{align}
\text{e}nd{subequations}

As previously, the steady-state phonon number $\bar{n}_\text{fin}$ can be computed using the position fluctuation spectrum
\begin{equation}
    S_q[\omega] = \abs{\chi_\text{mec}^\text{eff}[\omega] }^2(S_\text{th}[\omega] + S^\text{squ}_\text{rp}[\omega]).\label{squeezing:Sq}
\text{e}nd{equation}
$S_\text{th}[\omega]$ is given by Eq.~\text{e}qref{Sth} and the radiation pressure noise spectrum reads
\begin{widetext}
    \begin{align}
        \hspace*{-0.25cm}S^\text{squ}_\text{rp}[\omega] =\,& S_\text{rp}[\omega]\left(1 +2\pi_\text{s}\frac{ n_\text{s}[\omega] (\kappa^2 + \Delta^2 + \omega^2) - \abs{m_\text{s}[\omega]}\left[( \Delta^2 - \kappa^2 - \omega^2)\cos(2\theta) + 2\Delta\kappa\text{s}in(2\theta)\right]}{\kappa^2 + (\omega + \Delta)^2}\right)\!\!\!\!
    \text{e}nd{align}
\text{e}nd{widetext}
where $S_\text{rp}[\omega]$ is the radiation pressure noise spectrum for unsqueezed vacuum noise (Eq.~\text{e}qref{Srp}).

\text{s}ection{Fano-mirror setup}\label{appendix:Fano}

\text{s}ubsection{Steady state and linearization}\label{appendix:Fano:lin}

We determine the semi-classical steady state of Eqs.~\text{e}qref{Fano:Langevin} and derive the Langevin equations~\text{e}qref{Fano:Langevin lin} like for the standard setup in Appendix~\ref{appendix:standard:lin}. The difference is the presence of the Fano-mirror mode, whose annihilation operator is split into $\hat{d} =\delta + \delta\hat{d}$. The steady state $(\alpha,  \bar{q}, \delta)$ is given by the system of non-linear equations~\text{e}qref{Fano:average values}.\\

In this case, the stability condition is $S_\text{RH}^\text{Fano} > 0$, with 

\begin{align}
    S_\text{RH}^\text{Fano} =\,& {\Delta}^{2} \Omega_\text{mec} {\gamma_d}^{2} + 4  \Omega_\text{mec} {\gamma_d}^{2} ({\kappa_0}^{2} -  {\Delta}   \text{s}qrt{{\kappa_0 \kappa_\text{L}}} + {\kappa_0} {\kappa_\text{L}}) \nonumber\\
    &  + 2   \Omega_\text{mec} {\gamma_d}(\kappa_\text{L} - {\kappa_0}) {\Delta_d}({\Delta} + 2  \text{s}qrt{\kappa_0{\kappa_\text{L}}} ) \nonumber\\
    & +  \Omega_\text{mec} {\Delta_d}^{2}\left({\Delta}^{2} +  ( {\kappa_0} + {\kappa_\text{L}})^2\right) 
     \nonumber\\
    & - 4  \big( {\gamma_d}^{2}(\Delta - 2 \text{s}qrt{{\kappa_0}{\kappa_\text{L}}})+ {\gamma_d} ({\kappa_\text{L}} - {\kappa_0}) {\Delta_d}  \nonumber\\
        & \hspace{0.75cm}+ {\Delta} {\Delta_d}^{2}\big) g^{2}.
\text{e}nd{align}
$S_\text{RH}^\text{Fano} > 0$ corresponds, up to a factor $\Omega_\text{mec}$, to the determinant of the $A$ matrix (Eq.~\text{e}qref{A_Fano}) of the Lyapunov equation (see Section \ref{appendix:Fano:Lyapunov}). 
As can be seen in Fig.~\ref{fig:appendix stability Fano}, the Fano-mirror setup for the four systems studied in this article fulfill this stability condition for all the considered values of the coupling $g$.

\begin{figure}[htb]
    \includegraphics[width=\linewidth]{fig_Routh-Hurwitz_Fano.pdf}
    \caption{\label{fig:appendix stability Fano}
         Stability condition $S_\text{RH}^\text{Fano}$ of the Fano-mirror setup as a function of the cooperativities considered in this article for each of the systems from Tab.~\ref{tab:Params}. $S_\text{RH}^\text{Fano}$ has been normalized by its value at zero coupling $S_{\text{RH}, 0}^\text{Fano}$.
    }
\text{e}nd{figure}

\text{s}ubsection{Solving the dynamics}\label{appendix:Fano:Langevin}

\text{s}ubsubsection{Solution of the Langevin equations}

Like in Appendix~\ref{appendix:standard:solution Langevin}, the Langevin equations~\text{e}qref{Fano:Langevin lin} can be solved in the frequency domain and the expression of the mechanical position is given by
\begin{align} \label{Fano:solution Langevin}
    \chi_\text{mec}^\text{eff}[\omega]^{-1}\delta\hat{q}[\omega] = \;
    & g \text{s}um_{\text{mec}u = \text{L},\text{R}} \Big\{(t_\text{mec}u[\omega] + t^*_\text{mec}u[-\omega])\hat{X}_{\text{in},\text{mec}u}[\omega] \nonumber \\
    & \hspace{1.25cm}+i(t_\text{mec}u[\omega] - t^*_\text{mec}u[-\omega])\hat{P}_{\text{in},\text{mec}u}[\omega] \Big\}\nonumber \\
    &+ \text{s}qrt{\gamma}\hat{\xi}[\omega], 
\text{e}nd{align}
with the input noise quadratures
$\hat{X}_{\text{in},\text{mec}u} = (\hat{a}_{\text{in},\text{mec}u} + \hat{a}^\delta\hat{a}gger_{\text{in},\text{mec}u})/\text{s}qrt{2}$ and $\hat{P}_{\text{in},\text{mec}u} = i(\hat{a}^\delta\hat{a}gger_{\text{in},\text{mec}u} - \hat{a}_{\text{in},\text{mec}u}) /\text{s}qrt{2}$.
We have defined $t_\text{mec}u[\omega]$ as
\begin{align}
    t_\text{L}[\omega] &= \frac{\text{s}qrt{2\kappa_\text{L}} - \text{s}qrt{2\gamma_d}\text{mec}athcal{G}\varepsilon_d[\omega]^{-1}}{\varepsilon_a[\omega] - \text{mec}athcal{G}^2\varepsilon_d[\omega]^{-1}},\\
    t_\text{R}[\omega] &= \frac{\text{s}qrt{2\kappa_0} }{\varepsilon_a[\omega] - \text{mec}athcal{G}^2\varepsilon_d[\omega]^{-1}},
\text{e}nd{align}
with $\varepsilon_a[\omega] = \kappa + i(\Delta - \omega)$ and $\varepsilon_d[\omega] = \gamma_d + i(\Delta_d - \omega)$. Note that in this section the $*$ exponent denotes the complex conjugate and not an optimal value like in the part on the squeezing.
The mechanical susceptibility is given by 
\begin{equation}
    \chi_\text{mec}^\text{eff}[\omega]^{-1} = \chi^0_\text{mec}[\omega]^{-1} + \chi^\text{eff}_\text{opt}[\omega]^{-1},
\text{e}nd{equation}
where $\chi^0_\text{mec}[\omega] = \Omega_\text{mec}\left(\Omega_\text{mec}^2 - \omega^2 - i\omega\gamma\right)^{-1}$ is the mechanical susceptibility of the bare resonator and $\chi^\text{eff}_\text{opt}[\omega]$ the optical contribution to the susceptibility, which reads
\begin{equation}
   \chi^\text{eff}_\text{opt}[\omega]^{-1} =  \frac{2ig^2 }{\varepsilon_a^*[-\omega] - \text{mec}athcal{G}^{*2}\varepsilon_d^*[-\omega]^{-1}} -\frac{ 2ig^2 }{\varepsilon_a[\omega] - \text{mec}athcal{G}^2\varepsilon_d[\omega]^{-1}}.
\text{e}nd{equation}
We then identify the optical contribution to the mechanical damping rate,
\begin{equation}
    \Gamma_\text{opt}[\omega] = -\frac{\Omega_\text{mec}}{\omega} \Im{\chi^\text{eff}_\text{opt}[\omega]^{-1}}.
\text{e}nd{equation}
Therefore, in the weak-coupling regime, the resonator's coupling to the effective cold bath, given by $\Gamma_\text{opt}[\Omega_\text{mec}]$, is different from the one of the standard setup.

\text{s}ubsubsection{Lyapunov equation}\label{appendix:Fano:Lyapunov}

Like in Section~\ref{appendix:standard:Lyapunov}, we write the Langevin equations for the quadratures and obtain the Lyapunov equation $\dv{V}{t} = A V+V A^T+ B$ for $\vec{Y} = \left( \delta\hat{X}_a, \delta\hat{P}_a, \delta\hat{q}, \delta\hat{p}, \delta\hat{X}_d, \delta\hat{P}_d  \right)$, with
{ \footnotesize   \noindent
   \begin{subequations}
        \begin{align}
       \!\! \!\!\!  A &=\! \begin{bmatrix}
            -{\kappa} &\! {\Delta} &\! 0 &\! 0 &\!\!-\text{s}qrt{{\gamma_d}{\kappa_\text{L}}} &\! \text{s}qrt{{\gamma_d}{\kappa_0}} \\
            -{\Delta} &\! -{\kappa} &\! 2  g &\! 0 &\! \!-\text{s}qrt{{\gamma_d}{\kappa_0}} &\! -\text{s}qrt{{\gamma_d}{\kappa_\text{L}}} \\
            0 &\! 0 &\! 0 &\! \Omega_\text{mec} &\! 0 &\! 0 \\
            2  g &\! 0 &\!\! -\Omega_\text{mec} &\! -{\gamma} &\! 0 &\! 0 \\
            -\text{s}qrt{{\gamma_d}{\kappa_\text{L}}} &\! \text{s}qrt{{\gamma_d}{\kappa_0}} &\! 0 &\! 0 &\!-{\gamma_d} &\! {\Delta_d} \\
            -\text{s}qrt{{\gamma_d}{\kappa_0}} &\! -\text{s}qrt{{\gamma_d}{\kappa_\text{L}}} &\! 0 &\! 0 &\! -{\Delta_d} &\! -{\gamma_d}
        \text{e}nd{bmatrix}\label{A_Fano}\!,\!\!\\
        \!\!\! \! \! B &=\!\begin{bmatrix}
            {\kappa} & 0 & 0 & 0 &\!\text{s}qrt{{\gamma_d}{\kappa_\text{L}}} \!& 0 \\
            0 & {\kappa} & 0 & 0 & 0 & \!\text{s}qrt{{\gamma_d}{\kappa_\text{L}}}\! \\
            0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & {\left(2 {\bar{n}_\text{mec}} + 1\right)} {\gamma} \!& 0 & 0 \\
            \text{s}qrt{{\gamma_d}{\kappa_\text{L}}} & 0 & 0 & 0 & {\gamma_d} & 0 \\
            0 & \text{s}qrt{{\gamma_d}{\kappa_\text{L}}} & 0 & 0 & 0 & {\gamma_d}
        \text{e}nd{bmatrix}\!.\!\!
    \text{e}nd{align}
   \text{e}nd{subequations}
}
These results are in agreement with the ones from the supplementary material of~\cite{Cernotik2019Jun}.
The matrices are now too big to get a tractable analytical solution so we solve the Lyapunov equation numerically instead
(using python and in particular \verb|scipy.linalg.solve_continuous_lyapunov(A, -B)|). As previously, $\bar{n}_\text{fin}$ is given by Eq.~\text{e}qref{neff_Lyapunov} and the efficiencies $\text{e}taC$ and $\text{e}taL$ by Eqs.~\text{e}qref{etaC_Lyapunov} and~\text{e}qref{etaL_Lyapunov}. $\text{e}taC'$ (Eq.~\text{e}qref{etaC_Fano}) is obtained from
\begin{equation}
    \text{e}taC' = \frac{-2 g\bar{V}_{14}}{2g\bar{V}_{23}- \text{s}qrt{\gamma_d\kappa_\text{L}}(\bar{V}_{15} + \bar{V}_{26})+\text{s}qrt{\gamma_d\kappa_0}(\bar{V}_{16} - \bar{V}_{25})}.
\text{e}nd{equation}

\text{s}ubsubsection{Beyond the white-noise approximation}\label{appendix:Fano:beyond}

Like for the standard setup in Appendix~\ref{appendix:standard:beyond white noise}, we can go beyond the white-noise approximation for the mechanical noise and compute the effective phonon number from position fluctuation spectrum. Using Eq.~\text{e}qref{Fano:solution Langevin}, we obtain 
\begin{equation}
    S_q[\omega] = \abs{\chi_\text{mec}^\text{eff}[\omega] }^2(S_\text{th}[\omega] + S^\text{Fano}_\text{rp}[\omega])\label{Fano:Sq}
\text{e}nd{equation}
where $S_\text{th}[\omega]$ is given by Eq.~\text{e}qref{Sth} and the radiation pressure noise spectrum reads
\begin{equation}
    S^\text{Fano}_\text{rp}[\omega] = 2g^2(\abs{t_\text{L}[\omega]}^2+\abs{t_\text{R}[\omega]}^2).\label{Fano:Srp}
\text{e}nd{equation}

\text{s}ubsubsection{Weak-coupling limit}\label{appendix:Fano:weak-coupling}

We study here the weak-coupling limit $g \ll \kappa, \Omega_\text{mec}$. In this limit, illustrated in Fig.~\ref{fig:Fano_and_squeezing}(c), the Stokes and anti-Stokes rates are given by 
\begin{equation}
    A^\text{Fano}_\pm = g^2(\abs{t_\text{L}[\text{mec}p\Omega_\text{mec}]}^2+\abs{t_\text{R}[\text{mec}p\Omega_\text{mec}]}^2). \label{Fano:Stokes rates}
\text{e}nd{equation}
Therefore, at $\Delta_d = \Omega_\text{mec}$ and the optimal $\gamma_d$, which is relevant for the sideband-resolved systems (1)-(3) considered in this article (see \Tref{tab:Params}), we have
\begin{align}
   \!\! \frac{A^\text{Fano}_-}{A_-} &= \frac{(4\Omega_\text{mec})^2}{\Gamma\kappa_\text{eff}},\\
    \!\!\frac{A^\text{Fano}_+}{A_+} &= \frac{2\Gamma\kappa_\text{eff}(4\Omega_\text{mec})^2(4\Omega_\text{mec}^2 + \kappa_\text{eff}^2)}{(4\Omega_\text{mec})^2\left((4\Omega_\text{mec})^2 + 2\Gamma\kappa_\text{eff}\right)^2 + \Gamma^2\kappa_\text{eff}^4},
\text{e}nd{align}
where $A_\pm$ are the rates (Eqs.~\text{e}qref{Stokes_rates_perturbative}) for the cavity in the standard setup of equivalent linewidth and associated loss rate $\kappa_\text{eff}$ (see Sec.~\ref{sec:II:Fano-mirror setup}). Namely, in systems (1) to (3), the use of the Fano-mirror setup reduces the rate of the detrimental Stokes process by 10 to 12 orders of magnitude while the rate of the anti-Stokes process that cools down the resonator is only decreased by 2 to 4 orders of magnitude. As a consequence, the temperature of the effective cold phonon bath is lower ($\bar{n}_\text{opt}$ is 8 to 10 orders of magnitude smaller) but $\Gamma_\text{opt}$ is also a couple of orders of magnitude smaller. This explains the shift of curves for the Fano-mirror setup in terms of cooperativity in Fig.~\ref{fig:neff}. Indeed, expressing $\Gamma_\text{opt}$ as a function of the cooperativity, $\Gamma_\text{opt} = \alpha_\text{opt} C$, and neglecting $\Gamma_\text{opt}\bar{n}_\text{opt}$ in Eq.~\text{e}qref{weak-coupling:neff}, we obtain $\bar{n}_\text{fin} \text{s}imeq \gamma \bar{n}_\text{mec} /(\gamma + \alpha_\text{opt} C)$. The reduction of the phonon number starts to be non-negligible when $C$ becomes larger than $\gamma / \alpha_\text{opt}$, but $\alpha_\text{opt}$ is a few orders of magnitude smaller for the Fano-mirror setup, hence the shift in cooperativity. For a system in the resolved-sideband regime, with $\Delta = \Omega_\text{mec}$, the cooperativity threshold to start seeing a decrease in $\bar{n}_\text{fin}$ is $C \text{s}im 1/\bar{n}_\text{mec}$ for the standard setup and $C \text{s}im \Gamma\kappa_\text{eff} /((4\Omega_\text{mec})^2\bar{n}_\text{mec})$ for the Fano-mirror setup.

\text{s}ubsubsection{Optimal detuning in the non-resolved-sideband regime}\label{appendix:Fano:non-sideband-resolved}

In the non-sideband-resolved regime, for the standard and squeezed-light setups, the best cooling is not obtained for a detuning $\Delta$ close to $\kappa$. However, for the Fano-mirror setup, the optimal choice is $\Delta_d \text{s}imeq 2.5\Omega_\text{mec}$, not $\Delta_d = \kappa_\text{eff}$. We can understand it by looking at the phonon number in the cold bath in the weak-coupling regime: 
$\bar{n}_\text{opt}^\text{Fano} = A^\text{Fano}_+/(A^\text{Fano}_- - A^\text{Fano}_+)$ (see Eq.~\text{e}qref{Fano:Stokes rates}). Indeed, the parameters of the Fano-mirror setup are such that $\gamma_d = 4\Omega_\text{mec}\zeta_0$ and $\zeta_0 = 4\Omega_\text{mec} / \kappa_\text{eff}$ (see Sec.~\ref{sec:II:experimental setups:Fano}), so in the limit $\Omega_\text{mec} \ll \kappa_\text{eff}$ and $\gamma_d \ll \Gamma$, we find that the minimum of $\bar{n}_\text{opt}^\text{Fano}$ is reached for 
\begin{align}
    \frac{\Delta_d}{\Omega_\text{mec}} &= \frac{{\left(3 \, \text{s}qrt{37} + 26\right)}^{\frac{2}{3}} + 2 \, {\left(3 \, \text{s}qrt{37} + 26\right)}^{\frac{1}{3}} + 7}{3 \, {\left(3 \, \text{s}qrt{37} + 26\right)}^{\frac{1}{3}}} \nonumber\\
    &\text{s}imeq 2.5.\\\nonumber
\text{e}nd{align} 

\text{s}ection{Comparison with absorption refrigerators} \label{appendix:comp absorption refrigerator}

Optomechanical cooling has also been studied in the context of absorption refrigerators \cite{mitchison_realising_2016, Naseem2020Jun}, which are, like the setups we study, autonomous thermal machines.
However, this framework is different from the situation we consider in this article as it typically assumes a three-body interaction, where three system parts interact with three heat baths \cite{Mitchison2018}. One of the baths has to be hotter than the other and can, for example, be realized by thermal light. 

Conversely, in our case, we study a cavity optomechanical system, which is a two-body system (cavity and mechanical resonator), where each subsystem is coupled to a bath. Those baths are typically at the same temperature. Furthermore, we analyze coherent or squeezed light incident to the cavity - which are both minimum uncertainty states - and therefore interpreted as work. The valve picture depicted in Fig.~\ref{fig:principle}(a) is quite intuitive to describe the standard optomechanical setup and can be straightforwardly adapted to the squeezed-light and Fano-mirror setups. In contrast, the result of a quantum absorption refrigerator is that one can cool with thermal light, which we do not consider in our work.

The first setup studied in Ref.~\cite{mitchison_realising_2016} couples the cavity mode to the motion of a trapped atom (mechanical degree of freedom) and to its internal state (two-level system) while the second couples the modes of two cavities to the motion of a trapped atom. Both cases are quite different from the standard optomechanical sideband cooling setup we study and therefore hard to compare. We can nevertheless note that their predictions using sunlight (Fig.~7 in Ref.~\cite{mitchison_realising_2016}) brings the mechanics close to the ground state but $\bar{n}_\text{fin}$ does not go well below 1 as in the levitated optomechanical system which is the most similar system we study (see system (3) in Fig.~\ref{fig:neff} and Ref.~\cite{Delic2020} for the experimental results).

The setup studied in Ref.~\cite{Naseem2020Jun} would correspond to our standard setup where the coherent laser drive has been replaced by thermal light at temperature $T_w$. But their results are hard to compare to ours as they do not study the final phonon number and renormalize their efficiency by the Carnot absorption refrigerator efficiency, which is ill-defined in our two-terminal system with possibly equal temperatures.

\text{e}nd{document}