Locally optimal symplectic control of multimode Gaussian states

TL;DR

This paper develops locally optimal quantum control strategies using Gaussian unitaries to accelerate the relaxation of multimode Gaussian states in lossy channels, highlighting the role of squeezing in optimizing entropy and temperature dynamics.

Contribution

It derives conditions for optimal control of Gaussian states in lossy channels, providing explicit relaxation times and emphasizing the importance of squeezing in the process.

Findings

Squeezing enhances the rate of entropy change during relaxation.

Exact relaxation times are computed for optimal heating and cooling.

Control strategies are tailored to the parameters governing entropy and temperature.

Abstract

The relaxation of a system to a steady state is a central point of interest in many attempts to advance control over the quantum world. In this paper, we consider control through instantaneous Gaussian unitary operations on the ubiquitous lossy channel, and find locally optimal conditions for the cooling and heating of a multimode Gaussian state subject to losses and possibly thermal noise. This is done by isolating the parameters that encode entropy and temperature and by deriving an equation for their evolution. This equation is in such a form that it grants clear insight into how relaxation may be helped by instantaneous quantum control. It is thus shown that squeezing is a crucial element in optimising the rate of change of entropic properties under these channels. Exact relaxation times for heating and cooling are derived, up to an arbitrarily small distance from the fixed point of the lossy channel with locally optimal strategies.

Equations

$$\Omega =\underset{i=1}{\overset{n}{⨁}}\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right).$$

$${\hat{\rho }}_G=\frac{e^{-\beta \hat{H}}}{\mathrm{Tr}[e^{-\beta \hat{H}}]},$$

$$\mathbf{d}=\mathrm{Tr}[\hat{\mathbf{r}}{\hat{\rho }}_G],$$

$$\sigma =\mathrm{Tr}[\{(\hat{\mathbf{r}}-\mathbf{d}),(\hat{\mathbf{r}}-\mathbf{d}{)}^{⊺}\}{\hat{\rho }}_G].$$

$$S\Omega S^{⊺}=\Omega .$$

$$\sigma \to S\sigma S^{⊺}.$$

$$H_C=\left(\begin{array}{cc}0& C\\ C^{⊺}& 0\end{array}\right).$$

$$\dot{\sigma }=A\sigma +\sigma A^{⊺}+D,$$

$$A=\frac{1}{2}\Omega C\Omega C^{⊺}\text{and}D=\Omega C{\sigma }_B{C}^{⊺}{\Omega }^{⊺},$$

$${\hat{H}}_C=\sqrt{\eta }(\hat{a}{\hat{b}}^{†}+{\hat{a}}^{†}\hat{b}),$$

$$\dot{\sigma }=-\sigma +\chi \mathbb{I},$$

$$\sigma =\chi \mathbb{I}+(\sigma (0)-\chi \mathbb{I})e^{-t}.$$

$$\sigma =SW{S}^{⊺},$$

$$W=\underset{i=1}{\overset{n}{⨁}}{\nu }_i{\mathbb{I}}_2,$$

$${\vartheta }_k[X]:=\sum _{E\in {\mathcal{E}}_k^m}\prod _{j\in E}{\lambda }_j,$$

$$E\in {\mathcal{E}}_k^m\subset P({\mathbb{N}}_m)\text{iff}|E|=k,$$

$${\vartheta }_3[X]={\lambda }_1{\lambda }_2{\lambda }_3+{\lambda }_1{\lambda }_3{\lambda }_4+{\lambda }_1{\lambda }_2{\lambda }_4+{\lambda }_2{\lambda }_3{\lambda }_4.$$

$${\vartheta }_{2k}[\Omega \sigma ]=\sum _{E\in {\mathcal{E}}_k^n}\prod _{j\in E}{\nu }_j^2.$$

$$c_k[X]$$

$${\vartheta }_k[X]$$

$$c_k[X]=-\frac{1}{k}\sum _{i=0}^{k-1}\mathrm{Tr}[X^{k-i}]c_i[X].$$

$${\bar{\vartheta }}_{2k}^i[\Omega \sigma ]:=\sum _{E\in {\mathcal{E}}_k^{n,i}}\prod _{j\in E}{\nu }_j^2,$$

$$E\in {\mathcal{E}}_k^{n,i}\subset P({\mathbb{N}}_n\setminus \{i\})\text{iff}|E|=k.$$

$${\bar{\vartheta }}_4^2[\Omega \sigma ]={\nu }_1^2{\nu }_3^2+{\nu }_1^2{\nu }_4^2+{\nu }_3^2{\nu }_4^2.$$

$${\dot{\vartheta }}_{2k}[\Omega \sigma ]=-2k{\vartheta }_{2k}[\Omega \sigma ]+\chi \mathrm{Tr}[S{V}_k{S}^{⊺}],$$

$$V_k=\underset{i=1}{\overset{n}{⨁}}\left({\nu }_i{\bar{\vartheta }}_{2(k-1)}^i[\Omega \sigma ]\right){\mathbb{I}}_2.$$

$$\underset{S\in \overline{\mathrm{Sp}}(2n,\mathbb{R})}{\sup }\mathrm{Tr}[SY{S}^{⊺}]=\sum _{i=1}^{n}2{\zeta }_{{\bar{z}}_i}^{+}{y}_i,$$

$$\mathrm{Tr}[SY{S}^{⊺}]$$

$$=\mathrm{Tr}[Z^2{R}_{2}Y{R}_2^{⊺}].$$

$$P_{kl}=\{\begin{array}{l}1,k\le n,l=2k-1,\\ 1,k>n,l=2(k-n)\end{array}$$

$$Q=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}{\mathbb{I}}_n& i{\mathbb{I}}_n\\ {\mathbb{I}}_n& -i{\mathbb{I}}_n\end{array}\right).$$

Full text

1 Introduction

As quantum technologies advance, it has become increasingly important to understand the dynamics of a system interacting with its environment. It is hoped that increased insight into such setups will provide a more realistic notion of how to control and exploit quantum behaviours. As well as considering the impact of the environment, one would also like to allow for an external controller who may influence the evolution of the state of the system. This beckons the introduction of open-loop quantum control theory, which is a mathematical framework exploring how time-varying Hamiltonians drive dynamics in the absence of feedback. Exploring these together in a particular scenario is the topic of this study.

