# Photon cooling: linear vs nonlinear interactions

A. Hovhannisyan1, V. Stepanyan2 and A.E. Allahverdyan3 1Institute of Applied
Problems of Physics, Yerevan, Armenia
2Physics Department, Yerevan State University, Yerevan, Armenia
3Alikhanian National Laboratory (Yerevan Physics Institute), Yerevan Armenia

###### Abstract

Linear optics imposes a relation that is more general than the second law of
thermodynamics: For modes undergoing a linear evolution, the full mean
occupation number (i.e. photon number for optical modes) does not decrease,
provided that the evolution starts from a (generalized) diagonal state. This
relation connects to noise-increasing (or heating), and is akin to the second
law and holds for a wide set of initial states. Also, the Bose-entropy of
modes increases, though this relation imposes additional limitations on the
initial states and on linear evolution. We show that heating can be reversed
via nonlinear interactions between the modes. They can cool—i.e. decrease the
full mean occupation number and the related noise—an equilibrium system of
modes provided that their frequencies are different. Such an effect cannot
exist in energy cooling, where only a part of an equilibrium system is cooled.
We describe the cooling set-up via both efficiency and coefficient of
performance, and relate the cooling effect to the Manley-Rowe theorem in
nonlinear optics.

##  I Introduction

Cooling is needed for noise-reduction and for capturing quantum degrees of
freedom. It has been studied during the past 100 years in various set-ups
Sheik-Bahae and Epstein (2007); Abragam and Goldman (1978); Walls and Milburn
(2007). Cooling processes are also fundamental for thermodynamics: they
sharpen the understanding of the second law, and are instrumental for the
third law E. B. Stuart (1970). An interesting example of this is the laser
cooling of solids via the anti-Stokes effect, which does have both quantum and
thermodynamic nature Sheik-Bahae and Epstein (2007). Much attention is
currently devoted to cooling processes in quantum thermodynamics Mahler
(2014); Silva et al. (2016); Wilming and Gallego (2017); Clivaz et al. (2019);
Freitas et al. (2018); Raeisi (2021); Gelbwaser-Klimovsky et al. (2015); Uzdin
et al. (2015); Liuzzo-Scorpo et al. (2016); Long and Liu (2016); Gonzalez-
Ayala et al. (2018); Singh et al. (2020); Raeisi and Mosca (2015); Taranto et
al. (2021); Allahverdyan et al. (2011, 2010). It is known that only a part of
a thermally isolated (initially equilibrium) system can be cooled in terms of
energy (or temperature), that cooling such systems costs high-graded energy
(work), hence the definition of the coefficient of performance (COP), and that
cooling is limited by energy spectra and complexity costs.

Here we consider bosonic (for clarity photonic) degrees of freedom (modes),
and show that linear transformations (e.g. linear optics) always increase the
full photon number of the system. This statement holds for a wide class of
initial states. For such states increasing the mean photon number relates to
increasing of noise (heating). The heating is more general than the second
law. To confirm this point, we studied the full Bose-entropy of modes. This
coarse-grained entropy is conditionally maximal at equilibrium, and can change
under a unitary evolution, in contrast to the fine-grained von Neumann
entropy. We show that the Bose-entropy can increase, but this relation (a
formulation of the second law) demands additional limitations both on the
initial states and linear evolution.

Heating can be reversed by nonlinear interactions. One can cool in this sense
an initially equilibrium system, which consists of two or higher number of
modes. This is not possible for energy cooling, where as demanded by the
second law, only subsystem’s energy can be decreased (cooled). Our cooling
set-up is characterized by two efficiency-like parameters: the coefficient of
performance (COP) and the efficiency. The former refers to energy costs of
cooling, while the latter normalizes the cooling result over the total changes
introduced in the system. Nonlinear interactions achieve cooling in near-
resonance regimes, where there is an effective conservation law in the number
of photons (Manley-Rowe theorem) Weiss (1957); Landau and Lifshitz (2013).
Thus, this cooling scenario uncovers a thermodynamic role of nonlinear optical
processes. We work in terms of photons, but our results hold for other bosons
(e.g. phonons).

This paper is organized as follows. The second section shows that the mean
boson (photon) number increases in linear evolution if the evolution starts
from a certain class of generalized diagonal initial states. This class is
sufficiently large and includes usual diagonal states (in the Fock basis),
independent states (over the modes) etc. Section also relates the increase of
the mean number to noise and formulates this as a heating (no-cooling)
principle for linear evolution. In Section III, we study the Bose-entropy for
modes and explain under what additional restrictions (compared with the mean
photon number increase) this entropy grows. This section also addresses the
physical meaning of the Bose-entropy. Section IV describes the optimal cooling
set-up for two modes and introduces the basic characteristics of cooling,
viz., efficiency and the coefficient of performance (or COP). This section
emphasizes the key feature of this cooling setup, namely: the global cooling
of an equilibrium system in terms of the mean photon number (and noise) is
possible. Section V demonstrates that cooling is possible also via a feasible
nonlinear two-mode interaction, works out a simple example of such
interactions, and establishes relations with the Manley-Rowe theorem, a known
result in nonlinear physics. The last section summarizes our results.

##  II No cooling for linear interactions

###  II.1 A single bosonic mode

Linear processes, which are described with Hamiltonians quadratic in
creation/annihilation operators, describe the lion’s share of boson dynamics
Caves (1982); Garrison and Chiao (2008). Consider the simplest example of such
processes: a single mode that underwent a linear evolution governed by a
quadratic Hamiltonian. In the Heisenberg picture, the general form of this
evolution connects initial $a=a(0)$ and final $b=a(t)$ annihilation operator
of the mode:

$b=Sa+Ra^{\dagger}+f,$ (1)

where $S$, $R$ and $f$ are complex c-numbers that characterize the evolution.
The initial state (density matrix) $\rho$ of the mode satisfy:

$\displaystyle\langle a\rangle\equiv{\rm tr}(a\rho)=\langle
a^{\dagger}\rangle=0.$ (2)

The commutation relation $[b,b^{\dagger}]=[a,a^{\dagger}]=1$ impose
$|S|^{2}-|R|^{2}=1$ in (1). Then we get from (1, 2)

$\langle b^{\dagger}b\rangle-\langle a^{\dagger}a\rangle=2|R|^{2}\langle
a^{\dagger}a\rangle+|R|^{2}+|f(t)|^{2}\geq 0,$ (3)

i.e. the mean photon number difference defined in the LHS of (3) can only
increase. In particular, this conclusion holds for linear amplifiers Caves
(1982). According to (1) also the dispersion of the photon number increases:

$\displaystyle\langle(b^{\dagger}b)^{2}\rangle-\langle
b^{\dagger}b\rangle^{2}\geq\langle(a^{\dagger}a)^{2}\rangle-\langle
a^{\dagger}a\rangle^{2}.$ (4)

We emphasize that the analogue of (1, 3) for a fermion mode does not generally
hold. One heuristic reason for this is that only for the bosonic mode the mean
(photon) number can be arbitrary large.

###  II.2 Many modes, relations with noise and heating

####  II.2.1 Linear Heisenberg evolution

Importantly, (3) extends to the completely general $N$-mode situation, where
instead of (1) we write for initial $a_{i}=a_{i}(0)$ and final
$b_{i}=a_{i}(t)$ Heisenberg operators

$b_{i}={\sum}_{j=1}^{N}\left(S_{ij}a_{j}+R_{ij}a^{\dagger}_{j}\right)+f_{i},\quad
i=1,\dots,N,$ (5)

where $S_{ij}$, $R_{ij}$ and $f_{i}$ are c-numbers; cf. (1)). Write (5) in
block-matrix form:

$\displaystyle\begin{pmatrix}b\\\
b^{\dagger}\end{pmatrix}=E\begin{pmatrix}a\\\
a^{\dagger}\end{pmatrix}+\begin{pmatrix}f\\\ f^{*}\end{pmatrix},\quad
E=\begin{pmatrix}S&R\\\ R^{*}&S^{*}\end{pmatrix},$ (6)

where $a=(a_{1},\dots a_{N})^{T}$, $a^{\dagger}$, $b$, $b^{\dagger}$, $f$, and
$f^{*}$ are $N$-columns, and where T and ∗ denote (resp.) transposition and
complex conjugation; below ${}^{\dagger}=^{*T}$ will denote hermitean
conjugation.

Now commutation relations
$[b_{i},b_{j}^{\dagger}]=[a_{i},a_{j}^{\dagger}]=\delta_{ij}$, where
$\delta_{ij}$ is the Kronecker’s delta, and $[b_{i},b_{k}]=[a_{i},a_{k}]=0$
lead from (6) to (resp.):

$\displaystyle SS^{\dagger}-RR^{\dagger}=I,\quad SR^{T}=RS^{T},$ (7)

where $I$ is the $N\times N$ unit matrix. Eqs. (7) imply

$\displaystyle E^{-1}=\begin{pmatrix}S^{\dagger}&-R^{T}\\\
-R^{\dagger}&S^{T}\end{pmatrix}.$ (8)

The reasoning that led to (7) is now applied to (8), since the same
commutation relations hold. Then we get in addition to (7) the following new
relations:

$\displaystyle S^{\dagger}S-R^{T}R^{*}=I,\quad S^{\dagger}R=R^{T}S^{*}.$ (9)

####  II.2.2 The initial state

Now assume that the initial state $\rho$ of $N$ modes holds the following two
conditions:

$\displaystyle\langle a_{j}\rangle\equiv{\rm tr}(\rho a_{j})=0,$ (10)
$\displaystyle\langle a_{i}a_{j}\rangle\equiv{\rm tr}(\rho a_{j}a_{i})=0,$
(11)

where $i,j=1,\dots,N$. Two interesting examples of (10, 11) are as follows.
First, (10) can refer to initially independent modes in states with $\langle
a_{i}\rangle=0$. Then (11) holds automatically due to the independence:

$\displaystyle\langle a_{i}a_{j}\rangle=\langle a_{i}\rangle\langle
a_{j}\rangle=0.$ (12)

Second, we can consider diagonal states $\rho_{\rm diag}$ that read in the
Fock basis

$\displaystyle\rho_{\rm
diag}=\sum_{\nu_{1},...,\nu_{N}=0}^{\infty}r_{\nu_{1},...,\nu_{N}}|{\nu_{1},...,\nu_{N}}\rangle\langle{\nu_{1},...,\nu_{N}}|,$
(13) $\displaystyle
a_{i}^{\dagger}a_{i}|{\mu_{1},...,\mu_{N}}\rangle=\mu_{i}|{\mu_{1},...,\mu_{N}}\rangle,$
(14)

where (14) defines the Fock basis, and where (13) ensures conditions (10, 11).
It should be clear that neither independence nor diagonality is necessary for
the validity of (10, 11); e.g. non-diagonal state holding (10, 11) can be
easily constructed starting from (13). To be concise, we will refer to the
states $\rho$ satisfying (10, 11) as generalized diagonal states.

####  II.2.3 Increase of the mean photon number

Using (10, 11) together with the first equation in (9) we find that the change
of the total occupation number is non-negative:

$\displaystyle{\sum}_{i=1}^{N}\left(\langle
b^{\dagger}_{i}b_{i}\rangle-\langle
a^{\dagger}_{i}a_{i}\rangle\right)={\sum}_{i=1}^{N}|f_{i}|^{2}+{\sum}_{i,j=1}^{N}|R_{ij}|^{2}$
$\displaystyle+2{\sum}_{i=1}^{N}{\rm tr}\left(Y_{i}\rho
Y_{i}^{\dagger}\right)\geq 0,$ (15)

where we defined

$\displaystyle Y_{i}\equiv{\sum}_{k=1}^{N}R^{*}_{ik}a_{k}.$ (16)

When deducing (15), condition (10) was needed for nullifying terms $\propto
f_{i}\langle a_{k}\rangle$ in (15), while (11) was needed for nullifying terms
$\propto(R^{\dagger}S)_{kl}\langle a_{k}a_{l}\rangle$.

