Gaussian Thermal Operations and the Limits of Algorithmic Cooling

A. Serafini, M. Lostaglio, S. Longden, U. Shackerley-Bennett, C.-Y. Hsieh, and G. Adesso Department of Physics & Astronomy, University College London, Gower Street, WC1E 6BT, London, United Kingdom ICFO Institut de Ciències Fotòniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels, Spain School of Mathematical Sciences and Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems, University of Nottingham, University Park Campus, Nottingham NG7 2RD, United Kingdom (Dated: January 22, 2020)

Introduction and Summary -The past few years have witnessed a resurgence of studies into the thermodynamics of quantum systems [1], which have lent novel insight into the nature of thermodynamic relations, as well as into the role of thermodynamic quantities such as temperature, entropy and work [2,3], set against the practical backdrop of realising superior thermal machines operating in the quantum regime [4].
A key ingredient to any attempt to analyse these questions beyond the limited scope of a specific model is the characterisation of a class of "thermal operations", i.e., of operations that can be realised with the aid of the surrounding environment [5,6] (for reviews, see [7,8]). Whilst the frameworks resulting from this approach may yield significant wisdom concerning the ultimate limitations that constrain thermal scenarios, they are at times fraught by a certain 'lack of realism', in that they include interaction Hamiltonians which are not necessarily encountered in practice. Also, they are limited to finite-dimensional settings. It is therefore desirable to single out and characterise subclasses of thermal operations with direct practical relevance.
To this aim, this paper shall consider the subclass of Gaussian thermal operations (GTOs), i.e., the class of operations on continuous variable systems obtained by considering energy-preserving bilinear interaction Hamiltonians between the system and a thermal environment. This subclass is extremely relevant in practice, given that quadratic Hamiltonians, which generate Gaussian unitaries, are very common and that system-bath interactions are often linear or may be linearised, especially in quantum optics and analogous set-ups. Indeed, various experimental platform relevant to quantum thermodynamics operate in the Gaussian regime. Examples include cavity optomechanics [9,10], Bose-Einstein condensates loaded into cavities [11], ions in harmonic traps [12,13]. In view of the same practical reasons, work extraction, storage and fluctuations [14][15][16][17], entropy production [18], heat transport [19], thermometry [20] and fluctuation-dissipation theorems [21] have been investigated in Gaussian scenarios. The capabilities of Gaussian operations in other, not necessarily thermodynamical, settings are also being considered [22,23].
The other defining feature of GTOs, alongside Gaussianity, is energy-preservation. As we will prove, this feature implies that GTOs may be physically reproduced through operations corresponding, in the optical picture, to passive optical elements (i.e., semi-reflectant mirrors and dielectric plates), even in presence of thermal noise (which can be represented as a beam splitter coupling the input with a thermal mode). Alternately, dynamics equivalent to GTOs can be obtained by contact with a Markovian thermal reservoir (giving rise to the so-called quantum optical master equation, see [24,25]), or in coupled resonant cavities or optomechanical systems with negligible counter-rotating terms [26].
In this paper, we shall achieve a compact, constructive characterisation of the most general GTO on any number of modes. We shall see that such a characterisation becomes particularly simple for systems with nondegenerate eigenfrequencies, where it can be cast as a single-mode property. We shall then derive necessary and sufficient conditions for state transformation on singlemode systems and then proceed to analyse the possi-bilities offered by algorithmic cooling in the Gaussian regime, through alternating GTOs and unitaries. We will prove that, at variance with the finite-dimensional case [27], in the absence of ancillas no such strategy can cool the system below the environmental temperature. Sideband-like strategies involving high-frequency ancillas or higher order interactions are necessary to such an aim. Also in view of the ubiquity of Gaussian evolutions as a complete toolbox for quantum technologies [28] and in the modelling of open quantum systems of harmonic systems, the fundamental limitations to cooling techniques we will establish in the Gaussian regime possess a direct practical interest.
Gaussian systems -We will consider bosonic continuous variables encoded into vectors of self-adjoint operatorsr = (x 1 ,p 1 , . . . ,x n ,p n ) T obeying the canonical commutation relations [r,r T ] = iΩ, where the commutators are taken between all pairs of elements ofr (as in an outer product) and form the non-degenerate, antisymmetric symplectic form Ω, with Ω = Ω ⊕n 1 and Ω 1 = 0 1 −1 0 [25]. A second-order HamiltonianĤ is one that may be written as a second-order polynomial ofr: for a symmetric Hamiltonian matrix H and a real vector d. Gaussian states are then defined as the ground and (Gibbs) thermal states of second-order Hamiltonians, and are completely characterised by a vector of first moments r = r and the covariance matrix (CM) σ = {(r − r), (r − r) T } where, again, the anticommutators {·, ·} are taken between each pair of operator entries to form the symmetric, real matrix σ, satisfying σ + iΩ ≥ 0 [25]. Gaussian unitary operations -ones that map Gaussian states into Gaussian states -are those generated by second-order Hamiltonians and admit a symplectic representation: their action on the second moments may be written as σ → SσS T , where S ∈ Sp 2n,R (i.e., S is such that SΩS T = Ω). It is well known that any positive-definite real matrix P may be put into 'normal modes' by congruence with a symplectic transformation: ∃ S : where the ν j 's are the 'symplectic eigenvalues' of P ; if P is a Hamiltonian matrix, the quantities ν j represent the eigenfrequencies of P (the frequencies of its normal modes). In the case of the CM of a quantum state, one has ν j ≥ 1 (an expression of the uncertainty principle). Bear in mind that the spectrum of a Gaussian state is entirely determined by its symplectic eigenvalues and that tensor products at the Hilbert space level translate into direct sums in the Gaussian and phase space descriptions.
In the following, a major role will be played by the set of orthogonal symplectic transformations, for which SS T = 1 (also known as "compact", or "passive" transformations, as they do not require any source of energy in standard optical implementations). Further specific notation will prove convenient: we shall adopt the shorthand notation S[σ] = SσS T and the symbol Tr b to de-note partial tracing of the bath's degrees of freedom in the phase space, which just corresponds to pinching out the relevant part of a CM, discarding the rest.
Let us also recall that the most general deterministic Gaussian CP-map, obtained by letting the system interact with an environment in a Gaussian state through a quadratic interaction Hamiltonian, is characterised, up to arbitrary displacements of the first moments, by the mapping σ → XσX T +Y , with Y ≥ −iXΩX T +iΩ [25]. The first aim of this paper will be characterising the subclass of deterministic Gaussian CP-maps that are also thermal. A particularly relevant class of single-mode channels, which will play a prominent role in what follows, is the so-called 'phase-covariant' ones, where X = x1 2 and Y = y1 2 , with y ≥ |1 − x 2 | [25,29], (throughout the paper, the symbol 1 d denotes the identity matrix in dimension d).
The class of Gaussian thermal operations -Given a second-order system HamiltonianĤ s and an inverse temperature β = 1/(kT ) (where k is the Boltzmann's constant and T is the environment's temperature), we shall define GTOs as the operations obtained by: • Preparing an environmental ancilla with arbitrary second-order HamiltonianĤ b in the Gibbs state e −βĤ b /Tr e −βĤ b .
• Letting system and bath interact through an energy preserving Gaussian unitaryÛ I such that [Û I ,Ĥ s +Ĥ b ] = 0.
The maps above arise naturally through contact with thermal reservoirs where the interactions are well described by polynomials of order two in the canonical operators, whose importance has been already remarked. Note that all energy preserving Gaussian unitaries can be written asÛ I = e iĤ I t for some t ≥ 0 andĤ I a Hamiltonian of order two in the canonical operators satisfying [Ĥ I ,Ĥ s +Ĥ b ] = 0. Note also that the definition above coincides with the well-established definition of thermal operations [5,6], once the restrictions to second-order operations are lifted. ArbitraryĤ s ,Ĥ b andĤ I of order two are parametrised by the symmetric Hamiltonian matrices H s , H b and H I and the vectors d s , d b and d I . We will further restrict the Hamiltonian matrices of system and environment to be strictly positive. Hamiltonian matrices with negative eigenvalues correspond to Hamiltonian operators that are not bounded from below, and thus do not even admit a well-defined Gibbs state, so their exclusion is not a restriction. Positive semi-definite, but not strictly positive, Hamiltonian matrices correspond to a set of measure zero within the Gaussian realm, with Gibbs states that are not regular, trace-class Gaussian states and thus do not give rise to Gaussian CP-maps. It might still be possible to obtain legitimate operations from non-positive system Hamiltonians, but we shall disregard such peculiar cases in this treatment.
First-order terms in the interaction Hamiltonian generate displacements (shifts in the first-moment vector r). Since no first-order term commutes with a strictly positive quadratic Hamiltonian (linear displacements do affect the energy of trapped systems), displacements must be severely limited if they are to give rise to thermal operations. Rather than complicating our treatment with the inclusion of first-order terms, which do not add anything conceptually remarkable, we defer such a discussion to the Supplemental Material (SM) [30], and set all first order terms d s , d b and d I to zero to present our main results. Simulating Gaussian thermalisations -Within the above restrictions, a GTO generally involves an arbitrary number of bath modes, as well as an arbitrary sequence of second order interactions between these and the system modes. A crucial question is then if there exists a simpler protocol able to reproduce every Gaussian thermalisation with less extensive resources. Our first main result answers this question in the affirmative, presenting a very simple scheme able to simulate exactly a general GTO (recall the shorthand notation whereby symplectics act by congruence): T H sr be a system Hamiltonian with normal form l ω l 1 2n l = S −1 H s S T−1 , where n l ∈ N is the mode degeneracy of the eigenfrequency ω l and S ∈ Sp 2n,R for n = l n l . The class of GTOs at background inverse temperature β act on the system CM σ as where the direct sum runs over the distinct eigenfrequencies and, setting ν l = e βω l +1 e βω l −1 : 1. Each Φ l are phase-covariant CP maps [29], acting on the l−th eigenfrequency space as Φ l (σ) = . W l and Z l are passive symplectic transformations acting on the system's set of modes associated with the l-th eigenfrequency.
Let us now unravel this statement and the restrictions it poses on the structure of GTOs, which will also allow us to sketch the main lines of its proof (whose full details are found in the SM [30]). The transformation S is just the one bringing the system Hamiltonian into normal modes, set by the given system quadratic Hamiltonian [31].
The first step towards the statement above is realising that, once both system and ancillas are cast into normal modes, all GTOs are obtained by letting the n l system modes pertaining to the same eigenfrequency ω l interact with an equal number n l of environmental normal modes at the same frequency: being a passive symplectic transformation acting on the system plus bath degenerate eigenfrequency subspace labelled by l (of dimension 2n l ). Very significantly, normal modes belonging to different eigenfrequency sectors do not interact during thermal operations (this holds regardless of any correlations that may exist between the physical bath modes).
The second step to obtain the compact characterisation above is that, due to the symmetries of the problem at hand, each O l admits a very simple structure: where, as already stated, W l and Z l are passive symplectic on the system, and M l is a set of beam splitters independently mixing each mode j = 1, ..., n l with a corresponding mode of the environment: kk denotes a beam splitter mixing system mode k (with ladder operatorâ k = (x k + ip k )/ √ 2) with bath mode k (with ladder operatorb k ); at the Hilbert space level,R Thus, in a GTO, each oscillator within the degenerate frequency sector is mixed with a correspondent thermal oscillator by means of a beam splitting operation. Tracing out the bath after such an interaction gives rise to the tensor product of phase-covariant channels that were denoted with Φ l . What is perhaps surprising is that every GTO can be simulated in this simple way, by independent interactions with the environmental modes. Besides, since the loss channels Φ l are Markovian [25], the most general Gaussian thermalisation can be generated by a simple Markovian master equation. Indeed, GTOs are the most common as well as easiest to implement transformations, corresponding, in the normal mode basis, to passive optics or loss to a thermal Markovian reservoir. Single-mode criteria -For each non-degenerate system eigenfrequency, a GTO reduces to a single-mode transformation. All single-mode passive transformations are phase shifters, and the transformation Z l may always be simplified by left-multiplication with phase shifters (see [30]) and may thus, on a single-mode, be reduced to the identity without loss of generality. Hence, the most general GTO on a non-degenerate eigenfrequency subspace takes a very simple form indeed: T SS Tr be a single-mode system Hamiltonian, then the class of GTOs is given by e βω −1 and D ϕ = cos ϕ sin ϕ − sin ϕ cos ϕ .
We can now spell out the full criterion for Gaussian state transformations through single-mode GTOs. That is, given an input CM σ i and an output CM σ f , is there a GTO mapping σ i into σ f ? Here, one should recall that the most general single-mode CM σ may be written as a rotated and squeezed thermal state: with ν b = (e βω + 1)/(e βω − 1). Note that the parameters ϕ i and ϕ f are irrelevant to the transformation criterion, which admits a simple geometrical representation: if one parametrises the class of single-mode Gaussian states (in the basis of normal modes ofĤ s and modulo phase shifters) in the twodimensional space (νz, ν/z), one can thermally map the states (ν i z i , ν i /z i ) only into states lying along the seg- [30]). Note that the squeezed states to which this criterion applies display quantum coherence (off-diagonal elements) in the energy eigenbasis. In the case with no squeezing, where the states have no coherence in the energy eigenbasis, the transformation criterion reduces to ν f ∈ [ν b , ν i ]. In physical terms, this is equivalent to stating that GTOs send an initial thermal state at temperature T i into a final thermal state at temperature T f falling between T i and the environment's temperature T . This complies with the thermo-majorisation and the many second laws criteria of [32,33] (see [30]), while the case with squeezing falls beyond the criteria's scope. Interestingly, the prediction that T f must fall between T i and T differs from what happens in qubit systems and turns out to be crucial for the task of cooling, to which we now turn. Algorithmic cooling -Let us now discuss the main repercussions of the characterisation derived above on the algorithmic cooling of Gaussian systems. In the spirit of heat-bath algorithmic cooling (HBAC) [34] one aims at cooling a system by alternating Gaussian unitaries and thermal operations which, if one allows for partial rather than complete thermalisations, may lead to improvements in the cooling of finite-dimensional systems [27,35]. For example, a single qubit can be cooled arbitrarily close to the ground state by applying to it Pauli x unitaries interspersed with thermal operations, without the need of extra ancillas. In fact, at low enough temperatures, the required thermal operations can be approximated by resonant Jaynes-Cummings couplings to a single, initially thermal oscillator. A natural question is then if a single system oscillator can be cooled below the environment temperature in a similar fashion; that is, by unitaries on the system acting between the GTOs. This would be particularly advantageous because it would only require standard quadratic interaction Hamiltonians.
Here we answer this question in the negative for singlemode systems: If U j are single-mode (not necessarily Gaussian) unitaries and T j arbitrary single-mode GTOs, cannot be cooled below the minimum between the environment's entropy and the initial system entropy. This is the case since the output entropy of phase-covariant, single-mode Gaussian channels at given input entropy is minimised by (Gaussian) thermal inputs (with respect to the normal mode Hamiltonian) [36], with optimal output entropy that is monotonic in the input entropy. Thus, the best the unitaries U j can do is put the state in normal form which, for given initial symplectic eigenvalue ν j , yields the output symplectic eigenvalue pν j + (1 − p)ν b ≥ min{ν j , ν b }, so that the minimum entropy is obtained by either shielding completely from the environment or by complete thermalisation. Notice that, rather remarkably, such an entropic bound holds for any unitary operation and any input state, not necessarily Gaussian. We also show in the SM [30] that the impossibility of lowering the system entropy below the environment's value is maintained if one extends the class of thermal operations to include single-mode squeezed baths, which are not encountered spontaneously in nature but may be engineered under certain controlled conditions [37][38][39][40][41].
Cooling opportunities open up if non quadratic interaction Hamiltonians or control over the energy levels' structure are allowed, as is commonly assumed for quantum refrigerators [42], or if some of the thermal ancillary modes can be manipulated by general Gaussian unitaries. In point of fact, these latter schemes, unless restricted by practical constraints, allow one to always cool any oscillator arbitrarily close to the ground state. To this aim, one may in principle include a thermal ancillary mode at high enough frequency so that its entropy is arbitrarily low, and then swap such a low entropy state into the system through a beam splitter acting in the unitary step (notice that such an interaction between modes at different frequencies would not be prohibited, since the unitary does not have to be a thermal operation in the general set-up we are considering). This is nothing but the discrete version of sideband cooling, where excitations are extracted from the system of interest (such as a mechanical oscillator) into a coupled oscillator (such as a mode of light, in optomechanical set-ups) at higher frequency, from where they leak to the environment.
Our no-go theorem complements the impossibility of engineering absorption refrigerators with Gaussian resources alone, pointed out in [43]. Our treatment is broader, relying on the general GTOs rather than on a specific time-evolution, and focuses on the system temperature, rather than heat transport between reservoirs. Conclusions and outlook -We presented a full characterisation of Gaussian thermal operations, implying that they are all generated by a simple, time-local master equation, determined necessary and sufficient conditions for transformation under GTO on a single-mode and proved that no algorithmic cooling acting on a singlemode system alone can ever lower the entropy below the background or initial ones, a fact which is relevant in practice given the broad applicability of noise models based on bilinear interactions with an environment. The latter finding is intimately related to the fact that GTOs are all Markovian. As such, any dynamical trajectory reaching the thermal state must terminate there. In fact, the cooling protocol for a single qubit presented in [27] relied precisely on the fact that system-bath correlations can be used to cross the thermal state and achieve temperatures lower than that of the environment. This possibility is precluded, for Gaussian systems, by our nogo result. Our framework, however, sets up the scene to explore transformation conditions and more articulate cooling schemes in multimode scenarios (we refer the reader to the final section of the SM [30] for a detailed discussion of future perspectives). Note added -During the completion of this article, we became aware of closely related work [44], where thermal transformations are constrained to passive unitaries by design and several multimode necessary conditions for state transformation are discussed.