Quantum control theory of closed systems is a highly active area of research [1, 2, 3, 4, 5] with its roots in the older field of mathematical control theory [6, 7, 8, 9]. Its introduction to open systems has taken place in different contexts, including multilevel discrete systems [10], systems with closed ‘feedback’ control [11], dissipating qubits [12] and in a drive to understand ‘quantum speed limits’ [13]. The current paper explores open system control theory for continuous variable quantum mechanics, and specifically for the class of Gaussian states. These are ubiquitous in physics and serve as a good model for electromagnetic radiation [14, 15, 16], optomechanical systems [17], trapped ions [18] and mesoscopic massive systems seeking the ‘gravitational quantum regime’ [19, 20]. They are fundamental in the study of continuous variable quantum information due to the ease with which they can be described, and also the natural way in which noise may be introduced into the dynamical equations. The evolution of Gaussian states in dissipative quantum channels has been explored before and it is this literature that we build upon [21, 22].

If one models the environment as an infinite thermal bath coupled to the system through a beam-splitter Hamiltonian, then the channel is known as lossy [23]. Lossy evolutions can be described as either ‘heating’ or ‘cooling’, depending on whether the environment has a higher or lower entropy than the initial state of the system. Here we consider the case of $n$ non-interacting bosonic modes, each coupled with a bosonic thermal environment with the same temperature, undergoing such a lossy evolution. We then ask for the locally optimal, open-loop control strategy to minimise the relaxation time to the steady state. In other words, we wish to maximise the rate of change of the entropy toward the fixed point of the dynamics. Control is considered in the idealised regime where only instantaneous transformations which act impulsively on the system are allowed, in particular we restrict the analysis to (possibly non-local) unitaries that are Gaussian preserving. The restriction to finite, instantaneous controls (corresponding to control Hamiltonians proportional to delta functions in time) is not a completely wild abstraction since experimental set-ups certainly exist where a unitary manipulation will typically take nanoseconds, or tens of nanoseconds, while the decoherence rates are in the order of $10-10^{3}$ ${\rm kHz}$, so that the former may fairly be regarded as instantaneous with respect to the time-scales on which the noise acts. In what follows, we also disregard the possibility of any feedback on the quantum system. Let us anticipate that, in the absence of specific constraints, the non-compact nature of the symplectic group allows for an arbitrarily high amount of squeezing to occur in these control operations. To make the analysis more realistic and reduce infinities, this will be capped.

A central result of this study is showing that squeezing is the main parameter of interest in the analysis of heating and cooling dynamics. A highly squeezed state is locally optimal for minimising the relaxation time in a heating channel. Interestingly, squeezing can also be used to drive the state away from the fixed point of a cooling channel. Conversely, undoing the squeezing is shown to be locally optimal for cooling strategies, which is a point of experimental interest. The minimum time it takes to reach the fixed point of the channel within an arbitrarily small distance is derived in both cases. The analysis rests on the derivation of a compact equation for the evolution of entropic quantities – Eq. (24) – whose variations might also be applicable to the evolution of entanglement of Gaussian states under lossy channels, as discussed in the conclusion.

2 Gaussian states

We begin with an outline of Gaussian states, referring to the introductions given in Refs. [14, 24].

Let $\hat{\mathbf{r}}=(\hat{x}_{1},\hat{p}_{1},\ldots,\hat{x}_{n},\hat{p}_{n})^{\intercal}$ be a vector of canonical operators such that $[\hat{x}_{j},\hat{p}_{j}]=i\delta_{jk}$, where $\delta_{jk}$ is the Kronecker delta function. From this we find that $[\hat{\mathbf{r}},\hat{\mathbf{r}}^{\intercal}]=i\Omega$ where

[TABLE]

Quadratic Hamiltonians are defined as those that can be written as $\hat{H}=\frac{1}{2}\hat{\mathbf{r}}^{\intercal}H\hat{\mathbf{r}}+\hat{\mathbf{r}}^{\intercal}\mathbf{a}$ where $H$ is a $2n\times 2n$, real, symmetric matrix and $\mathbf{a}$ is a vector of real numbers. The set of Gaussian states can be defined as the ground and thermal states of positive definite quadratic Hamiltonians:

[TABLE]

where $\beta$ is the inverse temperature of the state (we set Boltzmann’s constant $k_{B}=1$). Such states are referred to as Gaussian due to their Wigner representation which takes a Gaussian form. The Wigner representation immediately suggests that such states should be totally specified by their first and second moments, defined respectively as

[TABLE]

[TABLE]

$\sigma$ is the $2n\times 2n$ covariance matrix of the Gaussian state and $\{\cdot,\cdot\}$ denotes the anticommutator.

Evolution under quadratic Hamiltonians transforms Gaussian states into Gaussian states. The second moment encodes all entropic and entanglement properties and so we restrict our attention, here and in the rest of the paper, to ignore the first moments. In this case it is possible to consider $\mathbf{d}=0$ in both the state and the quadratic Hamiltonian transformation. Once this is done the transformations generated by the remaining Hamiltonians, and their combinations, form the symplectic group $\operatorname{Sp}(2n,\mathbb{R})$. This is the set of $2n\times 2n$ matrices $S$ such that

[TABLE]

Covariance matrices transform under a finite dimensional representation of the symplectic group by conjugation,

[TABLE]

3 Open Diffusive dynamics

The details of Markovian evolution for Gaussian states are presented in Refs. [14, 22] where, following a standard approach in quantum optics, the weak coupling to an infinite bath is replaced with the coupling to a finite number of continuously refreshed environmental modes. This coupling is provided by the Hamiltonian,

[TABLE]

The evolution of the Gaussian state has a representation in terms of the associated covariance matrix which takes the form

[TABLE]

where

[TABLE]

with $\sigma_{B}$ denoting the bath covariance matrix.

Lossy channels occur by setting the coupling Hamiltonian to induce a swapping of field excitations between the system and the bath. We focus on the case of $n$ non-interacting bosonic modes, each coupled with a bosonic thermal environment via an individual exchange Hamiltonian of the form:

[TABLE]