Eq. (5) can describe absorption (attenuation) of photons from a few selection
target modes, at the expense of their overall increase. For the particular
case of Gaussian initial states, (15) follows from the result of Ref.
Hovhannisyan et al. (2020) on the maximal work. Thus, according to (15) the
full mean photon number can only increase under linear evolution.

Where these additional photons come from? Answering this question is
contingent on realization of the linear transformation. For example, the
genesis of additional photons is relatively clear when the increase of the
mean number of photons is accompanied by an increase in the overall mean
energy; cf. (3). This energy increase comes from external sources that realize
the linear dynamics. In particular, this is the case when the $N$ modes start
their evolution from the overall vacuum state, because then the mean energy
can only increase. More generally, the relation between the mean energy
increase and the mean photon number increase in a linear dynamics is absent:
the latter is more general than the former; see (32) for clarification. In
such cases the genesis of additional photons should be prescribed to the
general fact that the mean photon number is not conserved within linear
dynamics.

####  II.2.4 Noise increase and heating

We emphasize that (15) can be interpreted as uncertainty increase. To this
end, let us note, for a mode with annihilation operator $a$, that $\langle
a^{\dagger}a\rangle$ characterizes the dispersion $\langle\Delta a^{2}\rangle$
of $a$ Caves (1982):

$\displaystyle\langle\Delta a^{2}\rangle\equiv\frac{1}{2}\langle
aa^{\dagger}+a^{\dagger}a\rangle-|\langle a\rangle|^{2}=\langle
a^{\dagger}a\rangle+\frac{1}{2}-|\langle a\rangle|^{2}~{}~{}~{}~{}~{}~{}$ (17)
$\displaystyle=\langle x^{2}\rangle-\langle x\rangle^{2}+\langle
y^{2}\rangle-\langle y\rangle^{2},\quad a=x+iy,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(18)

where $x=(a+a^{\dagger})/2$ and $y$ are Hermitian operators. Eq. (17) is the
definition of dispersion for non-hermitian $a$, while (18) shows how it can be
measured via its Hermitian components $x$ and $y$. Note from (17) that for
$\langle a\rangle=0$, the dispersion $\langle\Delta a^{2}\rangle$ reduces to
the mean photon number $\langle a^{\dagger}a\rangle$ 111This quantity also
controls the shot noise in photodetection Garrison and Chiao (2008)..

For considered initial states (10), we have $\langle a_{i}\rangle=\langle
b_{i}\rangle=0$, and then (15, 18) imply that also the sum of uncertainties
(17) increases together with the photon number:

$\displaystyle{\sum}_{i=1}^{N}\left(\langle\Delta
b^{2}_{i}\rangle-\langle\Delta
a^{2}_{i}\rangle\right)={\sum}_{i=1}^{N}\left(\langle
b^{\dagger}_{i}b_{i}\rangle-\langle a^{\dagger}_{i}a_{i}\rangle\right)\geq 0,$
(19)

i.e. as the mean photon number rises, so does the total dispersion. Eq. (19)
holds due to initial conditions (10, 11) and will be interpreted as heating.
Likewise, the decrease of both quantities in (19)—that is possible due to
nonlinear interactions—will mean cooling; see below.

We close this part by stressing that the relation between
${\sum}_{i=1}^{N}\left(\langle\Delta b^{2}_{i}\rangle-\langle\Delta
a^{2}_{i}\rangle\right)$ and ${\sum}_{i=1}^{N}\left(\langle
b^{\dagger}_{i}b_{i}\rangle-\langle a^{\dagger}_{i}a_{i}\rangle\right)$ is not
automatic. For example, linear dynamics under a particular condition $f_{i}=0$
in (5) will hold (15) under condition (11) only, i.e. (10) is now not needed.
And then if $\langle a_{i}\rangle\not=0$, then generically also $\langle
b_{i}\rangle\not=0$ and ${\sum}_{i=1}^{N}\left(\langle\Delta
b^{2}_{i}\rangle-\langle\Delta a^{2}_{i}\rangle\right)\geq 0$ does not hold,
though ${\sum}_{i=1}^{N}\left(\langle b^{\dagger}_{i}b_{i}\rangle-\langle
a^{\dagger}_{i}a_{i}\rangle\right)\geq 0$ still holds due to $f_{i}=0$.

##  III Entropic formulation of the second law for bosons

###  III.1 When Bose-entropy increases for a linear dynamics?

####  III.1.1 Definition of Bose-entropy

Eq. (15) shows that for initial conditions (10, 11) the total mean number of
photons can only increase. In the context of this unidirectional change it is
natural to ask whether one can find a suitable entropy function that also
increases under linear dynamics. As we show below, the answer to this question
is positive provided the initial states and the type of the linear dynamics
are restricted.

First of all, we need to define the entropy function: as always with the
unitary dynamics the von Neumann entropy $-{\rm tr}(\rho\ln\rho)$ (with $\rho$
being the density matrix) is not suitable for defining the second law, since
it is conserved. We need a more coarse-grained (i.e. less microscopic)
definition of entropy. A good choice is the time-dependent Bose entropy

$\displaystyle S(t)={\sum}_{k=1}^{N}s(n_{k}(t)),\quad n_{k}(t)\equiv\langle
a^{\dagger}_{k}(t)a_{k}(t)\rangle,$ (20) $\displaystyle
s(n_{k})\equiv(1+n_{k})\ln[1+n_{k}]-n_{k}\ln[n_{k}].$ (21)

Eq. (20) is deduced for an ideal Bose gas from the microcanonic distribution
Landau and Lifshitz (2013). If $s(n_{k})$ from (21) is maximized for a fixed
mean energy $\hbar\omega_{k}n_{k}$ of the mode $k$ with frequency $k$, one
obtains the thermal expression for the mean occupation (photon) number.
Indeed, making the Lagrange function $s(n_{k})-\beta\hbar\omega_{k}n_{k}$,
where $\beta$ is the Lagrange multiplier (inverse temperature) one obtains
$n_{k}=(e^{\beta\hbar\omega_{k}}-1)^{-1}$. Eq. (20) also increases in time
within kinetic equations for weakly interacting bosons; see Madeira et al.
(2020) for a recent discussion.

####  III.1.2 Increase of Bose-entropy

To study the behavior of $S$ in time for our situation, we need to add an
additional initial condition in (10, 11)

$\displaystyle\langle a^{\dagger}_{i}a_{j}\rangle=\delta_{ij}\langle
a^{\dagger}_{i}a_{i}\rangle,$ (22)

where (22) holds for examples (12, 13). Without (22), i.e. staying with (10,
11) only, we cannot express $n_{k}(t)$ via $n_{k}(0)$. Together with (22) this
task is possible from (5):

$\displaystyle n_{i}(t)=\langle
b^{\dagger}_{i}b_{i}\rangle={\sum}_{k=1}^{N}(|S_{ik}|^{2}+|R_{ik}|^{2})n_{k}(0)$
(23) $\displaystyle+{\sum}_{k=1}^{N}|R_{ik}|^{2}+{\sum}_{i=1}^{N}|f_{i}|^{2},$
(24)

where (7) and (9) imply

$\displaystyle{\sum}_{k=1}^{N}(|S_{ik}|^{2}+|R_{ik}|^{2})=1+2{\sum}_{k=1}^{N}|R_{ik}|^{2}\geq
1,$ (25)
$\displaystyle{\sum}_{i=1}^{N}(|S_{ik}|^{2}+|R_{ik}|^{2})=1+2{\sum}_{i=1}^{N}|R_{ik}|^{2}\geq
1.$ (26)

Let us assume that (consistently with (25, 26)) there exists a double
stochastic matrix $\Theta_{ik}$, i.e. a matrix holding

$\displaystyle\Theta_{ik}\geq
0,\quad{\sum}_{i=1}^{N}\Theta_{ik}=1,\quad{\sum}_{k=1}^{N}\Theta_{ik}=1,$ (27)

such that 222For the validity of (29) we in fact need instead of (27) a
seemingly weaker condition, where ${\sum}_{i=1}^{N}\Theta_{ik}=1$ in (29) is
replaced by ${\sum}_{i=1}^{N}\Theta_{ik}\geq 1$. However, this condition
together with ${\sum}_{k=1}^{N}\Theta_{ik}=1$ and $\Theta_{ik}\geq 0$ leads to
${\sum}_{i=1}^{N}\Theta_{ik}=1$.

$\displaystyle|S_{ik}|^{2}+|R_{ik}|^{2}\geq\Theta_{ik}.$ (28)

Matrices $|S_{ik}|^{2}+|R_{ik}|^{2}$ that satisfy (28) are called double-
superstochastic Marshall et al. (2011); Marshall and Olkin (1979). Once (27,
28) are assumed, the derivation of the second law in the Bose-entropic
formulation becomes straightforward from noting that $s(n_{k})$ from (21) is a
positive, increasing and concave function:

$\displaystyle
S(t)={\sum}_{i=1}^{N}s(n_{i}(t))\geq{\sum}_{i=1}^{N}s\left[{\sum}_{k=1}^{N}\Theta_{ik}n_{k}(0)\right]$
$\displaystyle\geq{\sum}_{i,k=1}^{N}\Theta_{ik}s[n_{k}(0)]={\sum}_{k=1}^{N}s[n_{k}(0)]=S[0].$
(29)

Thus initial conditions (10, 11, 22) and dynamic restriction (28) are
sufficient for the second law (29).

####  III.1.3 Validity of inequality (28)

Note that (28) trivially holds for $|R_{ik}|^{2}=0$. We emphasize that (28)
implies (25, 26), but the converse does not hold. To avoid confusions note
that ${\sum}_{k=1}^{N}(|S_{ik}|^{2}+|R_{ik}|^{2})\leq 1$ and
${\sum}_{i=1}^{N}(|S_{ik}|^{2}+|R_{ik}|^{2})\leq 1$ do imply
$|S_{ik}|^{2}+|R_{ik}|^{2}\leq\Theta_{ik}$ for some double-stochastic matrix
$\Theta_{ik}$ Marshall et al. (2011); Marshall and Olkin (1979).

Inequality (28) holds for $N=2$; see Appendix A which also discusses the
simplest counter-example of (28) for $N=3$. A constructive necessary and
sufficient condition for the validity of (28) was found in Cruse (1975):

$\displaystyle{\sum}_{i\in{\cal I},\,k\in{\cal
J}}(|S_{ik}|^{2}+|R_{ik}|^{2})\geq|{\cal I}|+|{\cal J}|-N,$ (30)

where (30) should hold for all subsets ${\cal I}$ and ${\cal J}$ of
$\\{1,...,N\\}$, and where $|{\cal I}|$ and $|{\cal J}|$ are the number of
elements in (resp.) ${\cal I}$ and ${\cal J}$. Conditions (30) are
straightforward to check at least for not very large $N$. The physical meaning
of (30) is that sufficiently small values of $(|S_{ik}|^{2}+|R_{ik}|^{2})$ are
to be excluded.

More general (but less constructive) sufficient dynamical conditions for (29)
can be stated as well. For example, whenever (28) does not hold, but still

$\displaystyle
s\left[{\sum}_{k=1}^{N}(|S_{ik}|^{2}+|R_{ik}|^{2})n_{k}(0)\right.$
$\displaystyle\left.+{\sum}_{k=1}^{N}|R_{ik}|^{2}+{\sum}_{i=1}^{N}|f_{i}|^{2}\right]$
$\displaystyle\geq{\sum}_{k=1}^{N}(|S_{ik}|^{2}-|R_{ik}|^{2})s[n_{k}(0)],$
(31)

holds for all $i$, we sum both parts of (31) over $i$, employ (26) and find
$S(t)\geq S(0)$.

###  III.2 Similarities and differences with the standard formulation of the
second law