FIRST-ORDER TERMS
Since the main text considers only the purely quadratic case, let us discuss here the effect on thermal operations of terms of the first-order in the canonical operators. First-order terms in the bath Hamiltonian are immaterial, since they can always be set to zero by a local unitary operation (a local phase-space displacement). They can therefore be disregarded without loss of generality, as it has been done in the paper.
Any system Hamiltonian with first-order terms, such asĤ s = 1 2 (r − r) T H s (r − r) can be written asĤ s = D † r 1 2r T H srDr , for the unitary displacement operatorD † r = e ir T Ωr , which indeed just displaces the canonical operators by real quantities. Thermal operations with respect to such a displaced Hamiltonian are therefore just given by where T is the thermal operation with respect to the corresponding centred Hamiltonian (with no first-order term), as derived in the main text. Clearly, the displacement does not generally commute with T , so that the net effect of a thermal operation will involve a finite displacement of the first moments (which would be very easy to evaluate in specific cases).
In the main text, we also stated without proof that no Hamiltonian with strictly positive Hamiltonian matrix commute with displacement operators: this is immediately apparent since any translation ofr inr T Hr will always produce a nonzero shift to the value of the operator if H > 0. This would not be the case for a semi-definite Hamiltonian, such as the free Hamiltonianp 2 , which is obviously invariant under translations of thex operator.