where $\hat{a}$ denotes the annihilation operator of the system and $\hat{b}$ the annihilation operator of the bath mode. This is equivalent to setting $C=\sqrt{\eta}\Omega$ which in turn sets $A=-\eta\mathbb{I}/2$ and $D=\eta\sigma_{B}$. We set the state of the bath as $\sigma_{B}=\chi\mathbb{I}$ where $\chi=2\bar{N}+1$ and $\bar{N}$ is the average photon number in each bath mode. We assume that the $n$ thermal baths have the same $\bar{N}$, i.e. the same temperature. Rescaling the time parameter by $\eta$, the relaxation rate of the lossy channel, we arrive at

[TABLE]

which represents the free, lossy evolution of the Gaussian state with solution

[TABLE]

Eq. (11) will be used extensively in the rest of the paper as we explore the evolution of some important parameters of $\sigma$.

4 Symplectic invariants

An $n$-mode covariance matrix has $n(2n+1)$ free parameters that compose its information content. $n$ of these parameters have the special status of being invariant under symplectic transformations and so play a unique role in describing the state. These symplectic invariants are important in that they encode the entropy of the state, relating to its temperature and mode frequency [25, 26]. We will explore two different ways of representing these invariants allowing the derivation of a new equation from Eq. (11) that isolates their evolution. This analysis will clear the route to an understanding of entropy evolution and will provide a platform on which to introduce control.

There is no unique way to define the symplectic invariants but the collection most commonly considered is the set of symplectic eigenvalues, expressible through Williamson’s theorem [24, 27]. This states that any $2n$-dimensional, real, symmetric, positive-definite matrix, for instance $\sigma$, may be decomposed as

[TABLE]

where

[TABLE]

and $\mathbb{I}_{2}$ is the $2\times 2$ identity matrix, $W>0$ and $S\in\operatorname{Sp}(2n,\mathbb{R})$. The elements of $W$ are unique up to reordering. Furthermore the uncertainty principle on $\sigma$ ensures that $\nu_{i}\geq 1$ [24]. The parameters $\nu_{i}$ are the symplectic eigenvalues which can also be calculated as moduli of the eigenvalues of $\Omega\sigma$.

Any function of symplectic invariants is also a symplectic invariant. Consider the $k$th elementary symmetric function of the eigenvalues $\lambda_{i}$ of some $m\times m$ matrix $X$, defined as [28]

[TABLE]

where the sum runs over all the possible $k$-subsets $E\in\mathcal{E}_{k}^{m}$ of the first $m$ natural integers. More explicitly,

[TABLE]

where $\mathbb{N}_{m}=\{1,\ldots,m\}$, $P(\cdot)$ denotes the power set and we put $\vartheta_{0}[X]:=1$. As an example, for a matrix $X$ with eigenvalues $\{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}\}$,

[TABLE]

Following Ref. [29] this allows for the definition of a new set of symplectic invariants of $\sigma$, given by $\vartheta_{2k}[\Omega\sigma]$ and related to the symplectic eigenvalues via

[TABLE]

The elementary symmetric functions of odd order vanish in this case because of the alternating signs of the eigenvalues of $\Omega\sigma$. Both the set of symplectic eigenvalues and this new set contain the entropic information of the state.

To further our analysis of these new invariants we recall some basic linear algebra [28]. Consider an $m\times m$ matrix $X$ and delete the same $m-r$ rows and columns. The remaining $r\times r$ submatrix is known as a principal submatrix of $X$. The determinant of a principal submatrix is known as a principal minor. Given the characteristic polynomial $\sum_{k=0}^{m}c_{k}\lambda^{m-k}$ of $X$, where $c_{0}=1$, we have that

[TABLE]

Noting that $c_{2k}[X]\equiv\vartheta_{2k}[X]$ we see that the study of principal minors and characteristic polynomial coefficients gives a new way to express and understand the symplectic invariants of a covariance matrix. In fact there is a way to recursively generate $c_{k}[X]$ using Fadeev-Le Verrier recursion.

Theorem 1** (Fadeev-Le Verrier recursion [30, 31, 32]).**

Let $X$ be an $m\times m$ real matrix. Let its characteristic polynomial be written $\det(X-\lambda\mathbb{I})=\sum_{k=0}^{m}c_{k}\lambda^{m-k}$ with $c_{0}=1$. It is possible to calculate the coefficients of the polynomial via the recursive formula,

[TABLE]

This recursive generation of $c_{k}[X]$ will be the guide to stripping Eq. (11) to only consider the evolution of the invariants. Before stating this evolution equation it is necessary to spend a little more time on notation, defining

[TABLE]

where

[TABLE]

This new object acts as a sort of reduced $\vartheta_{2k}[\Omega\sigma]$ where we remove the terms involving the $i$th symplectic eigenvalue. Again we define $\bar{\vartheta}_{0}^{i}[\Omega\sigma]:=1$. To illustrate this new function consider some covariance matrix $\sigma$ with symplectic eigenvalues $\{\nu_{1},\nu_{2},\nu_{3},\nu_{4}\}$. Here, we would have

[TABLE]

Using Fadeev-Le Verrier recursion to express the characteristic polynomial coefficients and their link to symplectic eigenvalues we are able to derive an evolution equation for these invariants under the open diffusive dynamical equation, Eq. (11).

Theorem 2**.**

For $\sigma$ evolving under $\dot{\sigma}=-\sigma+\chi\mathbb{I}$, the evolution of the symplectic invariants, defined by $\vartheta_{2k}[\Omega\sigma]$, obeys

[TABLE]

where $\sigma=SWS^{\intercal}$ by the Williamson decomposition and

[TABLE]

**Proof ** See A.

Eq. (24) directly reveals the evolution of the $n$ invariants under Eq. (11), meaning that we are able to focus on the parameters that interest us and ignoring the $2n^{2}$ remaining ones.

Open loop control refers to the enactment of predetermined operations without feedback on a quantum state, which would typically be unitary. The importance of the symplectic group is that it is the Gaussian analogue of the unitary group, in that the elements do not alter the entropic, or Gaussian, properties of the state. We consider symplectic control acting instantaneously at a given time in the evolution, allowing the alteration of $S$ on the right hand side of Eq. (24). Notice that this represents an impulsive action altering the state by a finite amount at a certain time and not a term generated by a Hamiltonian proportional to a step function in time (the $S$ in Eq. (24) is a time-dependent property of the evolving state). This is with the aim of optimising the rate of the change of $\vartheta_{2k}$, where the argument of $\vartheta_{2k}$ is assumed to be $\Omega\sigma$ unless otherwise stated.