We found two unidirectional relations inherent in linear dynamics for bosons:
inequality (15) states on a increase of the mean photon number, while (29) is
about the increase of the Bose-entropy. It is useful to compare these
relations with the standard (Thomson’s) formulation of the second law Lindblad
(2001); Allahverdyan and Nieuwenhuizen (2002): a unitary dynamics does not
decrease the mean energy of a quantum system that started its evolution from a
Gibbsian equilibrium (or at least passive) state. The unitary dynamics is
realized via time-dependent, cyclically changing Hamiltonian; the cyclic
condition is needed for ensuring that the initial and final Hamiltonians are
equal Lindblad (2001); Allahverdyan and Nieuwenhuizen (2002).

Similarities:

– Eqs. (15, 29) and Thomson’s formulation refer to unidirectional changes
inherent in a unitary evolution. All of them hold for specific initial states.

– For the single-mode situation (15) [i.e. (3)] refers to the basically same
quantity as the Thomson’s formulation, since the mean photon number is just
proportional to the mean energy.

– Eq. (15) relates to noise increase; cf. (17, 18). The same holds for the
entropic formulation (29) that refers to the Bose-entropy (20). Thomson’s
formulation has a similar bridge, since it also tells about the broadening of
the energy distribution in the final state as compared to the initial state.
This broadening is quantified by the entropy of the energy probability
distribution Lindblad (2001); Allahverdyan and Nieuwenhuizen (2002).

Differences:

$\bullet$ The second law holds for any unitary evolution, while (15) is
restricted to a linear evolution of boson modes. Inequality (29) assumes even
more restriction; see (28) and (22).

$\bullet$ The direct relation between the energy and photon number is broken
for the multimode situation, i.e. the analogue of (15) for energy does not
hold: the mean energy change

$\displaystyle{\sum}_{i=1}^{N}\omega_{i}\left(\langle
b^{\dagger}_{i}b_{i}\rangle-\langle a^{\dagger}_{i}a_{i}\rangle\right),$ (32)

need not have a definite sign for initial conditions (10, 11). For (32), the
derivation that led to (15) breaks down at the point when after the summation
over index $i$, one needs to employ the first equation in (9). The same holds
for (31): it does not apply to the mean energy. In other words, (31) states
that the Bose-entropy must increase without simultaneously increasing the mean
energy (or at least keeping it constant).

$\bullet$ Applicability domain: the second law demands equilibrium (e.g.
Gibbsian), or at least passive initial state Lindblad (2001); Allahverdyan and
Nieuwenhuizen (2002), while (10,11) and (2) allow initial states that need not
be equilibrium or passive; cf. (13). Recall that a passive state has a density
matrix $\rho$ that a non-increasing function of the Hamiltonian $H$ Lindblad
(2001); Allahverdyan and Nieuwenhuizen (2002). For a (Gibbsian) equilibrium
state this function is specific: $\rho=e^{-\beta H}/{\rm tr}e^{-\beta H}$ with
$\beta>0$ being the inverse temperature Lindblad (2001); Allahverdyan and
Nieuwenhuizen (2002). Thus, (15) is more general than the second law in the
context of initial states, but at the same time it is less general in the
context of dynamics, as it is restricted to linear evolution. Inequality (29)
assumes more restriction on the initial state; see (22).

##  IV Cooling two equilibrium modes

###  IV.1 Set-up

Once (19) is understood to define heating for linear dynamics with initial
conditions (10, 11), it is natural to ask whether non-linear processes can
cool, i.e. decrease the initial number of photons. To facilitate the
thermodynamic meaning of this question, we shall consider two initially
Gibbsian equilibrium bosonic modes at the same temperature $T$. Now a single
equilibrium mode cannot be cooled by any unitary (generally nonlinear)
operation, since the mean occupation number is proportional to the energy, and
the energy decrease for such a situation is prohibited by the second law.
However, two initially equilibrium modes at different frequencies can be
cooled, in terms of the mean full occupation number, via specific non-linear
interactions. Hence, we shall first determine the optimal cooling, and then
turn to non-optimal but feasible scenario from the viewpoint of experimentally
realizable nonlinear interactions.

Consider the initial state of two modes with frequencies $\omega_{1}$ and
$\omega_{2}$ at temperature $T$:

$\displaystyle\rho=\xi\,e^{-\beta\sum_{i=1}^{2}\omega_{i}\hat{n}_{i}},\quad\xi=(1-e^{-\beta\omega_{1}})(1-e^{-\beta\omega_{2}}),$
(33) $\displaystyle\hat{n}_{i}\equiv a^{\dagger}_{i}a_{i},\quad
i=1,2,\quad\hat{n}\equiv{\sum}_{i=1}^{2}a^{\dagger}_{i}a_{i},$ (34)

where $\hbar=1$, $\beta=1/(k_{\rm B}T)$ and ${\hat{n}}_{i}$ are the occupation
number operator for each mode. The two-mode system undergoes a unitary process
that aims at cooling:

$\displaystyle\rho(t)=U\rho\,U^{\dagger},\qquad UU^{\dagger}=1.$ (35)

###  IV.2 COP and efficiency

Besides targeting the mean occupation number, we characterize the cooling via
two efficiency-like quantities. Since $\rho$ in (33) is an equilibrium state,
the final average energy found from (35) is larger than the initial one, which
is the second law:

$\displaystyle{\sum}_{i=1}^{2}\omega_{i}\Delta n_{i}\geq 0,\quad\Delta
n_{i}\equiv\text{tr}(\rho[U^{\dagger}\hat{n}_{i}U-\hat{n}_{i}]).$ (36)

Eq. (36) defines the energy cost of cooling and it motivates the usual
definition of coefficient of performance (COP) Allahverdyan et al. (2010),
where the achieved cooling $-{\sum}_{i=1}^{2}\Delta n_{i}>0$ is divided over
the energy cost ${\sum}_{i=1}^{2}\omega_{i}\Delta n_{i}$.

Let us define the frequency ratio as

$\displaystyle\alpha\equiv\frac{\omega_{2}}{\omega_{1}}<1.$ (37)

We use the dimensionless COP (coefficient of performance) conventionally
defined as:

$K=-\frac{\Delta n_{1}+\Delta n_{2}}{\Delta n_{1}+\alpha\Delta n_{2}},$ (38)

where a larger $K$ means e.g. a better cooling with a smaller energy cost. In
(38) we took $\alpha<1$ without loss of generality. Hence, the fact of cooling
$-{\sum}_{i=1}^{2}\Delta n_{i}>0$ implies via (36) and $\alpha<1$

$\displaystyle 0\leq\alpha(-\Delta n_{2})\leq\Delta n_{1}\leq(-\Delta n_{2}).$
(39)

Now (39) motivates us to define $\Delta n_{1}-\Delta n_{2}=|\Delta
n_{1}|+|\Delta n_{2}|$ as the total number of occupation changes introduced in
the system. This is consistent with thinking about the cooling as photon
conversion: some amount of low energy photons ($\Delta n_{2}<0$) transform
into a smaller amount of higher energy photons ($\Delta n_{1}>0$). The sum of
low energy photons given and high energy photons received will be the total
number of occupation changes. Only a fraction $\eta$ of those lead to cooling:

$\eta=-\frac{\Delta n_{1}+\Delta n_{2}}{\Delta n_{1}-\Delta n_{2}}.$ (40)

We call $\eta$ the efficiency of cooling. It is similar to other quantum
efficiencies employed in optics Garrison and Chiao (2008); Walls and Milburn
(2007). Using (36, 39) we get a bound where temperatures are replaced by
frequencies:

$\eta\leq\frac{\Delta n_{1}+\Delta n_{2}}{\Delta n_{2}}\leq
1-\frac{\text{min}[\omega_{1},\omega_{2}]}{\text{max}[\omega_{1},\omega_{2}]},$
(41)

i.e. cooling is impossible for $\omega_{1}=\omega_{2}$. Note that (41) is more
similar to the Otto efficiency than to the Carnot efficiency of heat-engines
Allahverdyan et al. (2008).

###  IV.3  Optimal cooling

Given (35,34), we look for the unitary which minimizes the mean of $\hat{n}$
in the final state:

$\displaystyle U_{\rm opt}={\rm argmin}_{U}\,[\,{\rm tr}(U\rho
U^{\dagger}\hat{n})],$ (42)

Noting the eigenresolutions [cf. (33, 34)]

$\rho={\sum}_{k=0}^{\infty}r_{k}|r_{k}\rangle\langle
r_{k}|,\quad\hat{n}={\sum}_{l=0}^{\infty}\nu_{l}|\nu_{l}\rangle\langle\nu_{l}|,$
(43)

we get from (35, 43, 42)

$\displaystyle{\rm tr}(U\rho
U^{\dagger}\hat{n})={\sum}_{k,l=0}^{\infty}r_{k}\nu_{l}z_{kl},~{}z_{kl}=|\langle\nu_{l}|U|r_{k}\rangle|^{2},~{}~{}$
(44)

where

$\displaystyle{\sum}_{k}z_{kl}={\sum}_{l}z_{kl}=1,$ (45)

i.e. $z_{km}$ is a doubly stochastic matrix; cf. (27). Such matrices form a
compact convex set with vertices being permutation matrices Marshall and Olkin
(1979). As (44) is linear over $z_{km}$, it reaches the minimum value on the
vertices, i.e. on permutation matrices $z_{kl}$. This implies from (44) that
$U_{\rm opt}$ can be chosen as a permutation matrix.

Thus $U_{\rm opt}$ is a permutation matrix, and its form is seen from (44,
43):

$\displaystyle\text{min}_{U}[\,{\rm tr}(U\rho
U^{\dagger}{\hat{n}})\,]={\sum}_{k=0}^{\infty}\nu_{k}^{\uparrow}r_{k}^{\downarrow},$
(46)
$\displaystyle\nu_{1}^{\uparrow}\leq\nu_{2}^{\uparrow}\leq\nu_{3}^{\uparrow}...,\qquad
r_{1}^{\downarrow}\geq r_{2}^{\downarrow}\geq r_{3}^{\downarrow}...,$ (47)

where in (47) [cf. (43)] the ordered (anti-ordered) eigenvalues of $\hat{n}$
($\rho$) refer to the final state in (35). We visualize the orderings of
eigenvalues in the initial state (33, 34):

$\hat{n}$ | $0$ | $1$ | $2$ | $3$ | …
---|---|---|---|---|---
| $(0,0)$ |  $(0,1)$, $(1,0)$, |  $(0,2)$, $(1,1)$, $(2,0)$ |  $(0,3)$, $(1,2)$, $(2,1)$, $(3,0)$ | …
$\rho$ | $1$ |  $y^{\alpha}$, $y$ |  $y^{2\alpha}$, $y^{\alpha+1}$, $y^{2}$ |  $y^{3\alpha}$, $y^{2\alpha+1}$, $y^{\alpha+2}$, $y^{3}$ | …
(48)

where $y\equiv e^{-\beta\omega_{1}}$. The first, second and third row in (48)
show the eigenvalues of (resp.) $\hat{n}$, $(\hat{n}_{1},\hat{n}_{2})$ and
$\rho$, with the prefactor $\xi$ is omitted; cf. (33). The unitary process
(35, 46) permutes the eigenvalues of $\rho$. Using (48) one calculates
averages of $\hat{n}$ and $\hat{n}_{i}=a_{i}^{\dagger}a_{i}$:

$\begin{split}\langle\hat{n}\rangle&=\xi\left(1y^{\alpha}+1y+2y^{2\alpha}+2y^{\alpha+1}+2y^{2}+\dots\right),\\\
\langle\hat{n}_{1}\rangle&=\xi\left(0y^{\alpha}+1y+0y^{2\alpha}+1y^{\alpha+1}+2y^{2}+\dots\right),\\\
\langle\hat{n}_{2}\rangle&=\xi\left(1y^{\alpha}+0y+2y^{2\alpha}+1y^{\alpha+1}+0y^{2}+\dots\right).\end{split}$
(49)