SYMPLECTIC RENDITION OF THERMAL OPERATIONS
We can work in the local system and bath symplectic bases where the local Hamiltonian matrices are in normal form, and then consider the most general interaction Hamiltonian matrix H I . For the bath, this can be done without loss of generality, since it just corresponds to a choice of basis of a subsystem which will be ultimately traced out. For the system, such an assumption will be relaxed by including the action of the symplectic S that brings the system Hamiltonian to normal modes.
In such bases, one has H s = l ω s,l 1 2n l and H b = l ω b,l 1 2m l . Bear in mind that, because we allow for the addition of ancillary modes with arbitrary Hamiltonians, the bath eigenfrequencies ω b,l and degeneracies m l are whatever we like them to be. In other words, the only input parameters determining the set of thermal operations are the system's eigenfrequencies ω s,l and degeneracies n l , as well as the inverse temperature β.
Notice now that a necessary condition for the Hamiltonian operatorĤ I to commute withĤ s +Ĥ b is that the unitary transformations generated by exponentiating iĤ I leaveĤ s +Ĥ b unchanged. In terms of quadratic Hamiltonians, this is equivalent to stating that the symplectic transformations e −ΩH I t [25] must belong to the subgroup of transformations that leave H s ⊕ H b = l ω s,l 1 2n l l ω b,l 1 2m l unchanged when acting by congruence. But such an isotropy group is easily characterised: Lemma 1 -Isotropy group of normal form matrices. The symplectic isotropy group of the transformation Y = l ω l 1 2d l is given by the direct sum of the compact symplectic subgroups K(2d l ) = Sp 2d l ,R ∩ SO(2d l ), each acting on the 2d l -dimensional subspace pertaining to a certain eigenfrequency ω l . Proof. Let K be a symplectic transformation part of the isotropy group. Then, by definition KY K T = Y and KΩK T = Ω. Recalling that K is invertible, it is easy to show that the previous two equations imply [K, Y Ω] = 0. If K is written in terms of 2 × 2 sub-blocks K jk , as per then the simple form of Y Ω allows one to reduce the commutation condition with K to a condition on the sub-blocks: (where Ω 1 is the 2 × 2 symplectic form on a single mode). Writing this yields the set of equations which, for ω j = ω k , imply K jk = 0. Therefore, the isotropy transformation K must be block-diagonal with respect to subspaces associated with distinct symplectic eigenvalues of Y , and must be a direct sum of symplectic orthogonal transformations on each such subspace (since any such transformation clearly preserves Y ).
Let us remark that one may show that all of these isotropy transformations are generated by Hamiltonians that commute with the Hamiltonian they preserve, so that each of them does indeed define a legitimate Gaussian thermal operation. The orthogonal symplectic transformations that form the isotropy subgroups are also referred to as "passive" in the quantum optics tradition, since they preserve the number of photons.
By virtue of the statement above, Gaussian thermal operations act separately on each of the system's phase space sectors pertaining to a different eigenfrequency. Besides the passive, symplectic transformations acting on such subspaces, which are obviously all thermal, less trivial examples of Gaussian thermal operations are obtained by appending to each degenerate subspace with eigenfrequency ω l a set of bath modes at the very same frequency ω l . Such modes are all prepared, before the unitary transformation, in the (Gaussian) thermal Gibbs state with covariance matrix ν l 1 2 , with ν l = (e βω l + 1)/(e βω l − 1), and we can add as many as we like (see Ref. [25] for the formula relating frequency and temperature to the symplectic eigenvalue).
In order to complete our characterisation of Gaussian thermal operations, we now set out to characterise the set of Gaussian CP-maps obtained by letting an input Gaussian state of n modes, with arbitrary covariance matrix σ, interact with an environment with covariance matrix ν1 2m , through a global passive symplectic transformation, for all integer m.