Note that $S$ is defined via the Williamson decomposition of $\sigma$ and so both terms on the right hand side of Eq. (24) could be affected by its alteration. The first term, however, is a symplectic invariant and so no amount of manipulation will alter its value. However, the second term does vary with $S$ and so if, at a given point in the evolution, we alter it then this term will change. Maximising or minimising the value of the trace term is therefore the route towards locally optimising the rate of change of $\vartheta_{2k}$. In Sec. 5 we find that this optimisation will allow for the study of decoupled dynamics where each mode evolves independently. This drastically simplifies the analysis and allows for the optimal relaxation times for heating and cooling. We will show that squeezing is the key parameter that alters this term in Eq. (24).

Notice also that all possible initial Gaussian states are allowed in our treatment, since all physical covariance matrices admit a Williamson decomposition in terms of some symplectic $S$.

5 Heating and cooling

Eq. (24) shows that varying $S$ to alter the trace term is the required action to optimise the rate of change of $\vartheta_{2k}$. First we derive the symplectic matrix that will maximise or minimise this term and then we turn our attention to the dynamics that this optimisation invokes. This latter analysis will provide locally optimal heating and cooling times for the evolution of the channel under symplectic control.

5.1 Trace optimisation

To explore the optimisation of the trace term in Eq. (24) we first require an understanding of the singular value decomposition of symplectic matrices. Any symplectic matrix may be decomposed as $S=R_{1}ZR_{2}$ where $R_{1},R_{2}\in\operatorname{OSp}(2n,\mathbb{R})=\operatorname{O}(2n)\cap\operatorname{Sp}(2n,\mathbb{R})$ and $Z={\rm diag}(z_{1},1/z_{1},\ldots z_{n},1/z_{n})$. The elements of $\operatorname{OSp}(2n,\mathbb{R})$ correspond to ‘passive’ elements in the lab such as phase shifters and beam splitters. The $Z$ component encodes the squeezing element of the transformation. Any symplectic that is singular value decomposed such that $Z\neq\mathbb{I}$ is referred to as ‘active’. Using this decomposition will allow us to deconstruct and understand the trace term of Eq. (24).

For generality and brevity of notation we will consider a trace term of the form ${\rm Tr}[SYS^{\intercal}]$ where $Y=\bigoplus_{i=1}^{n}y_{i}\mathbb{I}_{2}$, with $y_{i}$ positive, and $S\in\operatorname{Sp}(2n,\mathbb{R})$. Note that the elements of $Z$ are unbounded and closely related to the energy required to enact the operation. Therefore the full optimisation will be done over the set $\overline{\operatorname{Sp}}(2n,\mathbb{R})$, defined as the restriction of the symplectic group to elements with maximum singular value $\overline{z}_{i}$ for each mode. This encapsulates the reasonable assumption that the energy that may be employed in each mode, rather than the sum of the energies, is bounded.

Proposition 1**.**

The supremum of the trace term over the set $\overline{\operatorname{Sp}}(2n,\mathbb{R})$ is given by

[TABLE]

where $\zeta^{\pm}_{z}:=\frac{z^{2}\pm 1/z^{2}}{2}$, $\bar{z}_{i}$ is the maximum local squeezing allowed on mode $i$, ordered such that $\bar{z}_{1}\geq\ldots\bar{z}_{n}$, and $y_{1}\geq\ldots\geq y_{n}$. The infimum is given by ${\rm Tr}[Y]$.

**Proof ** The singular value decomposition provides

[TABLE]

It is possible to change to a basis in which the symplectic form becomes $\Omega\to\begin{pmatrix}0_{n}&\mathbb{I}_{n}\\ -\mathbb{I}_{n}&0_{n}\end{pmatrix}$, where $0_{n}$ denotes the $n\times n$ zero matrix and $\mathbb{I}_{n}$ the $n\times n$ identity matrix. This is a commonly used basis in Gaussian state theory and is enacted with

[TABLE]

acting by similarity on the components of the trace term [16]. After enacting this we will then consider a further similarity transformation with

[TABLE]

This transforms symplectics into a basis which highlights the isomorphism between $\operatorname{OSp}(2n,\mathbb{R})$ and $\operatorname{U}(n)$, as will become clear. Acting on each matrix in the trace by $P$ and then $Q$ we find that the transformation of $Z^{2}$ gives

[TABLE]

where $\Gamma^{\pm}(z):=\operatorname{diag}(\zeta^{\pm}_{z_{1}},\ldots,\zeta^{\pm}_{z_{n}})$. $Y$ transforms as

[TABLE]

where $\Upsilon=\operatorname{diag}(y_{1},\ldots,y_{n})$ and $R_{2}$ transforms as

[TABLE]

where $U$ is some general unitary matrix. Note that $*$ here denote the complex and not the Hermitian conjugate. Altogether

[TABLE]

where $\alpha(z)$ is a vector of the diagonal elements of $\Gamma^{+}(z)$ and $\beta$ is a vector of the diagonal elements of $\Upsilon$. $P_{ij}=|U_{ij}|^{2}$ making it a general unistochastic matrix, which is a subset of the bistochastic matrices [33]. Bistochastic matrices are those that have have non-negative entries and whose rows and columns sum to 1. It will suffice to show that, for the supremum case (at $z$ fixed),

[TABLE]

and for the infimum case that

[TABLE]

because all permutations matrices are unistochastic [33]. Note that $\downarrow$ refers to rewriting the elements of the vector in descending order and the reverse for $\uparrow$. This is shown in Ref. [34] and the proof is reproduced in B. In the supremum case we further have to take the maximum over $z_{i}$ which is set at $\bar{z}_{i}$. Similarly, the infimum taken over $z_{i}$ is $z_{i}=1$. These two situations provide the final result stated in Eq. (26).

Thus we see that the maximisation of the trace term is found by setting $S={\bar{Z}}={\rm diag}({\bar{z}}_{1},1/{\bar{z}}_{1},\ldots{\bar{z}}_{n},1/{\bar{z}}_{n})$ for the maximum squeezing value. The minimum is obtained by setting $S=\mathbb{I}$. Note that neither of these provide the unique maximum or minimum. We could just as well maximise by setting $S=R{\bar{Z}}$ for some $R\in\operatorname{OSp}(2n,\mathbb{R})$ and ${\bar{Z}}$ maximum. Minimisation will occur for any $S=R$. Proposition 1 shows that the instantaneous control must either squeeze or unsqueeze depending on the desired effect.