The eigenvalues of $\rho$ in (48) are organized in columns. Whenever the
maximal element $y^{k\alpha}$ of $k$’th column is larger than the minimal
element $y^{l}$ of $l$’th column ($l<k$), we interchange them and achieve some
cooling. Formally, we should iterate till all elements in the third row are
arranged in descending order; cf. (47). Thus the optimal cooling increases the
probability of eigenstates of $\hat{n}$ with lower photon number. Note from
(48, 49) that we can interchange elements within each column without changing
$\Delta n$.

When $\Delta n$ is fixed, the descending order of $\rho$’s eigenvalues in the
final state yields simultaneously the minimum value of $\Delta n_{1}$ and the
maximum value of $\Delta n_{2}$. This is because the eigenvalues of
$\hat{n}_{1}$ ($\hat{n}_{2}$) in (48) are arranged in ascending (descending)
order. Eqs. (38, 40) show that thereby also $\eta$ and $K$ reach their maximal
values at the optimal $\Delta n$. The rule (49) stays intact and can be used
after permutations.

The exact calculation of (46) is out of reach, since $\rho$ has an infinite
number of eigenvalues. But we can develop a useful bound for it by focusing on
permutations between nearest-neighbour columns. Define from (37):

$\displaystyle m\equiv\lceil\frac{\alpha}{1-\alpha}\rceil,$ (50)

where $\lceil c\rceil$ is the smallest integer $\geq c$. Looking at (48) it is
seen that for $k\geq m$, the maximal element of the $(k+1)$’th column is
larger than the minimal element of $k$’th column. Permuting them will
contribute to $\Delta n$ calculated via (46). Likewise, for $k\geq m+2$, the
next to maximal element of the $(k+1)$’th column is larger than the next to
minimal element of $k$’th column. To visualize this situation consider a part
of (48) between columns $m+p$ and $m+p+1$ ($p\geq 0$):

(53) (54)

where we omitted the second row of (48). Continuing this logic, we see that a
new permutation appears for each even $p$, and that we can cover all nearest-
neighbor permutations. Hence a bound [cf. (33, 46)]:

$\displaystyle 0<-\Delta n_{\rm
opt}\equiv{\sum}_{k=0}^{\infty}(n_{k}r_{k}-n_{k}^{\uparrow}r_{k}^{\downarrow})$
$\displaystyle\geq\xi{\sum}_{l=0}^{\infty}y^{l(\alpha+1)}\times\sum_{k=m}^{\infty}(y^{(k+1)\alpha}-y^{k})$
$\displaystyle=\frac{(1-y)y^{\alpha(m+1)}-(1-y^{\alpha})y^{m}}{1-y^{\alpha+1}}.$
(55)

According to (55), cooling is possible for any $0\leq\alpha<1$, i.e. (55) is
positive and grows with $y=e^{-\beta\omega_{1}}$ changing from $0$ (at $y=0$)
to $\frac{1-\alpha}{1+\alpha}$ at $y=1$. Appendix B studies the optimal
cooling numerically. In particular, it shows numerical plots for the optimal
$K_{\rm opt}$ (COP) and $\eta_{\rm opt}$ (efficiency).

Now assume that $m$ given by (50) satisfies $m\gg 1$. Then the bound (55) gets
small, and becomes nearly exact, since the relative error between $\Delta
n_{\rm opt}$ and (55) scales as ${\cal O}(y^{2m})$. This estimate follows from
the contribution of next to nearest-neighbor permutations and is confirmed in
Appendix C. We report here the limiting values of $K$ and $\eta$ only, which
are obtained as described above [cf. (38, 40, 41)]:

$\displaystyle\alpha\to 1:\quad K_{\rm opt}\to\infty,\quad\eta_{\rm opt}\to
0,$ (56) $\displaystyle\alpha\to 0:\quad K_{\rm opt}\to\infty,\quad\eta_{\rm
opt}\to 1,$ (57)

where $\alpha=\omega_{2}/\omega_{1}\to 0$ in (57) is understood in the sense
of a large $\omega_{1}$ and a small $\omega_{2}$. It is also important to note
that both $\Delta n_{1}$ and $\Delta n_{2}$ are functions of $\alpha$ and in
the limit of $\alpha\to 1$ both tend to zero. In the last limits of (56, 57)
$\eta$ coincides with Otto bound. In both limits the energy costs of cooling
are negligible: $K_{\rm opt}\to\infty$. In the more general case of large
$\omega_{1}$ and fixed $\omega_{2}$, we studied $\eta_{\rm opt}$ and $K_{\rm
opt}$ in Appendix D.

##  V Feasible interaction Hamiltonian for cooling

How is a permutation unitary $U_{\rm opt}$ realized? This relates to one of
major questions of quantum control; see e.g. Wu et al. (2015). Any Hamiltonian
that is a polynomial of a fixed degree over $a_{1}$, $a_{1}^{\dagger}$,
$a_{2}$ and $a_{2}^{\dagger}$ can be realized via sufficiently many linear
operations plus a single-mode non-linearity Lloyd and Braunstein (1999).
However, realizing the permutation $U_{\rm opt}$ should be difficult in
practice, since it refers to a Hamiltonian that is a highly non-linear over
$a_{1}$, $a_{1}^{\dagger}$, $a_{2}$ and $a_{2}^{\dagger}$.

Now we focus on a feasible non-linear interaction and determine its cooling
ability. The feasiblity comes at a cost: now cooling will be possible mostly
next to nonlinear resonances: $\omega_{2}\gtrsim 2\omega_{1}$ or
$2\omega_{2}\lesssim\omega_{1}$. This will also connect to the Manley-Rowe
theorem, a known relation of nonlinear optics Weiss (1957); Landau and
Lifshitz (2013).

The simplest $\chi^{2}$ nonlinear interactions can be realized in an
anisotropic (e.g. crystalline) medium. Here the medium polarization $\vec{P}$
is quadratic in electric field $\vec{E}$ Landau and Lifshitz (2013); New
(2011); Hillery (2009); Walls and Milburn (2007):
$\vec{P}=\chi^{(1)}\vec{E}+\vec{E}\chi^{(2)}\vec{E}$, where $\chi^{(1)}$ and
$\chi^{(2)}$ are susceptibilities. Neglecting the polarization degree of
freedom for the electric field, its quantum operator representation is
$\vec{E}\to a^{\dagger}+a$ Hillery (2009); Walls and Milburn (2007). Hence,
the nonlinear interaction can be written as

$H_{I}=(a_{1}^{\dagger}+a_{1})(a_{2}^{\dagger}+a_{2})^{2}+(a_{1}^{\dagger}+a_{1})^{2}(a_{2}^{\dagger}+a_{2}),$
(58)

with the full Hamiltonian of the system being

$H=\omega_{1}a_{1}^{\dagger}a_{1}+\omega_{2}a_{2}^{\dagger}a_{2}+gH_{I}=H_{0}+gH_{I}$
(59)

where $g$ is the interaction constant.

Yet another scenario for (59) is realized in the optomechanics. In addition to
its applications in quantum technologies Aspelmeyer et al. (2014), this field
emerged as a potential basis for quantum gravity and foundations of quantum
mechanics Belenchia et al. (2016); Armata et al. (2017). In the optomechanical
setting, the interaction between a laser and a mechanical oscillator is such
that the resonance frequency $\omega_{1}(x)$ of the laser depends on the
position $x$ of the mechanical oscillator. Hence their joint Hamiltonian
reads: $H=\omega_{1}(x)a_{1}^{\dagger}a_{1}+\omega_{2}a_{2}^{\dagger}a_{2}$
Aspelmeyer et al. (2014). Here $a_{1}$ and $a_{2}$ are the annihilation
operators for (resp.) the laser and the mechanical oscillator. Keeping up to
the linear term of the Taylor expansion of $\omega_{1}(x)$ and using
$x=a_{2}^{\dagger}+a_{2}$ we get

$\displaystyle
H=\omega_{1}a_{1}^{\dagger}a_{1}+\omega_{2}a_{2}^{\dagger}a_{2}+(\partial_{x}\omega_{1})a_{1}^{\dagger}a_{1}(a_{2}^{\dagger}+a_{2}),$
(60)

which closely relates to (59).

To employ (58, 59) in (36) we introduce the free Heisenberg interaction
Hamiltonian $H_{I}(t)=e^{iH_{0}t}H_{I}e^{-iH_{0}t}$ and represent
$\rho(t)={\rm e}^{-itH}\rho\,{\rm e}^{itH}$ in (35) via chronological exponent
$\overleftarrow{\rm e}$:

$\displaystyle\rho(t)={\rm
e}^{-iH_{0}t}\widetilde{U}\rho\,\widetilde{U}^{\dagger}{\rm
e}^{iH_{0}t},~{}~{}\widetilde{U}=\overleftarrow{\rm
e}^{-i\int_{0}^{t}ds\,gH_{I}(s)}.$ (61)

Now expand $\widetilde{U}$ into Dyson series

$\displaystyle\widetilde{U}$ $\displaystyle=$ $\displaystyle
1-ig\int_{0}^{t}dsH_{I}(s)$ (62) $\displaystyle-$ $\displaystyle
g^{2}\int_{0}^{t}ds_{1}\int_{0}^{s_{1}}ds_{2}H_{I}(s_{1})H_{I}(s_{2})+...$

Using $e^{iH_{0}s}a_{k}e^{-iH_{0}s}=e^{-i\omega_{k}s}a_{k}$ ($k=1,2$) in
$H_{I}(t)$, one can show that the order of magnitude estimate of the $k$’th
term in (62) reads

$\displaystyle g^{k}\Omega^{-k}\,\sin^{k}\left({\Omega t}/{2}\right),$ (63)
$\displaystyle\Omega={\rm
min}[\omega_{1},\omega_{2},|2\omega_{1}-\omega_{2}|,|2\omega_{2}-\omega_{1}|].$
(64)

Thus, for a suitable $g$, $\omega_{1}$ and $\omega_{2}$ we can keep in (62)
the first three terms. Within this weak-coupling approximation we calculated
(36) in Appendix E showing that sufficiently large cooling $\Delta n<0$ is
possible only for

$\displaystyle\omega_{2}\gtrsim 2\omega_{1}\quad{\rm or}\quad
2\omega_{2}\lesssim\omega_{1},$ (65)

i.e. for two possible near-resonance conditions. Restricting ourselves with
the latter case $\alpha\equiv\omega_{2}/\omega_{1}\lesssim 0.5$ we note that
terms $a_{1}a_{2}^{\dagger 2}+a_{1}^{\dagger}a_{2}^{2}$ in (58) oscillate much
slower than other terms. Hence within the rotating wave approximation we can
take in (58):

$\displaystyle H_{I}\simeq\overline{H}_{I}\equiv a_{1}a_{2}^{\dagger
2}+a_{1}^{\dagger}a_{2}^{2}.$ (66)

The approximation is studied in Appendix E, where we also work out (58). Now
$\overline{H}_{I}$ in (66) leads to an exact operator conservation:

$\displaystyle 2\hat{n}_{1}+\hat{n}_{2}={\rm
const},~{}~{}\hat{n}_{k}=a^{\dagger}_{k}a_{k},~{}~{}k=1,2.$ (67)

This conservation is the Manley-Rowe theorem for the considered nonlinear
system Weiss (1957); Landau and Lifshitz (2013). The theorem does not
generally hold for the complete interaction Hamiltonian (58). However, the
cooling necessitates $\alpha\lesssim 0.5$ (or $\alpha\gtrsim 2$) and is
accompanied by an approximate conservation law (67) (or
$\hat{n}_{1}+2\hat{n}_{2}={\rm const}$). Using (66, 67) we get from (62, 61,
36) keeping there the first three terms only (the order of $g^{2}$):