UNITARY REPRESENTATION OF THE COMPACT SUBGROUP
It is well known that, by adopting a representation in terms of annihilation and creation operators, passive symplectic transformations in dimension 2d may be represented as U 0 0 U * , where U ∈ U (d) (in the field theory tradition, this is known as the 'Bogoliubov' representation of passive symplectic operations) [25]. Such an isomorphism between K(2d) and U (d) will be very advantageous in describing arbitrary passive symplectic acting on the degenerate eigenfrequency sectors of the system plus bath Hamiltonian.
In this notation, which, in each eigenfrequency sector, corresponds to taking the basis of operators (â s 1 , . . . ,â s n l ,â b 1 , . . . ,â b m l ,â †s 1 , . . . ,â †s n l ,â †b 1 , . . . ,â †b m l ) T (with s denoting the system and b the bath), the global, initial CM describing system and bath may be written as where σ a l a † l is an n l × n l hermitian matrix reporting the values of the symmetrised covariances of all pairs of system annihilation and creation operators (one each), whilst σ aa contains the covariances of pairs of annihilation operators.
The blocks ν l 1 m l correspond to the covariances of the initial thermal state of the m l bath modes.
Inspection of the initial CM above reveals that the CP-map obtained by letting such an initial state evolve through a global passive represented by U is invariant under right multiplication of U by an arbitrary bath unitary Z m l . Besides, one can also left-multiply U by another, generally different, bath unitary W m l , since the bath is ultimately traced out (corresponding, in the CM formalism, to pinching the relevant part of the matrix). These symmetries will be key to what follows.
Note also that, under such a choice of basis, 'standard' beam splitters may be represented as real two-dimensional rotations, which we shall denote with the letter R below: R = cos θ sin θ − sin θ cos θ for θ ∈ [0, 2π[ (acting on the relevant components, which will be specified through indexes below).

SIMPLIFYING THE UNITARY MATRIX
As we just saw, a global (n l + m l )-dimensional unitary U , that determines the thermal Gaussian CP-map by acting on system and bath in a certain eigenfrequency sector, may be simplified by acting on the left and right through a local, bath unitary, as in the lemma below.
Because of the previous lemma we can restrict, without loss of generality, to baths with the same number of modes as the system (n = m) [49]. Now, a U with off-diagonal blocks of the form above can be decomposed as follows.
Lemma 3 -Cosine-sine decomposition. Let U be a 2n × 2n unitary matrix, then where W , X, Z and Y are n×n unitary matrices, while R jj is a (real) beam splitter between the j-th and the (n+j)-th mode.
This is a standard decomposition of unitary matrices, which follows from taking the singular value decomposition of the two n × n off-diagonal blocks through the local unitaries and then apply the unitarity conditions (see, e.g., [45]).
Let us notice, en passant, that minor variations of the lemma above may be employed to obtain an explicit proof of the well known results that (i) any unitary may be decomposed into two-level unitaries, and (ii) any passive symplectic transformation is the product of beam splitters and phase shifters.
Lemma 3 is incredibly revealing to the purpose of simplifying Gaussian thermal operations: indeed it is telling us that, in each eigenfrequency subsector, and up to an initial and final passive symplectic acting on the system (Z and W , respectively), the action of a thermal map boils down to mixing each system normal mode with a bath mode, independently, through a standard beam splitter. The local unitary transformations on the bath X and Y can be completely disregarded: the former because it acts at the very end, just before the bath is traced out, the latter because the initial bath state, given by a thermal state on modes with degenerate normal frequency and hence with CM proportional to the identity, is invariant under passive transformations.
Note also that the decomposition above is slightly redundant, as it involves 4n 2 + n real degrees of freedom (n 2 per unitary, plus n for the n mixing angles of the beam splitters), against the 4n 2 degrees of freedom of a 2n-dimensional unitary. In fact, one of the four unitaries is not completely arbitrary, but can be simplified by multiplication on a side by any diagonal matrix of complex phases (corresponding to a tensor product of single-mode phase shifters in physical set-ups). It is easy to see that such a multiplication may be absorbed by redefining the other unitaries without affecting the singular values of the off-diagonal blocks (which, effectively, set the beam splitters's angles). To our purpose, it will be convenient to simplify the matrix Z, although W might also have been chosen. FIG. 1. Schematics of a GTO acting in the normal-mode basis: (a) a 5-mode system, with degenerate eigenfrequencies ω1, pertaining to two modes, and ω2, pertaining to three modes, undergoes the initial passive symplectics Z1 and Z2, followed by a tensor product of phase-covariant channels Φ1 and Φ2; in turn, each Φj is the tensor product of phase-covariant channels θ jk , each acting on a mode separately; finally, the passive symplectics W1 and W2 act separately on the degenerate eigenspaces; (b) each phase-covariant channel θ jk is shown to result from the mixing of the system mode at a beam splitter, whose transmittivity sets the parameter θ jk (here, for simplicity, the parameter θ jk also denotes the single-mode channel itself).

PARAMETRISATION OF GENERAL GAUSSIAN THERMAL OPERATIONS
All the above was derived for the system normal modes, whose local Hamiltonian matrix we shall denote hereafter with ω s . The most general local Hamiltonian matrix is therefore Sω s S T , where S is any local symplectic transformation on the system. Above, we determined and simplified all of the global symplectic transformations S I that preserve ω S ⊕ H b , where H b is the bath Hamiltonian matrix. It follows that the whole set of global symplectic that preserves a general quadratic Hamiltonian not in normal form, as given above, is just This fact, along with the decomposition (12) and the basic piece of knowledge that a beam splitting interaction with an environmental mode gives rise to the phase-covariant CP-map Φ, that maps a single-mode CM σ according to Φ(σ) = cos 2 θσ + sin 2 θν b 1 2 , leads directly to the general characterisation of Gaussian thermal operations given by Theorem 1, which is illustrated in Fig. 1.
As explained in the previous section, there is some residual freedom in the constructive characterisation of Theorem 1. Because of the residual ambiguity in the cosine-sine decomposition, whilst the operations W l may be taken as completely arbitrary passive symplectic transformations, the transformations Z l are passive symplectic operations that can be simplified by the action of a tensor product of phase shifters acting on them from the left: each of them thus bear n 2 l − n l free parameters (recalling that n 2 l is the number of parameters in an arbitrary passive symplectic transformation). Therefore, up to the transformation S, a GTO acting on a degenerate eigenfrequency sector comprising n l modes is parametrised by 2n 2 l + 1 parameters (one of them being the inverse temperature β).
The single-mode case For n l = 1, which covers all systems with non-degenerate eigenfrequencies, the only local passive transformation is the phase shifter D ϕ given, in the (x,p) basis, by D ϕ = cos ϕ sin ϕ − sin ϕ cos ϕ .
As discussed above, the passive transformation Z l entering Eq. (2) of the main text can be simplified through leftmultiplication by a phase shifter, and may thus be reduced to the identity without loss of generality in the single-mode case. Setting p = cos 2 θ, one is therefore left with the expression reported in Proposition 1 and Eq. (3) of the main text.

SINGLE-MODE STATE TRANSFORMATIONS
Let us restate the most general thermal mapping for a non-degenerate (single-mode) system frequency: which has been written in terms of the initial and final CMs σ i and σ f in view of our next objective, which is characterising allowed mappings between pairs of states at given ν b (temperature). Clearly, one can re-write the initial and final CMs in the normal basis to obtain a condition independent from S. Formally, one can act on the left and right hand sides with S and obtain a condition for the transformed input and output σ i,f = S −1 σ i,f S −1T : Single-mode Gaussian states are particularly simple, as can be seen by applying the symplectic singular value decomposition to the Williamson form of a state [25]. Their most general form is where ν i,f are the initial and final symplectic eigenvalues (which determines any entropy in the single-mode case), D i,f are single-mode rotations and Z i,f = diag(z i,f , z −1 i,f ), and we can assume z i,f ≥ 1 without loss of generality (since phase space rotations allow one to invert z i,f ).
Since thermal mappings are rotationally invariant in phase space, one can always match the optical phases of input and output, and we can therefore disregard the rotations altogether. One is then left with the following necessary and sufficient conditions for state transformations:

Isotropic states
In the absence of squeezing (z i,f = 1), the situation is very simple to depict, as the conditions above lead to the necessary and sufficient condition that ν f must lie between ν b and ν i .
Note that, for single-mode Gaussian states, the free energy F in the normal mode basis (at eigenfrequency ω) may be easily expressed as (see [25] for a formula expressing the von Neumann entropy of a Gaussian state as a function of the symplectic eigenvalue ν b ) For z = 1 and at given β, such a function of ν b has a single minimum at the environmental value ν b = e βω+1 e βω −1 . Therefore, the transformation criterion ν f ∈ [ν b , ν i ] (regardless of the ordering of ν b and ν i ) tells us that, even in the absence of squeezing, the decrease in the free energy is necessary (as it always is, since thermal operations have thermal fixed points) but not sufficient for two Gaussian states to be thermally connectable through an environment at inverse temperature β. Transformation criterion for single-mode systems. The shaded area contains all single-mode Gaussian states which, up to rotations and first moments, are parametrised by the symplectic eigenvalue ν ≥ 1 and squeezing parameter z ≥ 1 (in the normal-mode basis of the system Hamiltonian). Non-squeezed states (which, for zero first-moments, are diagonal in the energy eigenbasis) lie on the z = 1 line that bisects the two axes. Given an environmental symplectic eigenvalue ν b (set by frequency and temperature), one has that an initial state parameterised by νi and zi may be mapped into a final state with ν f and z f if and only if the point (ν f z f , ν f /z f ) belongs to the segment connecting (νizi, νi/zi) to (ν b , ν b ) (represented in red between diamonds on the graph).
Notice also that, since such states are diagonal in the energy eigenbasis, the hierarchy of free energy criteria pointed out in [33] will apply to them. However, under the additional assumptions of a single-mode system in a Gaussian state, all such thermal transformation criteria coalesce to a single one, since all Renyi entropies are determined by a single quantity.
General squeezed states Solving the system above for p yields whose boundedness (0 ≤ p ≤ 1) gives the necessary and sufficient conditions for possible transformations. Direct inspection of (17) reveals the whole geometric nature of such a necessary and sufficient condition, illustrated in Fig. 2. Given ν b , as well as the input ν i and z i , it is convenient to parametrise the possible output state in the space ν f z f and ν f /z f , for z f ≥ 1. For a thermal mapping to be possible, it is necessary that such variables belong to the interval [ν b , ν i z i ] and [ν b , ν i /z i ] (denoting, up to the proper ordering, the interval between the two values). The necessary and sufficient condition is that (ν f z f , ν f /z f ) belong to the diagonal of such an interval, joining (ν b , ν b ) to (ν i z i , ν i /z i ). Notice that for z i = 1 the interval becomes a square and the conditions reduce to z f = 1 and ν f ∈ [ν b , ν i ]. The effect of the initial squeezing is precisely to make such a square oblong.
Simple necessary conditions about ν f and z f may also be obtained as follows. Taking the ratio of the two equations in (17), one gets which can be shown by observing that the derivative of the function above with respect to p is the always positive [(1−p)ziν b +pνi] 2 (recalling that z i ≥ 1). Clearly, mixing with a non-squeezed state cannot increase the squeezing.