5.2 Two-mode example

Let us now apply our findings to the iconic case of an initial two-mode squeezed state, such as the output of a degenerate parametric down conversion process, or the state that would result by driving an optomechanical set-up on a blue sideband, characterised by the following covariance matrix [14, 16]:

[TABLE]

with $\gamma\geq 1$. The parameter $\gamma$ allows us to range from pure states (for $\gamma=1$) to mixed states, with entropy growing for increasing $\gamma$. The parameter $r$ is sometimes referred to as the ‘two-mode squeezing parameter’ and, in a sense, quantifies the correlations between the two modes in this class of states. In practice, $r$ is determined by the strength of parametric interaction between the two modes and by the interaction time, if one assumes the state to be generated unitarily from the vacuum. In the following, we will assume the realistic value $r=0.4$.

To fix ideas, let us consider in what follows a loss rate of $100$ ${\rm KHz}$, and a $\chi$ parameter equal to $1.000013$, corresponding to the thermal noise at room temperature experienced by a mode of visible radiation at $450$ ${\rm THz}$ (a noisier situation could be construed but, as we shall see, this typical optical circumstances will suffice to illustrate our methods). For simplicity, we will quantify the evolution of the entropy of the relaxing states $\varrho_{t}$ through their purity $\mu={\rm Tr}[\varrho_{t}^{2}]$ which, for Gaussian states, turns out to be a function of the invariant $\vartheta_{4}[\Omega\sigma]$ alone (which is just the determinant of $\sigma$): $\mu=1/\sqrt{\vartheta_{4}[\Omega\sigma]}$ [14]. The steady-state purity with parameters set as above is $0.99997$.

The analytic optimisation expressed by Theorem 2 and Proposition 1 may now be applied to compare the uncontrolled lossy relaxation of the initial states given above with the relaxation of optimally adjusted states under instantaneous, symplectic control.

In the case of cooling, for $\gamma=2$ (such that the purity of the initial state is $0.25$), the optimal strategy to speed up relaxation is to undo the two-mode squeezing operation and let the vacuum state evolve through the channel. This allows one to attain a purity of $0.9$ after $30$ $\mu{\rm s}$, as opposed to the $35$ $\mu{\rm s}$ required by the uncontrolled evolution; similarly, a purity of $0.99$ is reached after $53$ $\mu{\rm s}$, as opposed to the $59$ $\mu{\rm s}$ of the uncontrolled evolution (clearly, the asymptotic purity is never perfectly reached, so one must resort to thresholds).

Let us now turn our attention to heating by setting $\gamma=1$ (perfectly pure initial state). Then, the fastest relaxation strategy would be to undo the two-mode squeezing of the initial state and apply maximum single-mode squeezing, as per Theorem 2 and Proposition 1. In order to avoid the introduction of arbitrary maximum squeezing parameters (that would, in a lab, be dictated by specific technical constraints), let us just contrast the uncontrolled evolution of the initial pure two-mode squeezed state with the evolution of an initial vacuum state, obtained by unitarily undoing the initial two-mode squeezing: This will provide us with a spectacular demonstration of the effect of squeezing on relaxation dynamics. In mathematical terms, the difference between the two cases is captured by the trace term in Eq. (24), for $n=2$ and $k=2$, which corresponds to the evolution of the determinant of a two-mode covariance matrix. In the case of heating, relaxation is hastened by an increase in the trace term: hence, whilst an initial vacuum state takes about $40$ $\mu{\rm s}$ to relax (to the significant digits reported above), the initial two-mode squeezed state takes as little as $4\times 10^{-4}$ ${\mu}{\rm s}$! It can be shown that an instantaneous symplectic applied at this point would then be capable of stabilising not only the purity, but the whole covariance matrix (and hence the state) to the steady-state values. Such a symplectic is the one that enacts the Williamson decomposition of the covariance matrix and can be efficiently evaluated with standard methods. To give a quantitative idea of the effect of such control on the fly, let us add that the determinant of the initial two-mode squeezed state would, without it, overshoot the asymptotic determinant after $4\times 10^{-4}$ ${\mu}{\rm s}$, reach a maximum (corresponding, in this instance, to a minimum purity of $0.855633$ at $7$ ${\mu}{\rm s}$), and then readjust to the steady state after a total time of about $140$ ${\mu}{\rm s}$.

5.3 Decoupling

Proposition 1, in conjunction with Eq. (24), allows us to understand the dynamics of $\vartheta_{2k}$ when $S$ is controlled to be optimal at a given moment of time. In the optimal limit, either maximal or minimal, the value of $S$ must be either ${\bar{Z}}$ or $\mathbb{I}$. Using Williamson’s theorem, as stated earlier, Eq. (11) may be rewritten as

[TABLE]

If we enact a control to set $S=Z$ or $S=\mathbb{I}$ at a given time then this equation will decouple so that it suffices to consider single mode evolution, as explored in Ref. [21]. For the single mode case Eq. (24) becomes

[TABLE]

where $\nu_{i}$ refers to the symplectic eigenvalue of the decoupled mode $i$, $z_{i}$ is its squeezing value obtained from the singular value decomposition and $\zeta_{z_{i}}^{+}$ is defined in Proposition 1. Note that in the optimal maximal limit $z_{i}\to\overline{z}_{i}$ and in the minimal limit $z_{i}=1$. We may talk about maximising or minimising the rate of change of $\vartheta_{2}$ by considering different limits of $z_{i}$. Recalling that $\vartheta_{2}=\nu_{i}^{2}$ we see that Eq. (38) is equivalent to

[TABLE]

which has solution

[TABLE]

where $\nu_{i0}$ is the initial value of the decoupled mode $\nu_{i}$. The value of $\nu_{i}$ directly determines the entropy of the Gaussian state which for a single mode is defined as [35]

[TABLE]

$\kappa$ monotonically increases with $\nu_{i}$ and so, for a single mode, $\nu_{i}$ is a good entropy measure, being equal to one for pure states and greater than one for mixed states. Gaussian states have an associated temperature encoded in $\beta$, and each mode has an associated mode frequency $\omega_{i}$. The value of $\nu_{i}$ depends on the product of these two parameters via [36]

[TABLE]