$\displaystyle\Delta
n_{1}=\frac{8g^{2}\sin^{2}\left(\frac{(2\omega_{2}-\omega_{1})t}{2}\right)}{\left(2\omega_{2}-\omega_{1}\right){}^{2}}\frac{\left(e^{\beta\omega_{1}}-e^{2\beta\omega_{2}}\right)}{\left(e^{\beta\omega_{1}}-1\right)\left(e^{\beta\omega_{2}}-1\right){}^{2}},$
(68) $\displaystyle\Delta n_{2}=-2\Delta n_{1},\,\Delta n=-\Delta n_{1},$ (69)

Hence the cooling at $\alpha\lesssim 0.5$ is described via $\eta=\frac{1}{3}$
and $K=\frac{1}{1-2\alpha}$; cf. (38, 40). Once $\eta$ is finite and $K$ is
large, we achieve cooling with a small energy cost.

Eq. (68) shows that a sizable cooling is achieved for sufficiently long times,
because $\sin^{2}\left(\frac{(2\omega_{2}-\omega_{1})t}{2}\right)$ maximizes
for ${|2\omega_{2}-\omega_{1}|t}\sim{\pi}$, while $|2\omega_{2}-\omega_{1}|$
is small; cf. (65). This relation resembles the third law for the ordinary
(energy) cooling, though more efforts are needed for its systematic
investigation; e.g. we need a more complete understanding of the evolution
generated by (59).

##  VI Summary

Our starting point was that linear transformations on boson modes (linear
optics) increase the overall mean photon number, provided that the initial
state is (generalized) diagonal; see (10, 11) and (15). This unidirectional
relation refers to the linear evolution, but applies for a wider set of
initial states (10, 11) than the second law does. Its similarities and
differences with respect to the second law are discussed in section III.2. In
its full generality this relation is formulated for the first time, though the
literature was close to its formulation several times Caves (1982);
Hovhannisyan et al. (2020). Given that the lion’s share of boson dynamics is
linear, this general result will hold for a number of fields including optics,
phononics etc. Importantly, we show explicitly that relation (15) connects to
increasing the overall noise in the system (though its subsystems can get a
noise reduction, as e.g. happen in squeezing Garrison and Chiao (2008)). Hence
we interpret it as heating.

It is interesting to ask how specifically the increase (15) of the overall
mean photon number for initial states (10, 11) relates to the second law. To
answer this question, we studied the behavior of the Bose-entropy (20) for
linear dynamics and for the same class of initial states (10, 11). The Bose-
entropy is conditionally maximized at equilibrium, and it can change during
unitary evolution in contrast to the (fine-grained) von Neumann entropy. We
show in section III that for a subclass of linear evolution the Bose-entropy
(20) increases, and this increase also demands more restricted initial states
(10, 11, 22) than the validity of (15). A precise definition of this subclass
relates to certain non-trivial problems in linear algebra. We thus confirm
that for linear evolution the increase (15) of the overall mean photon number
is a more general unidirectional relation than the second law.

We show that the inverse of the heating in terms of the mean photon number
(i.e. cooling) is possible within nonlinear (inter-mode) interactions. The
cooling interpretation is not arbitrary and is characterized by efficiency and
coefficient of performance (COP). The former holds Otto’s bound of the heat-
engine efficiency (i.e. Carnot efficiency with temperatures replaced by
frequencies). For the COP we anticipated, but so far did not identify, a
general relation similar to Carnot’s bound for the refrigeration COP
Allahverdyan et al. (2010).

We studied feasible nonlinear processes (e.g. $\chi^{2}$ Landau and Lifshitz
(2013); New (2011); Hillery (2009); Walls and Milburn (2007)) on two modes
with different frequencies $\omega_{1}$ and $\omega_{2}$. Then the cooling in
terms of the mean photon number happens (mostly) in the vicinity of nonlinear
resonances. We also studied the optimal cooling, which is possible for any
$\omega_{1}\not=\omega_{2}$, but is demanding from the viewpoint of dynamic
realization.

###### Acknowledgements.

We are grateful to Karen Hovhannisyan for important remarks and to David
Petrosyan for discussions. This work was supported by SCS of Armenia, grant
No. 20TTAT-QTa003.

## References

  * Sheik-Bahae and Epstein (2007) M. Sheik-Bahae and R. I. Epstein, Nature Photonics 1, 693 (2007).
  * Abragam and Goldman (1978) A. Abragam and M. Goldman, Reports on Progress in Physics 41, 395 (1978).
  * Walls and Milburn (2007) D. F. Walls and G. J. Milburn, _Quantum optics_ (Springer Science & Business Media, 2007).
  * E. B. Stuart (1970) A. B. E. B. Stuart, B. Gal-Or, ed., _A Critical Review of Thermodynamics_ (Mono Book Corporation, New York, 1970).
  * Mahler (2014) G. Mahler, _Quantum thermodynamic processes:
Energy and information flow at the nanoscale_ (CRC Press, 2014).

  * Silva et al. (2016) R. Silva, G. Manzano, P. Skrzypczyk, and N. Brunner, Physical Review E 94, 032120 (2016).
  * Wilming and Gallego (2017) H. Wilming and R. Gallego, Physical Review X 7, 041033 (2017).
  * Clivaz et al. (2019) F. Clivaz, R. Silva, G. Haack, J. B. Brask, N. Brunner, and M. Huber, Physical review letters 123, 170605 (2019).
  * Freitas et al. (2018) N. Freitas, R. Gallego, L. Masanes, and J. P. Paz, in _Thermodynamics in the Quantum Regime_ (Springer, 2018), pp. 597–622.
  * Raeisi (2021) S. Raeisi, Physical Review A 103, 062424 (2021).
  * Gelbwaser-Klimovsky et al. (2015) D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki, Advances In Atomic, Molecular, and Optical Physics 64, 329 (2015).
  * Uzdin et al. (2015) R. Uzdin, A. Levy, and R. Kosloff, Physical Review X 5, 031044 (2015).
  * Liuzzo-Scorpo et al. (2016) P. Liuzzo-Scorpo, L. A. Correa, R. Schmidt, and G. Adesso, Entropy 18, 48 (2016).
  * Long and Liu (2016) R. Long and W. Liu, Physica A: Statistical Mechanics and its Applications 443, 14 (2016).
  * Gonzalez-Ayala et al. (2018) J. Gonzalez-Ayala, A. Medina, J. Roco, and A. C. Hernández, Physical Review E 97, 022139 (2018).
  * Singh et al. (2020) V. Singh, T. Pandit, and R. S. Johal, Physical Review E 101, 062121 (2020).
  * Raeisi and Mosca (2015) S. Raeisi and M. Mosca, Physical review letters 114, 100404 (2015).
  * Taranto et al. (2021) P. Taranto, F. Bakhshinezhad, A. Bluhm, R. Silva, N. Friis, M. P. Lock, G. Vitagliano, F. C. Binder, T. Debarba, E. Schwarzhans, et al., arXiv preprint arXiv:2106.05151 (2021).
  * Allahverdyan et al. (2011) A. E. Allahverdyan, K. V. Hovhannisyan, D. Janzing, and G. Mahler, Physical Review E 84, 041109 (2011).
  * Allahverdyan et al. (2010) A. E. Allahverdyan, K. Hovhannisyan, and G. Mahler, Physical Review E 81, 051129 (2010).
  * Weiss (1957) M. T. Weiss, Proc. IRE 45, 1012 (1957).
  * Landau and Lifshitz (2013) L. Landau and E. Lifshitz, _Electrodynamics of continuous media_ , vol. 8 (Elsevier, 2013).
  * Caves (1982) C. M. Caves, Physical Review D 26, 1817 (1982).
  * Garrison and Chiao (2008) J. Garrison and R. Chiao, _Quantum optics_ (OUP Oxford, 2008).
  * Hovhannisyan et al. (2020) K. V. Hovhannisyan, F. Barra, and A. Imparato, Physical Review Research 2, 033413 (2020).
  * Madeira et al. (2020) L. Madeira, A. D. García-Orozco, F. E. A. Dos Santos, and V. S. Bagnato, Entropy 22, 956 (2020).
  * Marshall et al. (2011) A. W. Marshall, I. Olkin, and B. C. Arnold, _Inequalities: Theory of Majorization and Its Applications_ (Springer Science & Business Media, 2011).
  * Marshall and Olkin (1979) A. W. Marshall and I. Olkin, _Inequalities: theory of majorization and its applications_ (Academic Press, NY, 1979).
  * Cruse (1975) A. B. Cruse, Linear Algebra and its Applications 12, 21 (1975).
  * Lindblad (2001) C. Lindblad, _Non-equilibrium entropy and irreversibility_ , vol. 5 (Springer Science & Business Media, 2001).
  * Allahverdyan and Nieuwenhuizen (2002) A. Allahverdyan and T. M. Nieuwenhuizen, Physica A: Statistical Mechanics and its Applications 305, 542 (2002).
  * Allahverdyan et al. (2008) A. E. Allahverdyan, R. S. Johal, and G. Mahler, Physical Review E 77, 041118 (2008).
  * Wu et al. (2015) R.-B. Wu, C. Brif, M. R. James, and H. Rabitz, Physical Review A 91, 042327 (2015).
  * Lloyd and Braunstein (1999) S. Lloyd and S. L. Braunstein, in _Quantum information with continuous variables_ (Springer, 1999), pp. 9–17.
  * New (2011) G. New, _Introduction to nonlinear optics_ (Cambridge University Press, 2011).
  * Hillery (2009) M. Hillery, arXiv preprint arXiv:0901.3439 (2009).
  * Aspelmeyer et al. (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * Belenchia et al. (2016) A. Belenchia, D. M. T. Benincasa, S. Liberati, F. Marin, F. Marino, and A. Ortolan, Phys. Rev. Lett. 116, 161303 (2016).
  * Armata et al. (2017) F. Armata, L. Latmiral, A. D. K. Plato, and M. S. Kim, Phys. Rev. A 96, 043824 (2017).

##  Appendix A  Examples and counter-examples for inequality (27)

We are given any $2\times 2$ matrix

$\displaystyle\begin{pmatrix}a&b\\\ c&d\end{pmatrix}.$ (70)

with non-negative elements and

$\displaystyle a+b\geq 1,\quad c+d\geq 1,\quad a+c\geq 1,\quad b+d\geq 1.$
(71)

Define

$\displaystyle c={\rm min}[a,b,c,d]$ (72)

and note that only $c<1$ is non-trivial, since otherwise (27) holds for (70)
and any double-stochastic matrix. Now

$\displaystyle\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\geq\begin{pmatrix}1-c&c\\\ c&1-c\end{pmatrix}=\Theta,$ (73)

where the latter matrix is double-stochastic. Thus (27) holds for $N=2$.

The simplest counter-examples for (27) at $N=3$ is the following matrix with
non-negative elements:

$\displaystyle\begin{pmatrix}a_{11}&a_{12}&a_{13}\\\ 0&a_{22}&a_{23}\\\
0&a_{32}&a_{33}\\\ \end{pmatrix},$ (74)

where $a_{21}=a_{31}=0$. We assume

$\displaystyle{\sum}_{i=1}^{3}a_{ik}\geq 1,\quad{\sum}_{k=1}^{3}a_{ik}\geq
1,\quad$ (75)

and additionally

$\displaystyle a_{22}+a_{32}<1.$ (76)

If a double-stochastic matrix holding (27) exists, then we have

$\displaystyle\begin{pmatrix}a_{11}&a_{12}&a_{13}\\\ 0&a_{22}&a_{23}\\\
0&a_{32}&a_{33}\\\ \end{pmatrix}\geq\begin{pmatrix}1&0&0\\\
0&\Theta_{22}&\Theta_{23}\\\ 0&\Theta_{32}&\Theta_{33}\\\ \end{pmatrix}.$ (77)

Now the latter matrix cannot be double-stochastic and hold (27), as the
condition (76) is violated, if we take in (30) ${\cal I}={2,3}$ and ${\cal
J}={1,2}$.

##  Appendix B  Numerical results for the optimal cooling

Recall our discussion after (48) of the main text. There we explained that the
optimal cooling—with respect to all involved quantities $\Delta n_{\rm opt}$
(photon number difference), $K_{\rm opt}$ (COP or coefficient of performance)
and $\eta_{\rm opt}$ (efficiency)—is achieved once all eigenvalues of the
final density matrix are arranged in the descending order; see the third row
in (48) of the main text. Numerically, this means that we need to take a
sufficiently long but a finite sequence of eigenvalues (starting from the
largest one) and ensure that the results are stable with respect to increasing
the length of this block.

Our numerical results are shown in Figs. 1, 2 and 3. First, recall that in
$K=-\frac{\Delta n_{1}+\Delta n_{2}}{\Delta n_{1}+\alpha\Delta n_{2}}$, the
achieved photon number decrease $\Delta n=\Delta n_{1}+\Delta n_{2}<0$ is
divided over the dimensionless energy cost $\Delta n_{1}+\alpha\Delta n_{2}$;
cf. (38) of the main text. It is seen from Fig. 1 that $K_{\rm opt}$ as a
function of $\alpha$ [cf. (37)] has (singular) local minima at points
$\alpha=\frac{k}{k+1}$ where $k\in\mathbb{N}$ is an integer. We checked that
these local minima of $K$ come mostly from the singular behavior of the energy
cost $\Delta n_{1}+\alpha\Delta n_{2}$; see Fig. 2. Now $\Delta n$ (not shown
in figures) shows weak singularities at those points $\alpha=\frac{k}{k+1}$,
but these singularities are much weaker than those of the energy cost $\Delta
n_{1}+\alpha\Delta n_{2}$.

The origin of these singularities for $K_{\rm opt}$ (and $\Delta
n_{1}+\alpha\Delta n_{2}$) can be clarified as follows. Recall that [cf. (50)]

$\displaystyle m\equiv\left\lceil\frac{\alpha}{1-\alpha}\right\rceil,$ (78)

refers to the the group of eigenvalues of the initial state $\rho$ starting
from which the eigenvalues of $\rho$ are not arranged in the descending order;
see the discussion after (48) of the main text. At points
$\alpha=\frac{k}{k+1}$ the index of the block from which the permutations
start undergoes a jump discontinuity of increasing by one.

Figure 1: The optimal coefficient of performance (COP) $K_{\rm opt}$ versus
$\alpha=\omega_{2}/\omega_{1}<1$ for the optimal cooling. Here
$y^{\alpha}=e^{-\beta\omega_{2}}=0.6$ and numerical calculations are done up
to the block number $300$; see (48) of the main text. Figure 2: The same as in
Fig. 1, but for the optimal dimensionless energy cost $\Delta
n_{1}+\alpha\Delta n_{2}$ versus $\alpha$ defined via (37). Figure 3: The same
as in Fig. 1, but for the efficiency $\eta_{\rm opt}$ versus $\alpha$.

Fig. 3 presents the numerical behavior of $\eta_{\rm opt}$ as a function of
$\alpha$. It is seen that $\eta_{\rm opt}$ also shows singularities at
$\alpha=\frac{k}{k+1}$, though these singularities are weaker than those for
$K_{\rm opt}$; cf. Fig. 1. In particular, these singularities do not change
the monotonous behavior of $\eta_{\rm opt}$ as a function of $\alpha$.

##  Appendix C  Asymptotic results for the optimal cooling: The limit
$\alpha\to 1$

Eq. (55) of the main text provides the nearest-neighbour approximation for
$\Delta n$. There we also indicated that (55) of the main text becomes close
to its exact value whenever $m$ defined via (78) is sufficiently large, or,
equivalently $\alpha\to 1$. The precision of this approximation relates to the
necessity of next-nearest-neighbour permutations. The largest value of $p$ in
(54) of the main text, where such permutations are necessary can be estimated
from the following diagram:

(81)

Now note from (78) that $y^{2m}<y^{(2m+2)\alpha}$, i.e. a next-nearest-
neighbor permutation is necessary. Hence the contribution from next-nearest-
neighbor permutation scales as ${\cal O}(y^{2m})$, and for $m\gg 1$ this is
smaller than what was retained in (55) of the main text. This estimate is
crude, since it did not account for permutations that already occurred (within
the nearest-neighbor approach) between the columns $2m$ and $2m+1$. However,
it is sufficient for our purposes. Indeed, Fig. 4 shows the relative error of
numerically exact calculation of $\Delta n$ and compares it with (55) of the
main text showing that it is well within the above bound ${\cal O}(y^{2m})$.

Figure 4: Here we calculated $\Delta n$ in numerically exact way by directly
arranging all eigenvalues of the final density matrix $\rho(t)$ in the
descending order; cf. (48) of the main text. We took
$y^{\alpha}=e^{-\beta\omega_{2}}=0.8$ and numerical calculations were done up
to block number 100 which is greater than $5m$ [cf. (78)] for all those
$\alpha$ values included in the graph. This quantity was denoted by $\Delta
n_{\rm exact}$. We denote via $\Delta n_{\rm nn}$ the neasrest neighbor
approximation given by (55) of the main text.
Blue curve: the relative error $\frac{|\Delta n_{\rm exact}-\Delta n_{\rm
nn}|}{\Delta n_{\rm nn}}$. Red curve: $|\frac{y^{2m(\alpha)}}{\Delta n_{\rm
nn}}|$. This curve is kinked, because so is $m(\alpha)$; see (78). It is seen
that the relative error is well within the announced range $|\frac{{\cal
O}(y^{2m(\alpha)})}{\Delta n_{\rm nn}}|$.

###  C.1  COP in the limit $\alpha\to 1$

For studying COP $K$, we can write the mean changes of $\hat{n}_{1}$ and
$\hat{n}_{2}$ in the approximation of nearest-neighbor permutations [cf. (54)
of the main text]:

$\begin{split}\Delta
n_{1}&=\xi{\sum}_{i=\tilde{n}}^{\infty}(y^{\alpha(i+1)}-y^{i}){\sum}_{j=0}^{\infty}y^{j(1+\alpha)}(i+j)\\\
\Delta
n_{2}&=-\xi{\sum}_{i=\tilde{n}}^{\infty}(y^{\alpha(i+1)}-y^{i}){\sum}_{j=0}^{\infty}y^{j(1+\alpha)}(i+j+1),\end{split}$
(82)

where $\xi=(1-e^{-\beta\omega_{1}})(1-e^{-\beta\omega_{2}})$ is the
normalization factor; cf. (33) of the main text.

Note that for obtaining (82) we do not make any permutation within columns
with the same eigenvalue of $\hat{n}$; cf. (48, 87) of the main text. Doing
such permutations will make the estimates in (82) closer to the minimal value
of $\Delta n_{1}$ and the maximal value of $\Delta n_{2}$ (both for a fixed
$\Delta n$). Hence (82) suffices for bounding $K$ from below:

$K\geq\frac{1}{(1-\alpha)(\frac{y^{1+\alpha}}{1-y^{1+\alpha}}+\frac{1}{1-y^{\alpha}})}\implies\lim_{\alpha\to
1}K\to\infty,$ (83)

which is also observed numerically.

Figure 5: Numerical results for the coefficient of performance (COP) $K$. Here
$y^{\alpha}=0.8$ and numerical calculations are done up to block number $5m$
for all those $\alpha$ values included in the graph.

##  Appendix D  Asymptotic results for the optimal cooling: The limit
$\alpha\to 0$

###  D.1  Error estimation

For $\alpha$ finite and sufficiently close to $0$, the action of an optimal
unitary results in (84)

$\begin{split}&\begin{tabular}[]{ c| c| c| c| c|
c|}$N$&$0$&$1$&$\dots$&$a^{\prime}$&$\dots$\\\ \hline\cr$R$&$1$&$y^{\alpha}$,
$y^{2\alpha}$&$\dots$&$\dots$, $y^{m_{1}\alpha}$ , $y$, $y^{(m_{1}+1)\alpha}$,
$\dots$&$\dots$\\\ \end{tabular}\\\ &\begin{tabular}[]{| c|
c}$a^{\prime}+i$&$\dots$\\\ \hline\cr$\dots$, $y^{(m_{2}-1)\alpha}$, $y^{2}$,
$y^{(m_{2}+1)\alpha}$ $y^{(m_{2}+2)\alpha}$, $\dots$&$\dots$\\\
\end{tabular},\end{split}$ (84)

where $m_{1}=\lfloor 1/\alpha\rfloor$, $m_{2}=\lfloor 2/\alpha\rfloor\geq
2m_{1}$, $a^{\prime}$ is determined from

$\displaystyle a^{\prime}(a^{\prime}+1)/2\leq
m_{1}\leq(a^{\prime}+1)(a^{\prime}+2)/2$ (85)

and $i$ from

$(a^{\prime}+i)(a^{\prime}+i+1)/2\leq
m_{2}\leq(a^{\prime}+i+1)(a^{\prime}+i+2)/2.$ (86)

We see that $i<3a^{\prime}$. Now we show that for the calculation of averages
of photon numbers one can use

$\begin{tabular}[]{ c| c| c| c| c| c}$N$&$0$&$1$&$2$&$3$&$\dots$\\\
\hline\cr$R$&$1$&$y^{\alpha}$, $y^{2\alpha}$&$y^{3\alpha}$, $y^{4\alpha}$,
$y^{5\alpha}$&$y^{6\alpha}$, $y^{7\alpha}$, $y^{8\alpha}$,
$y^{9\alpha}$&$\dots$\\\ \end{tabular}.$ (87)

instead of (84), as in the limit $\alpha\to 0$ corresponding error terms
vanish. We denote by $n_{1}^{(0)}$ and $n_{1}^{(*)}$ the average $\hat{n}_{1}$
calculated with (84) and (87) respectively, and by $\Delta n_{1e}$ the error
term $n_{1}^{(0)}-n_{1}^{(*)}$. Firstly, we write the contribution from the
$a^{\prime\text{th}}$ block to the error term

$\begin{split}\Delta
n^{a^{\prime}}_{1e}&=(y-y^{(m_{1}+1)\alpha})(m_{1}+1-u)\\\
&+(y^{(m_{1}+1)\alpha}-y^{(m_{1}+2)\alpha})(m_{1}+2-u)+\dots\\\
&+(y^{(u+a^{\prime}-1)\alpha}-y^{(u+a^{\prime})\alpha})a^{\prime},\end{split}$
(88)

where $u=a^{\prime}(a^{\prime}+1)/2$. $\Delta n^{a^{\prime}}_{1e}$ can be
estimated from above

$\begin{split}\Delta
n^{(a^{\prime})}_{1e}&\leq(y-y^{(m_{1}+1)\alpha})a^{\prime}+(y^{(m_{1}+1)\alpha}-y^{(m_{1}+2)\alpha})a^{\prime}+\\\
&+\dots+(y^{(u+a^{\prime}-1)\alpha}-y^{(u+a^{\prime})\alpha})a^{\prime}=\\\
&=(y-y^{(u+a^{\prime})\alpha})a^{\prime}.\end{split}$ (89)

Similarly, one can estimate the contribution from $(a^{\prime}+1)^{\text{th}}$
block

$\Delta
n^{(a^{\prime}+1)}_{1e}\leq(y^{(u+a^{\prime})\alpha}-y^{(u+2a^{\prime}+2)\alpha})(a^{\prime}+1).$
(90)

Summing up all contributions we get

$\begin{split}&\Delta
n_{1e}\leq\big{[}y+y^{(u+a^{\prime})\alpha}+y^{(u+2a^{\prime}+2)\alpha}+y^{(u+3a^{\prime}+5)\alpha}+\dots\big{]}\\\
&+\big{[}(y^{2}-y^{(m_{2}+1)\alpha})(a^{\prime}+i)+(y^{3}-y^{(m_{3}+1)\alpha})(a^{\prime}+i^{\prime})\\\
&+(y^{4}-y^{(m_{4}+1)\alpha})(a^{\prime}+i^{\prime\prime})+\dots\big{]},\end{split}$
(91)

where $m_{3}=\lfloor 3/\alpha\rfloor\geq 3m_{1}$, $m_{4}=\lfloor
4/\alpha\rfloor\geq 4m_{1}$. $i^{\prime}$ and $i^{\prime\prime}$ in (91) are
determined from conditions similar to (86):

$\begin{split}&\frac{(a^{\prime}+i^{\prime})(a^{\prime}+i^{\prime}+1)}{2}\leq
m_{3}\leq\frac{(a^{\prime}+i^{\prime}+1)(a^{\prime}+i^{\prime}+2)}{2},\\\
&\frac{(a^{{}^{\prime}}+i^{{}^{\prime\prime}})(a^{{}^{\prime}}+i^{{}^{\prime\prime}}+1)}{2}\leq
m_{4}\leq\frac{(a^{{}^{\prime}}+i^{{}^{\prime\prime}}+1)(a^{{}^{\prime}}+i^{{}^{\prime\prime}}+2)}{2}\end{split}$
(92)

and result in $i^{\prime}<4a^{\prime}$ and $i^{\prime\prime}<5a^{\prime}$ Now,
we can estimate (91) further

$\begin{split}&\Delta n_{1e}\leq\big{[}y+y^{2}+y^{3}+y^{4}+\dots\big{]}\\\
&+\big{[}y^{2}3a^{\prime}+y^{3}4a^{\prime}+y^{4}5a^{\prime}+\dots\big{]}\leq\\\
&\leq\frac{y(1+y)}{1-y}+a^{\prime}\frac{y\ln{y}}{(1-y)^{2}}.\end{split}$ (93)

Now, note that in the limit $\alpha\to 0$, which is $\omega_{1}\to\infty$,
$a^{\prime}$ goes to infinity as $\sqrt{\frac{\omega_{1}}{\omega_{2}}}$. As,
$0\leq\Delta n_{1e}\leq\frac{y(1+y)}{1-y}+a^{\prime}\frac{y\ln{y}}{(1-y)^{2}}$
we conclude that (remember that $y=e^{-\omega_{1}\beta}$)

$\lim_{\alpha\to 0}\Delta n_{1e}\to 0.$ (94)

###  D.2  Asymptotic expressions and their integral representations

In our further calculations we use (87). Using the same procedure as in (49),
we find from (87) the following expressions for $\Delta n_{1,2}$
($\epsilon=\alpha\beta\omega_{1}$)

$\begin{split}&\Delta
n_{1}=\xi\sum_{a=0}^{\infty}e^{-\epsilon\frac{a(a+1)}{2}}\sum_{b=0}^{a}e^{-\epsilon
b}b-\xi_{1}\sum_{a}ae^{-\beta\omega_{1}a}\\\
&=\xi\sum_{a=0}^{\infty}\frac{e^{2\epsilon}}{(e^{\epsilon}-1)^{2}}e^{-\left(\frac{a}{2}+1\right)(a+1)\epsilon}\left(a\left(e^{-\epsilon}-1\right)+e^{a\epsilon}-1\right)\\\
&-\frac{1}{e^{\beta\omega_{1}}-1},\\\ &\Delta
n_{2}=\xi\sum_{a=0}^{\infty}e^{-\epsilon\frac{a(a+1)}{2}}\sum_{b=0}^{a}e^{-\epsilon(a-b)}b-\xi_{2}\sum_{a}ae^{-\epsilon
a}\\\
&=\xi\sum_{a=0}^{\infty}\frac{e^{\epsilon}}{(e^{\epsilon}-1)^{2}}e^{-\frac{a(a+1)}{2}\epsilon}(a(e^{\epsilon}-1)+e^{-a\epsilon}-1)\\\
&-\frac{1}{e^{\epsilon}-1},\end{split}$ (95)

where we defined $\xi_{1}=1-e^{-\beta\omega_{1}}$ and
$\xi_{2}=1-e^{-\beta\omega_{2}}$; hence $\xi=\xi_{1}\xi_{2}$. Before studying
(95) numerically, we apply Hubbard-Stratonovich transformation:

$\displaystyle
e^{-\frac{a^{2}}{2}\epsilon}=\sqrt{\frac{1}{2\pi\epsilon}}\int_{-\infty}^{\infty}dve^{-\frac{v^{2}}{2\epsilon}-iav},$
(96)

for faster and more accurate calculations:

$\begin{split}\\\
&\sum_{a=0}^{\infty}e^{-\left(\frac{a}{2}+1\right)(a+1)\epsilon}\left(a\left(e^{-\epsilon}-1\right)+e^{a\epsilon}-1\right)\\\
&=\sqrt{\frac{1}{2\pi\epsilon}}e^{-\epsilon}\int_{-\infty}^{\infty}dve^{-\frac{v^{2}}{2\epsilon}}\bigg{(}(e^{-\epsilon}-1)\frac{e^{-(iv+\frac{3}{2}\epsilon)}}{(1-e^{-(iv+\frac{3}{2}\epsilon)})^{2}}\\\
&+\frac{1}{(1-e^{-(iv+\frac{1}{2}\epsilon)})}-\frac{1}{(1-e^{-(iv+\frac{3}{2}\epsilon)})}\bigg{)},\end{split}$
(97) $\begin{split}\\\
&\sum_{a=0}^{\infty}e^{-\frac{a(a+1)}{2}\epsilon}(a(e^{\epsilon}-1)+e^{-a\epsilon}-1)\\\
&=\sqrt{\frac{1}{2\pi\epsilon}}\int_{-\infty}^{\infty}dve^{-\frac{v^{2}}{2\epsilon}}\bigg{(}(e^{\epsilon}-1)\frac{e^{-(iv+\frac{1}{2}\epsilon)}}{(1-e^{-(iv+\frac{1}{2}\epsilon)})^{2}}\\\
&+\frac{1}{(1-e^{-(iv+\frac{3}{2}\epsilon)})}-\frac{1}{(1-e^{-(iv+\frac{1}{2}\epsilon)})}\bigg{)}.\end{split}$
(98)

The results of numerical calculations for $\eta$ and $K$ are depicted in Fig.
6 and Fig. 7. As seen from figures, $K\to\infty$ and $\eta\to 1$ in the limit
$\epsilon\ll 1$. Below, we show analytically, that indeed, $K$ and $\eta$
reach these limits.

Figure 6: Numerically calculated efficiency ($\eta$) versus $\epsilon$ in the
limiting case $\omega_{1}\gg\omega_{2}$ using (97) and (98). Here, we set
$\beta\omega_{1}=10$ and the smallest value of $\epsilon$ is 0.007.  Figure 7:
The same as in Fig. 6 but for COP (K).

###  D.3  Asymptotic results via the Euler-Maclaurin formula

To study the asymptotics of $\eta$ and $K$ in the limit $\epsilon\ll 1$ we
apply the Euler-Maclaurin formula for the sums in (95)

$\begin{split}&\sum_{a=0}^{\infty}\frac{e^{2\epsilon}}{(e^{\epsilon}-1)^{2}}e^{-\left(\frac{a}{2}+1\right)(a+1)\epsilon}\left(a\left(e^{-\epsilon}-1\right)+e^{a\epsilon}-1\right)\\\
&\equiv\sum_{a=0}^{\infty}f_{1}(a)=\int_{0}^{\infty}dxf_{1}(x)+I_{1}=S_{1}+I_{1}\end{split}$
(99)

where

$\begin{split}I_{1}&=\frac{f_{1}(0)+f_{1}(\infty)}{2}\\\
&+\sum_{k=1}^{\left\lfloor\frac{p}{2}\right\rfloor}\frac{B_{2k}}{(2k)!}\left(f_{1}^{(2k-1)}(\infty)-f_{1}^{(2k-1)}(0)\right)+R_{p},\end{split}$
(100)

$R_{p}\leq\frac{2\bm{\zeta}(p)}{(2\pi)^{p}}\int_{0}^{\infty}dx|f_{1}^{(p)}(x)|,$
(101)

$B_{2k}$ are Bernoulli numbers, $\bm{\zeta}(p)$ is the Riemann’s zeta function
and $f^{(p)}(x)$ is the $p^{\text{th}}$ order differential. $p$ in (99) takes
different integer values $p\geq 2$ and we use $p=2$, because this is the
simplest case amenable to estimates. Similarly,

$\begin{split}&\sum_{a=0}^{\infty}\frac{e^{\epsilon}}{(e^{\epsilon}-1)^{2}}e^{-\frac{a(a+1)}{2}\epsilon}(a(e^{\epsilon}-1)+e^{-a\epsilon}-1)\equiv\sum_{a=0}^{\infty}f_{2}(a)\\\
&=\int_{0}^{\infty}dxf_{2}(x)+I_{2}=S_{2}+I_{2},\end{split}$ (102)

The leading diverging terms in (99) and (102) when $\epsilon\to 0$ are $S_{1}$
and $S_{2}$ and we omit $I_{1}$ and $I_{2}$. Using (99) and (102) for the
efficiency and COP we get the following relations

$\begin{split}K&\approx-\frac{\xi S_{1}-n_{1i}+\xi S_{2}-n_{2i}}{\xi
S_{1}-n_{1i}+\alpha(\xi S_{2}-n_{2i})},\\\ \eta&\approx-\frac{\xi
S_{1}-n_{1i}+\xi S_{2}-n_{2i}}{\xi S_{1}-n_{1i}-\xi S_{2}+n_{2i}},\end{split}$
(103)

where $n_{1i}$ and $n_{2i}$ are initial average occupation numbers. The limits
$\lim_{\alpha\to 0}\xi S_{1,2}/n_{2i}$ can be studied analytically, and we get

$\lim_{\alpha\to 0}\xi S_{1,2}/n_{2i}=0.$

Thus, for the $K_{\rm opt}$ and $\eta_{\rm opt}$ we obtain

$\displaystyle K_{\rm opt}\to\infty,\quad\eta_{\rm opt}\to 1.$ (104)

##  Appendix E  Perturbative treatment of the full nonlinear Hamiltonian

Let us return to the full—i.e. without the rotating-wave
approximation—nonlinear Hamiltonian given by (58) of the main text:

$\displaystyle
H_{I}=(a_{1}^{\dagger}+a_{1})(a_{2}^{\dagger}+a_{2})^{2}+(a_{1}^{\dagger}+a_{1})^{2}(a_{2}^{\dagger}+a_{2}).$
(105)

See (59) of the main text for the complete Hamiltonian. Here we shall employ
(105) in the second-order of Dyson’s series given by (62) of the main text;
see in this context (61) of the main text. For simplicity we shall scale out
the factor $\beta$, i.e. we denote $\beta g\to g$,
$\beta\omega_{1,2}\to\omega_{1,2}$ and $t/\beta\to t$.

Using (62, 61) of the main text we get

$\begin{split}&{\rm
tr}(\rho(t)\hat{n}-\rho(0)\hat{n})=\mathcal{O}(g^{3})+g^{2}\times\\\
&\times\bigg{[}\int_{0}^{t}dsH_{I}(s)\hat{n}\int_{0}^{t}dsH_{I}(s)-\\\
&-\int_{0}^{t}ds_{1}\int_{0}^{s_{1}}ds_{2}H_{I}(s_{1})H_{I}(s_{2})\hat{n}-\\\
&-\hat{n}\int_{0}^{t}ds_{1}\int_{0}^{s_{1}}ds_{2}H_{I}(s_{2})H_{I}(s_{1})\bigg{]}.\end{split}$
(106)

Formally the same equation holds for $\hat{n}_{k}=a_{k}^{\dagger}a_{k}$, where
$k=1,2$ and $\hat{n}=\hat{n}_{1}+\hat{n}_{2}$.

Substituting (105) into (LABEL:eq:dysoncooling) we get

$\begin{split}&\Delta n_{k}={\rm
tr}\left(\rho(t)\hat{n}_{k}-\rho(0)\hat{n}_{k}\right)=\\\
=&g^{2}\bigg{[}A_{k}\Phi(\omega_{1}+2\omega_{2})+B_{k}\Phi(\omega_{1}-2\omega_{2})+C_{k}\Phi(\omega_{1})\bigg{]}+\\\
+&g^{2}\bigg{[}D_{k}\Phi(\omega_{2}+2\omega_{1})+E_{k}\Phi(\omega_{2}-2\omega_{1})+F_{k}\Phi(\omega_{2})\bigg{]},\end{split}$
(107)

where $k=1,2$,

$\Phi(x)\equiv\frac{4\sin^{2}{(\frac{1}{2}xt)}}{x^{2}},$ (108)
$\begin{split}&A_{1}=\frac{2\left(e^{\omega_{1}+2\omega_{2}}-1\right)}{\left(e^{\omega_{1}}-1\right)\left(e^{\omega_{2}}-1\right)^{2}},\quad
A_{2}=\frac{4\left(e^{\omega_{1}+2\omega_{2}}-1\right)}{\left(e^{\omega_{1}}-1\right)\left(e^{\omega_{2}}-1\right)^{2}},\\\
&B_{1}=\frac{2\left(e^{\omega_{1}}-e^{2\omega_{2}}\right)}{\left(e^{\omega_{1}}-1\right)\left(e^{\omega_{2}}-1\right)^{2}},\quad
B_{2}=\frac{-4\left(e^{\omega_{1}}-e^{2\omega_{2}}\right)}{\left(e^{\omega_{1}}-1\right)\left(e^{\omega_{2}}-1\right)^{2}},\\\
&C_{1}=\frac{4e^{\omega_{2}}}{\left(e^{\omega_{2}}-1\right)^{2}},\quad\qquad\qquad
C_{2}=0.\end{split}$ (109)

Now $D_{1}$, $E_{1}$ and $F_{1}$ are obtained from (resp.) $A_{2}$, $B_{2}$
and $C_{2}$ upon swapping $\omega_{1}$ and $\omega_{2}$. Likewise, $D_{2}$,
$E_{2}$ and $F_{2}$ are obtained from (resp.) $A_{1}$, $B_{1}$ and $C_{1}$
upon swapping $\omega_{1}$ and $\omega_{2}$.

For a representative pair of frequencies $\omega_{1}$ and $\omega_{2}$, Fig. 8
demonstrates to which extent $\Delta n=\Delta n_{1}+\Delta n_{2}$ calculated
via (107) predicts cooling, i.e. $\Delta n<0$. As announced in the main text,
cooling happens in near-resonance conditions $\omega_{2}\gtrsim 2\omega_{1}$
or $2\omega_{2}\lesssim\omega_{1}$, which is seen in Fig. 8; see also Fig. 9
for additional information.

Now the essence of rotating-wave approximation in (107) is that e.g. for
$|2\omega_{2}-\omega_{1}|\ll{\rm
min}[\omega_{1},\omega_{2},|2\omega_{1}-\omega_{2}|]$, we can take
$\Phi(\omega_{1}-2\omega_{2})$ in (107) much larger than other terms. This
reverts to (69) of the main text.

Figure 8: The photon number difference $\Delta n$ obtained from (107) and
(109) for $\omega_{1}=0.35$ and $t=10\pi$, where
$\alpha=\omega_{2}/\omega_{1}$. It is seen that near the resonating
frequencies $\alpha\lesssim 0.5$ and $\alpha\gtrsim 2$ the interaction
Hamiltonian (105) results in cooling. We see that $\Delta n>0$ (no cooling)
for other values of $\alpha$. Figure 9: The photon number difference $\Delta
n$ obtained from (107) and (109) for $\omega_{1}=0.6$ and $t=6\pi$, where
$\alpha=\omega_{2}/\omega_{1}$. We see that at the non resonance point
$\alpha\approx 2.8$ there are small values of $\Delta n<0$ (cooling). Hence,
although very small, cooling can also be achieved far from the resonance
frequencies.

###  E.1  Estimation of the higher-order terms in Dyson’s series

Using (62) one can show that the terms $\mathcal{O}(g^{l})$ in Dyson’s series
(cf. (LABEL:eq:dysoncooling) of the main text are based on the following
structure:

$\begin{split}&\prod_{i=1}^{k}g\int_{0}^{s_{i-1}}ds_{i}H_{I}(s_{i})\times\rho(0)\times\\\
&\prod_{i=1}^{k^{\prime}}g\int_{0}^{s_{i-1}}ds_{i}H_{I}(s_{k^{\prime}-i+1})\hat{n},\end{split}$
(110)

where $k+k^{\prime}=l$ and $s_{0}=t$. To get from (110) the term
$\mathcal{O}(g^{l})$ in Dyson’s series we should take the trace of (110) and
sum it as $\sum_{k=1,\,k^{\prime}=1;k+k=l}^{l}$.

To study (110), let us take its leftmost multiplier

$\begin{split}&\prod_{i=1}^{k}g\int_{0}^{s_{i-1}}ds_{i}H_{I}(s_{i})=\\\
=&\sum_{\alpha_{1}=1}^{8}...\sum_{\alpha_{k}=1}^{8}\prod_{i=1}^{k}g\int_{0}^{s_{i-1}}ds_{i}h_{\alpha_{i}}(s_{i}).\end{split}$
(111)

Here $\\{h_{i}\\}_{i=1}^{8}$ is the set of all monomials in the interaction
Hamiltonian (59):

$\begin{split}\\{h_{i}\\}_{i=1}^{8}=\\{&a_{1}a_{2}^{2},a_{1}^{\dagger}a_{2}^{2},a_{1}a_{2}^{\dagger}{}^{2},a_{1}^{\dagger}a_{2}^{\dagger}{}^{2},\\\
&a_{1}^{2}a_{2},a_{1}^{\dagger}{}^{2}a_{2},a_{1}^{2}a_{2}^{\dagger},a_{1}^{\dagger}{}^{2}a_{2}^{\dagger}\\}.\end{split}$
(112)

Let us also define the frequency set $\\{W_{i}\\}_{i=1}^{8}$

$\begin{split}\\{W_{i}\\}_{i=1}^{8}=\\{&\omega_{1}+2\omega_{2},-\omega_{1}+2\omega_{2},\omega_{1}-2\omega_{2},-\omega_{1}-2\omega_{2},\\\
&2\omega_{1}+\omega_{2},-2\omega_{1}+\omega_{2},2\omega_{1}-\omega_{2},-2\omega_{1}-\omega_{2}\\}.\end{split}$
(113)

Keeping in mind the equation $a_{1,2}(s)=a_{1,2}e^{-is\omega_{1,2}}$ let us
take one term from the sum (111) corresponding to some
$\alpha_{1}...\alpha_{k}$:

$\begin{split}&\prod_{i=1}^{k}g\int_{0}^{s_{i-1}}ds_{i}h_{\alpha_{i}}(s_{i})=\prod_{i=1}^{k}gh_{\alpha_{i}}\times\\\
&\times\int_{0}^{s_{i-1}}ds_{i}\exp{\bigg{[}-is_{i}W_{\alpha_{i}}\bigg{]}}=\\\
&=\prod_{i=1}^{k}gh_{\alpha_{i}}\frac{1}{k!}\prod_{i=1}^{k}\int_{0}^{t}ds_{i}\exp{\bigg{[}-is_{i}W_{\alpha_{i}}\bigg{]}}.\end{split}$
(114)

The last step uses the fact that we have $k!$ ways to order $k$ different
items and that after taking the operator part out of the integration we get
integration of complex valued functions which do not change with ordering.
Similarly for the rightmost multiplier of (110):

$\prod_{i=0}^{k^{\prime}-1}gh_{\alpha^{\prime}_{k^{\prime}-i}}\frac{1}{k^{\prime}!}\prod_{i=1}^{k^{\prime}}\int_{0}^{t}ds_{i}\exp{\bigg{[}-is_{i}W_{\alpha^{\prime}_{i}}\bigg{]}}.$
(115)

Straightforward calculation shows that the integral terms in (114) and (115)
result in

$\prod_{i=1}^{k}\frac{(1-e^{-iW_{\alpha_{i}}t})}{(iW_{\alpha_{i}})},\;\prod_{i=1}^{k^{\prime}}\frac{(1-e^{-iW_{\alpha^{\prime}_{i}}t})}{(iW_{\alpha^{\prime}_{i}})}$
(116)

Now we can write (110) as

$\begin{split}&\sum_{\alpha_{1}=1}^{8}...\sum_{\alpha_{k}=1}^{8}\sum_{\alpha^{\prime}_{1}=1}^{8}...\sum_{\alpha^{\prime}_{k^{\prime}}=1}^{8}\frac{g^{k}}{k!}\frac{g^{k^{\prime}}}{k^{\prime}!}\times\\\
&\times\bigg{[}\prod_{i=1}^{k}h_{\alpha_{i}}\prod_{i=0}^{k^{\prime}-1}h_{\alpha^{\prime}_{k^{\prime}-i}}\times\\\
&\times\prod_{i=1}^{k}\frac{(1-e^{-iW_{\alpha_{i}}t})}{(iW_{\alpha_{i}})}\prod_{i=1}^{k^{\prime}}\frac{(1-e^{-iW_{\alpha^{\prime}_{i}}t})}{(iW_{\alpha^{\prime}_{i}})}\bigg{]}.\end{split}$
(117)

and the equation for $\Delta n_{1,2}$ will be

$\begin{split}&\Delta
n_{1,2}=\sum_{k=0}^{\infty}\sum_{k^{\prime}=0}^{\infty}\delta_{k+k^{\prime}}^{0}\sum_{\alpha_{1}=1}^{8}...\sum_{\alpha_{k}=1}^{8}\sum_{\alpha^{\prime}_{1}=1}^{8}...\sum_{\alpha^{\prime}_{k^{\prime}}=1}^{8}\frac{g^{k}}{k!}\frac{g^{k^{\prime}}}{k^{\prime}!}\times\\\
&\times\prod_{i=1}^{k}\frac{(1-e^{-iW_{\alpha_{i}}t})}{(iW_{\alpha_{i}})}\prod_{i=1}^{k^{\prime}}\frac{(1-e^{-iW_{\alpha^{\prime}_{i}}t})}{(iW_{\alpha^{\prime}_{i}})}\times\\\
&\times{\rm
tr}\bigg{(}\prod_{i=1}^{k}h_{\alpha_{i}}\rho(0)\prod_{i=0}^{k^{\prime}-1}h_{\alpha^{\prime}_{k^{\prime}-i}}\,\hat{n}_{1,2}\bigg{)}.\end{split}$
(118)

Here the sum
$\sum_{\alpha_{1}=1}^{8}...\sum_{\alpha_{k}=1}^{8}\sum_{\alpha^{\prime}_{1}=1}^{8}...\sum_{\alpha^{\prime}_{k^{\prime}}=1}^{8}$
will have $8^{l}$ elements for any $l$, so the amount of terms of order
$\mathcal{O}(g^{l})$ is $8^{l}(l+1)$. This may put doubt in the claim that the
higher order $\mathcal{O}(g^{l})$ terms of $\Delta n$ can be neglected.
However, we believe that it can be done because $8^{l}$ is a huge
overestimation; for most
$\alpha_{1}...\alpha_{k},\alpha^{\prime}_{1}...\alpha^{\prime}_{k^{\prime}}$
the trace

${\rm
tr}\bigg{(}\prod_{i=1}^{k}h_{\alpha_{i}}\rho(0)\prod_{i=0}^{k^{\prime}-1}h_{\alpha^{\prime}_{k^{\prime}-i}}\,\hat{n}_{1,2}\bigg{)}$
(119)

is zero. Moreover, direct algebraic calculation shows that (119) is nonzero
only if the operator

$\Theta=\prod_{i=1}^{k}h_{\alpha_{i}}\prod_{i=0}^{k^{\prime}-1}h_{\alpha^{\prime}_{k^{\prime}-i}}$
(120)

is Hermitian. For example for $l=2$ from 192 terms we get 18 nonzero terms.