If we fix $\omega_{i}$ then $\nu_{i}$ monotonically increases with the temperature, thus we refer to rising $\nu_{i}$ as heating and lowering $\nu_{i}$ as cooling.

The fixed point of the evolution in Eq. (12) is at $\nu_{i}=\chi$. The time to reach this point diverges, hence we must fix a distance within which we are satisfied that we have arrived close enough to the target. Take this to be

[TABLE]

The case $\nu_{i0}<\chi$ means that the channel is heating. The optimal strategy in this case is to increase the squeezing of each of the decoupled modes as much as possible. In other words, we choose $z_{i}=\bar{z}_{i}$ for all modes and the minimum amount of time for all to come within distance $\epsilon$ of the fixed point is

[TABLE]

Note that this time is in general finite and goes to zero for large values of $\bar{z}_{i}$ (i.e., $\zeta_{\bar{z}_{i}}^{+}$ large). Furthermore it is possible to set $\epsilon=0$ and so we are able to reach exactly the fixed point in an infinitesimal amount of time for $\bar{z}_{i}\to\infty$ and a finite time for finite $\bar{z}_{i}$.

The case $\nu_{i0}>\chi$ means that the channel is cooling. In this case we optimally set $z_{i}=1$ for all the decoupled modes. The minimum amount of time to come within distance $\epsilon$ of the fixed point of the dynamics is

[TABLE]

$T_{\textit{cool}}$ diverges for $\epsilon=0$ and so to reach the true fixed point still requires an infinite amount of time. However, in most cases it suffices to cause the state to come arbitrarily close to the thermal state of the bath. We see that in order to minimise the time to achieve this it is necessary to reduce the mode squeezing to zero, explicitly we wish to reduce the squeezing measure $\xi=\max{\rm eig}[S^{\intercal}S]-1$ to zero. This ‘unsqueezing’ will need to be enacted to optimise cooling.

As it is shown in C, our analysis reproduces the results of Ref. [21] for the special case of an unsqueezed bath and a single mode system. Referring to equations from Ref. [21] with primes it is easy to show that $T_{\textit{heat}}$ coincides with $T^{\mathrm{heat}}_{\mathrm{fast}}$ (Eq. (67’)) and $T_{\textit{cool}}$ with $T^{\mathrm{cool}}_{\mathrm{fast}}$ (Eq. (63’)) provided that $\epsilon\simeq 2\sqrt{\chi^{2}-1}\sqrt{\epsilon^{\prime}}$.

Furthermore we note that optimal control can prevent a cooling channel from ever decreasing the temperature of the state if it is sufficiently, periodically squeezed. This can be seen by noting that $\nu_{i}$ increases under Eq. (39) when $\nu_{i}<\chi\zeta_{z_{i}}^{+}$ even when $\nu_{i}>\chi$, as seen in Fig. 1.

We would like to compare $T_{\textit{heat}}$ and $T_{\textit{cool}}$ with the free relaxation time of the state, in the absence of control. Computing this for a system of $n$ modes would require solving Eq. (12) with some analogue of Eq. (43), which is not easy to achieve analytically. Ref. [21], however, provides results for the relaxation of single-mode Gaussian states and so a comparison can be made in this regime. As shown in C, the free decoherence time for a single mode can be written as:

[TABLE]

where $z_{0}$ is the initial squeezing of the mode. For the case when a single mode is involved, therefore, we can draw exactly the same conclusions as in Ref. [21]. In particular, for vanishing tolerable errors, we get that the advantage of quantum control appears negligible if the channel is cooling while it allows for an exponential speed-up with respect to the uncontrolled dynamics for heating channels. This former statement is somehow reminiscent of the third law of thermodynamics according to which zero temperature is unattainable in finite time [37]. One can also study a measure of performance of the control procedure for the worst case scenario of possible initial conditions and obtain the same results discussed in Ref. [21].

Therefore we see that squeezing and heating are intrinsically related. Squeezing itself does not affect the temperature of the state but under a lossy channel it is the key factor in increasing or decreasing its rate of change, here shown for any number of modes.

6 Conclusion and Discussion

The number of degrees of freedom for an $n$-mode Gaussian state scales polynomially in $n$. Therefore when we consider the evolution of such a state it is important that we isolate the parameters that are most relevant. The information encoded in the $n$ symplectically invariant parameters allows one to extract the key property of entropy. In this paper we derived an equation for the evolution of such invariants and calculated the locally optimal rates for heating and cooling under a lossy channel. We have shown that among the time-optimal control schemes, one involves decoupling the dynamics into single modes. From here optimal cooling and heating come about by respectively unsqueezing and squeezing the decoupled modes. As a result of the decoupling, one can extend the single-mode analysis of [21] to any number of degrees of freedom, which represents our main result. We have demonstrated that on the one hand heating can be accomplished arbitrarily well via optimal symplectic control in an arbitrarily short time if enough squeezing is allowed while, on the other hand, controlled cooling cannot be achieved equally well and fast as heating.

This analysis rested on the derivation of Eq. (24) but the techniques required can also provide insight into a property that is not a symplectic invariant: entanglement. The positivity of the partial transpose is necessary for the separability of a bipartite $p+q$ mode state, and sufficient if either $p$ or $q$ is equal to $1$ [38, 14]. Separability is equivalent to the condition that $\tilde{\sigma}=T\sigma T$ obeys the uncertainty relation

[TABLE]

where $T=\bigoplus_{1}^{p}\sigma_{z}\oplus\bigoplus_{1}^{q}\mathbb{I}_{2}$, where the Pauli matrix $\sigma_{z}={\rm diag}(1,-1)$. This condition is equivalent to the smallest symplectic eigenvalue of $\tilde{\sigma}$ being greater than or equal to 1 which in turn gives the separability condition [38]

[TABLE]

$T$ acts by similarity and is symmetric meaning that $\tilde{\sigma}$ remains positive definite, implying that it also has a Williamson decomposition: $\tilde{\sigma}=\tilde{S}\tilde{W}\tilde{S}^{\intercal}$. Furthermore it obeys $\dot{\tilde{\sigma}}=T\dot{\sigma}T=-\tilde{\sigma}+\chi\mathbb{I}$ since $T$ is time invariant. This allows us to mimic the full analysis of A providing

[TABLE]

where $p+q=n$ and $\tilde{V}_{k}$ has the same definition as before except for the symplectic eigenvalues of $\tilde{\sigma}$. This provides us with an evolution equation for an entanglement measure of the $p+q$ mode system, directly from the analysis employed to derive Eq. (24). The difference here, however is that $\vartheta_{2k}[\Omega\tilde{\sigma}]$ is not a symplectic invariant and so further analysis will require more sophistication.

This paper brought together techniques from linear algebra to explore the evolution of multimode Gaussian states evolving under lossy channels. As exemplified by the calculation concerning entanglement, sketched above, such techniques may be extended to provide wider analytical insight. We restricted ourselves to the specific question of finding locally optimal cooling, or heating, for states undergoing lossy channel evolution. This has been answered and shown to be intrinsically related to squeezing. We showed that the study of locally optimal trajectories can be explored using decoupled dynamics and so the single mode case provides the lower bound on the relaxation time.

Acknowledgments

AS acknowledges financial support from EPSRC through the grant EP/K026267/1.

Appendix A Evolution of symplectic invariants

A.1 Recursive formulas

Reiterating Eq. (11) the evolution of our system obeys

[TABLE]

Using this and the Taylor expansion in $\rm{d}t$ we obtain

[TABLE]

where $F:=(1-{\rm d}t)\Omega\sigma$ and $G:=\chi\Omega$.

Lemma 1**.**

The Taylor expansion of the following trace term to first order is

[TABLE]

**Proof ** Expanding $(F+G{\rm d}t)^{2k-i}$ to first order in ${\rm d}t$ we obtain a single term of the form $F^{2k-i}$ and $(2k-i)$ terms of the form $F^{a}G{\rm d}tF^{b}$ where $a+b=2k-i-1$. The cyclicity of the trace allows us to reorder these elements to obtain ${\rm Tr}[F^{2k-i}]+(2k-i){\rm d}t{\rm Tr}[F^{2k-i-1}G]+o({\rm d}t)$. Now it remains to show that to first order this expression is zero for $i$ odd.

It suffices to show that

[TABLE]

To prove Eq. (53) we use the invariance of the trace under cyclic permutations and transposes giving

[TABLE]

where we used $\Omega^{\intercal}=-\Omega$. Eq. (54) is found using a similar argument. Putting these together we prove the proposition.

Recall the recursive definition for the coefficients of the characteristic function given in Sec. 4:

[TABLE]

Recalling that $c_{2k}[\cdot]\equiv\vartheta_{2k}[\cdot]$ we change notation to start considering these symmetric functions. Using the recursive formula and Lemma 1 we see that

[TABLE]

Lemma 2**.**

Taylor expanding $\vartheta_{2k}[F+G{\rm d}t]$ we arrive at

[TABLE]

**Proof ** From Eq. (56) we can explicitly show that

[TABLE]

where we define $\gamma:=\Omega\sigma$ for brevity in the proceeding proof and where we have also used Lemma 1. From here we will proceed with an inductive proof. We assume that Eq. (57) holds for some $k$ and then show that if this is true then it holds for $k+1$.

Using Eq. (56) we may write the expansion out for $k+1$ and then use Lemma 1 and the definition of $F$ to arrive at

[TABLE]

Substituting Eq. (57) into Eq. (59) we extract the first two terms that look like they would fulfill the proof plus a final one that we would hence like to show is zero:

[TABLE]

Dropping $(2(k+1))^{-1}$ we proceed to examine the final piece, referring to it as $L$, and rewriting it as

[TABLE]

Note that we have relabelled $k$ as $k-1$ to shorten the expression but it will not alter the analysis. The ${\rm d}t$ at the front reminds us that everything should be expanded to zeroeth order inside the sum. To prove the lemma it is necessary to show that $L\equiv 0$.

Expanding $\vartheta_{2i}[F]$ to introduce another sum we arrive at

[TABLE]

From here note that for a general sum with elements $Y_{ij}$ we have

[TABLE]

where $i^{\prime}=i-(j+1)$ and $i^{\prime\prime}=k-i^{\prime}-j-2$. The first equality of Eq. (63) can be seen with observation. The second involves a redefinition of the sums where we split them into two halves and then redefine the labels such that one is descending whilst the other ascends. When such a summation redefinition is applied to Eq. (62) it will be clear that $L\equiv 0$.

Thus we prove that if $\vartheta_{2k}[F+G{\rm d}t]$ is given in Eq. (57) then this also holds for $k\to k+1$. From Eq. (58) we see that it is true for $k=1$ and so, inductively it is true for all $k$. To write it in the form stated one must replace $\gamma$ with $\Omega\sigma$.

Lemma 3**.**

Taylor expanding $\vartheta_{2k}[F]$ using $F:=(1-{\rm d}t)\Omega\sigma$ we find that

[TABLE]

**Proof ** Expanding out the recursive formula and again defining $\gamma:=\Omega\sigma$ we get a product of sums of the form

[TABLE]

By only keeping terms that are less than second order in ${\rm d}t$ we get a smaller sum

[TABLE]

$X$ consists of the remaining terms which come in pairs. Take for instance the first pair which is generated by choosing the $+2i_{1}{\rm d}t$ coefficient in the first line of Eq. (65), with everything else at zeroeth order, and secondly the $-2i_{1}{\rm d}t$ coefficient in the second line, with everything else at zeroeth order. The pairs will each cancel to become zero. The final piece comes without a partner but has coefficient $i_{k}=0$, and so does not contribute. Therefore $X\equiv 0$ and the lemma is proven.

Lemma 4**.**

Using the previous two Taylor expansions we may write the rate of change of $\vartheta_{2k}$ as

[TABLE]

**Proof ** The first Taylor expansion came from Lemma 2 stating that

[TABLE]

Using Lemma 3 and the same reasoning as in Eq. (51) we rewrite the above as

[TABLE]

Dividing through by ${\rm d}t$ we prove the proposition, recalling that $G:=\chi\Omega$.

A.2 Telescoping the series

Eq. (67) provides a differential equation describing the rate of change of $\vartheta_{2k}$. We now work towards a simplification of this equation using Williamson’s theorem [27, 24] and noticing that the series ‘telescopes’ to provide a simpler form.

Given Williamson decomposition, $\sigma=SWS^{\intercal}$, as described earlier, and the symplectic property: $S^{\intercal}\Omega S=\Omega$, we may rewrite the trace term that appears in the sum of Eq. (67):

[TABLE]

where we used $(W\Omega)^{2}=W^{2}\Omega^{2}$, $\Omega^{2}=-\mathbb{I}$ as well as cyclic properties of the trace. Thus Eq. (67) becomes

[TABLE]

Using definitions given in Sec. 4 we may now prove the theorem stated in the main body:

Theorem** (2).**

For $\sigma$ evolving under $\dot{\sigma}=-\sigma+\chi\mathbb{I}$, the evolution of the symplectic invariants, defined by $\vartheta_{2k}[\Omega\sigma]$, obeys

[TABLE]

where

[TABLE]

**Proof ** Beginning with the final term of Eq. (68) and dropping the $\chi$ factor we may take the sum inside the trace to give

[TABLE]

which allows us to consider the internal sum before tracing. Recalling that the diagonal elements of $W$ occur in pairs we only need to consider $n$ symplectic eigenvalues denoted $\nu_{q}$, $q=1,\ldots,n$. The $i$th term of the sum in Eq. (71) takes the form

[TABLE]

and the $i+1$th takes the form

[TABLE]

Recall that $\bar{\vartheta}_{a}^{b}$ is a sort of reduced symplectic invariant which is similar to $\vartheta_{a}$ but where we remove all terms involving the $b$th symplectic eigenvalue. By splitting up Eq. (71) to consider each $\nu_{q}$ individually and also by spltting each $i$th term we see that telescoping is going to occur. This is the situation when the latter piece of the $i$th term cancels the former term of the $i+1$th term. Noting that this cancellation is going to occur between Eq. (72) and Eq. (73) we should be left with the very first and last pieces of the entire series. The first piece is equal to zero and so we are just left with the final term

[TABLE]

and therefore the expressions of Eq. (71) are equal to ${\rm Tr}\left[SV_{k}S^{\intercal}\right]$ where

[TABLE]

Appendix B Bistochastic limits

Here, we prove an optimisation result regarding bistochastic matrices. We begin with the reiteration of the notation that, for some vector $v$, $v^{\downarrow}$ denotes a new vector of elements of $v$ written in descending order. $v^{\uparrow}$ is similar but for ascending order.

Lemma 5**.**

Let $\alpha$ and $\beta$ be two real vectors of length $m$. The supremum of the following inner product is given by

[TABLE]

and the infimum by

[TABLE]

**Proof ** Taken from Ref. [34]. We begin with a consideration of the supremum and then turn our attention to the infimum. The permutation matrices form a subset of the bistochastic matrices and so we may alter $X$ to make $\alpha$ and $\beta$ descending. Then we fix $X$ to some specific bistochastic matrix to define

[TABLE]

where $a_{i}$ and $b_{i}$ are the elements of $\alpha^{\downarrow}$ and $\beta^{\downarrow}$ respectively. Consider $X\neq\mathbb{I}$. Let $k$ be the smallest index $i$ such that $X_{ii}\neq 1$. Note that for $i<k$, $X_{ii}=1$ and therefore $X_{ij}=0$ if $i<k$ and $i\neq j$, or if $j<k$ and $i\neq j$. Since $X_{kk}<1$, then for some $l>k$, $X_{kl}>0$. Likewise, for some $p>k$, $X_{pk}>0$. These imply that $X_{pl}\neq 1$.

The inequalities above mean that we can choose $\epsilon>0$ such that the matrix $X^{\prime}$ is bistochastic where

[TABLE]

and $X^{\prime}_{ij}=X_{ij}$ in all other cases. Now define

[TABLE]

Recalling that $l>k$ and $p>k$, so that $a_{k}<a_{p}$ and $b_{k}<b_{l}$,

[TABLE]

which means that the term $\sum a_{i}b_{j}X_{ij}$ is not decreased. $\epsilon$ may be chosen to reduce an off-diagonal term in $X$ to zero without affecting the bistochasticity of $X$ and without decreasing the value of $\chi^{\prime}-\chi$. As this process is iterated $X$ may be brought to identity $\mathbb{I}$ without decreasing $\chi$, achieving its maximum value.

The argument for the infimum is completely analogous except we begin by considering

[TABLE]

where $c_{i}$ are the elements of $\alpha^{\uparrow}$. The analysis is similar except that we will be considering $\chi^{\prime}-\chi\leq 0$, and so transforming $X$ into $\mathbb{I}$ does not increase the value of $\chi$. Thus we find $X=\mathbb{I}$ again but for the situation where the vectors are oppositely ordered.

Appendix C The case of a single mode

Here we compare $T_{\textit{heat}}$ and $T_{\textit{cool}}$ with $T^{\mathrm{heat}}_{\mathrm{fast}}$, $T^{\mathrm{cool}}_{\mathrm{fast}}$ and $T_{\mathrm{free}}$ derived in Ref. [21] for the case of a single mode coupled to an unsqueezed bath. For clarity, equations quoted from Ref. [21] will be primed. Both papers begin with Eq. (8) which allows us to identify

[TABLE]

with symbols from Ref. [21] on the left and ours on the right. Eqs. (38’)-(39’) show the fixed point of the dynamics to be

[TABLE]

Using Eq. (33’) and noting that $\mathrm{Tr}(\sigma)=2\nu\zeta_{z}^{+}$ we identify

[TABLE]

and maximal squeezing in $\cosh 2r_{M}\leftrightarrow\zeta_{\bar{z}}^{+}$, while for the initial conditions we identify

[TABLE]

where $z_{\textit{0}}$ is the initial squeezing value. Ref. [21] uses the fidelity Eqs. (50’)-(51’) for the tolerable error while we use Eq. (43). To avoid confusion, we rename $\epsilon^{\prime}$ the allowed error in Eq. (50’). In particular, we can rewrite Eqs. (50’) and (66’) for the heating and cooling cases, respectively, as $\mu_{\mathrm{Th/Tc}}\leftrightarrow\frac{1}{\chi}\left[1(+/-)2\frac{\sqrt{\epsilon^{\prime}}}{\chi}\sqrt{\chi^{2}-1}\right]$. With all of the above taken into account, it is then very easy to verify that Eq. (67’) and Eq. (63’) coincide with our optimal times for decoupled modes, i.e. that

[TABLE]

provided that $\epsilon\simeq 2\sqrt{\chi^{2}-1}\sqrt{\epsilon^{\prime}}$. Finally, the free decoherence time Eq. (54’) is translated as Eq. (46).