Thermodynamics and optimal protocols of multidimensional quadratic Brownian systems

We characterize finite-time thermodynamic processes of multidimensional quadratic overdamped systems. Analytic expressions are provided for heat, work, and dissipation for any evolution of the system covariance matrix. The Bures-Wasserstein metric between covariance matrices naturally emerges as the local quantifier of dissipation. General principles of how to apply these geometric tools to identify optimal protocols are discussed. Focusing on the relevant slow-driving limit, we show how these results can be used to analyze cases in which the experimental control over the system is partial.


Introduction
The minimisation of dissipation is a central goal in finite-time thermodynamics [1][2][3]. In most applications, one is interested in finding the optimal time variation of some control parameters, e.g. magnetic or electric fields, in order to achieve a desired task while minimising the amount of energy dissipated to the environment. Such tasks could range from the design of a cycle for a thermal engine [4,5] to the erasure of information in an information processing device [6][7][8]. Finding optimal protocols in finite time is however often a very challenging task, as it requires a functional optimisation over all possible paths in the control parameter space, as well as a perfect understanding of the non-equilibrium dynamics resulting from such control. In the regime of small mesoscopic systems, remarkable progress on this topic has been achieved in the last decades with the development of the field of stochastic thermodynamics [4,[9][10][11][12][13][14]. Optimal drivings are nowadays known for overdamped [15][16][17][18][19] and underdamped systems [20,21], as well as driven single-level quantum dots [22]. However, such explicit solutions only exist for one-dimensional systems and are, in general, computationally hard to scale up.
Other solutions are known for situations in which the control parameter varies slowly compared to the system relaxation time, as the optimisation admits a geometric formulation [23][24][25][26][27][28][29][30] and the problem considerably simplifies. Indeed, the space of control parameters can be endowed with a Riemannian metric in such a way that geodesic paths correspond to minimally dissipative thermodynamic processes. While the geodesic equations might be hard to solve, the important realisation is that the number of coupled equations is given by the number of control parameters, and independent of the size of the system of interest (by comparison, a full out-of-equilibrium solution of the dynamics needs a number of equations that scales exponentially with the number of components of the system). This has enabled finding optimal driving protocols in such regime for complex systems such as a two dimensional Ising model [31,32], nanomagnets [33], and quantum spin chains [34]. Optimal protocols for different classes of slowly driven heat engines have also been developed by such a geometric approach [30,[35][36][37][38][39][40][41][42]. Besides the slow driving regime, the optimization problem can also be simplified in the opposite, fast-driving, regime [43][44][45][46].
Beyond the slow or fast driving limits, a general connection has been established between optimal transport and minimally-dissipative thermodynamic protocols in the overdamped limit [47,48]. This connection was recently exploited to show that the minimal dissipation in any process governed by a Langevin equation is directly related to the L 2 -Wasserstein distance between initial and final states [49] (see also [50][51][52]). However in general full control on the system's Hamiltonian is needed to saturate this bound. To address the relevant case of partial experimental control, one therefore requires a different approach, that is able to quantify dissipation on non-optimal trajectories.
In this paper, we study thermodynamic transformations for many-body quadratic overdamped systems. We derive general expressions for the flux of work and heat, and we show that the dissipation is governed by the Bures-Wasserstein (BW) distance between covariance matrices, which coincides with the L 2 -Wasserstein distance between Gaussian distributions of [49][50][51][52]. Our derivation allows for a direct generalisation of the well-known single-body case, studied by Schmiedl and Seifert for a single-particle overdamped system [15], as well as new insights on the form of optimal drivings. In particular, we provide an integral analytic expression for the dissipation valid for any response trajectory of the system, not necessarily minimally dissipating. This also naturally enables the study of partial control. That is, the situation where the limited number of control parameters does not allow for exploring the whole space of states, so that the fundamental lower bounds of [47,49] might not be reachable. This is a common scenario in complex systems, where experimentally only a few degrees of freedom are controllable. In order to illustrate application of our results, and to show the difference between partial and global control, we analyse a system of two interacting particles and a particle confined in a 2-dimensional squeezing potential with different control limitations.

Model: Many-body overdamped spring
We consider a system of N overdamped Brownian particles described by the position vector x and mutually interacting via the time-dependent potential or, equivalently, via the force field F (x) = −∇V (x) = −Kx (when possible, we omit writing the time argument from now on). Each particle i = 1, . . . , N might have a different number of degrees of freedom The potential (1) accounts for both self-energy of the individual particles and interactions between the particles. The stiffness matrix K is symmetric and positive definite K ≥ 0 (that is, the potential is confining). Assuming that all the particles have the same friction coefficient γ, the system dynamics obeys the set of Langevin equations [53] where the Gaussian noise η obeys η = 0, its components η i (t)η j (t ) = δ ij δ(t−t ), and T is the temperature of the thermal environment, which we assume is fixed throughout (isothermal). From now on, we will use natural units in which γ = 1, k B = 1. For the general case in which different particles have different friction coefficients, see Sec. 4. Departing from an arbitrary normalized initial distribution, the state of the system at time t is represented by a Gaussian probability density function (PDF) [53] ‡ (3) ‡ The PDF describing the state under the dynamics (2) is always Gaussian given Gaussian initial conditions, or after a negligible initial transient [53].
Here Σ(t) = xx (t) denotes the covariance matrix at time t. The PDF (3) has zero mean x (t) = 0 (see Sec. 4 and App. Appendix B for the more general case).
The distribution p(x, t) is therefore defined by its covariance matrix, whose dynamicsΣ = xẋ + ẋx can be obtained from Eq. (2) and reads [53] where implicitly T = T 1. In case the response dynamics Σ(t) is given, and one wants to drive the potential K(t) accordingly (i.e. K(t) is the control protocol generating the response dynamics Σ(t)), the solution of Eq. (4) for K(t) is Notice that instantaneous quenches of K(t) can be added at the beginning and at the end of the protocol, without affecting the dynamics of Σ(t). For example, to end the transformation in equilibrium, one can add a final quench to K = T Σ −1 .
Remark. We stress here that the explicit evaluation of Eq. (5) can be performed analytically. To be consistent with the notation, throughout the text we use the operator I(A, B) = ∞ 0 dν e −νA Be −νA expressed in its integral matrix form; at the same time, in the basis that diagonalizes A, i.e. A ij = δ ij a i , the components of this operator can be easily expressed as

Thermodynamics of quadratic systems
The average energy of a system described by a multidimensional probability distribution (3) in the potential (1) reads where Tr[AB] = ij A ij B ji . The variation of energy is split canonically in a work contribution, originating in the variation of the external potential, and a heat contribution, stemming from the evolution of the system induced by the dissipative environment [4,54]. I.e. the work (W ) and heat (Q) fluxes entering the system are defined asẆ In a similar fashion to the seminal work by Schmiedl and Seifert [15] we can write the work input of a finite time transformation of duration τ as where in the second equality we integrated by parts and used Eq. (5). In the following we will indicate as Q(τ ) − Q(0) := ∆Q the variation of any quantity Q during a transformation. Given that 1 2 log det Σ| τ 0 = − dx p ln p | τ 0 = ∆S is the variation of Von Neumann entropy of the system (3), it is possible to rewrite Eq. (8) as This expression identifies the dissipated work, W irr , of an arbitrary response trajectory Σ(t). This is our first main result. The above derivation represents a natural multidimensional generalisation of the one-dimensional result of [15]. The irreversible work (9) turns out to be the integral in time of a quadratic form that coincides with the Bures-Wasserstein (BW) metric on positive-definite matrices [55,56].
with the latter being the infinitesimal BW squared distance. This metric has been intensely studied as it appears in problems of statistical inference and metrology in quantum information [55,[57][58][59], as well as in the theory of optimal transport [56,60]. For fixed endpoints, the lower bound for W irr is obtained for the response trajectorȳ Σ(t) that minimizes the integral of the quadratic form in Eq. (9). That is where the BW-geodesic length between the initial and final points Σ(0) = Σ 1 , Σ(τ ) = Σ 2 , is given by (see Appendix A, or [55,56]) The corresponding geodesic, i.e. the minimally dissipating response trajectory, is given by (cf. Appendix A) with s = t/τ and thus 0 ≤ s ≤ 1 independently on the total duration t of the protocol. The appearance of the distance (13) in (12) is no coincidence: it was realised recently [49] that the optimal transport problem is connected to the irreversible entropy production in diffusive dynamics, and its value is minimized by the L 2 -Wasserstein distance between the initial and final distributions [47,49,61]. In the case of Gaussian distributions, the L 2 -Wasserstein distance coincides with the above BW distance between the covariance matrices (13).
Besides it being straightforward, one key advantage of our derivation is that expression (9) is valid for any response trajectory and allows to compute W irr also when the transformation does not saturate the lower bound (12). In particular, it can be used for the realistic case of partial experimental control, when the system is constrained to explore only a subset of the distributions space (see the following Sec. 3 and examples in Sec. 5).

Total control vs partial control.
In experiments, the system is typically controlled by varying K(t). The optimal control parameter protocolK corresponding to the geodesic (14) is determined by substitutinḡ Σ into (5). Perfect implementation ofK would then saturate the minimal dissipation bound (13). However, this assumes thatK is experimentally feasible. This might not be the case in general. Performing the minimization over a restricted region of control parameters limits the system response to a submanifold of allowed states. In general, this results in a case-dependent minimum value strictly larger than the global minimum, W irr > D 2 BW /τ , e.g., see Example 5.1 below. In other cases, the initial and final point of the transformation might not even be fixed (e.g., when optimizing the performance of a cyclic heat engine). To show the consequences of fixed/unfixed boundary states Σ(0) and Σ(τ ), consider that the variationΣ =Σ d +Σ r can be divided into a diagonal contribution and a nondiagonal, rotating contribution. That is, given Σ . From Eqs. (9) and (10) we know thatẆ irr = g Σ (Σ,Σ). It is easy to check that g Σ (Σ d ,Σ r ) = 0 which implies that the irreversible work naturally decouples into a diagonal and a rotating contribution:Ẇ irr are positive, which means that the dissipation generated in a non-commuting transformation for Σ (Ẇ In the fully commuting case, it is easy to see that Eqs. (6-9) simplify and we recover (K being diagonal in the same basis of Σ, with eigenvalues k i ) When reduced to a single mode, this is exactly the result found by Schmiedl and Seifert in [15], which we thus see being extensive in the eigenmodes ω i of Σ: that is, all the modes {ω i , k i } can be treated as effectively independent in the commuting case. As an instance of a transformation with fixed boundaries that force non-commutation, see Example 5.2. As mentioned above, controlling the potential (1) in time defines the evolution of the state via the dynamical equation (4). Conversely, a given response trajectory Σ(t) is translated to its generating control K(t) through Eq. (5). This means that fixing the boundary values of Σ is non-trivially related to fixing boundary controls. The results of optimisation thus strongly depend on the imposed constraints [62]. At the same time, for the purpose of typical applications to isothermal processes (cf. Sec. 5), in which the goal is to minimize work dissipation, we can consider the slow-driving limit of the dynamics [63]. In this limit the potential K(t) is modified slowly, more precisely we assumeK ∼ 1/τ with τ much larger than the relaxation timescale of the system τ γ/|K|, and it is sufficient to solve the dynamical equation (4) up to the first order in 1/τ . The zeroth order corresponds to the quasistatic limit τ → ∞. This allows to expand any state-dependent quantity Q around its equilibrium value Q (0) , keeping only the leading correction term Q (1) ∼ O(1/τ ). In our specific setting, the covariance matrix can be expanded as with (4) and (5)). However, the irreversible work (9) is already of order O(1/τ ). To express the dissipation in the slow regime, it is therefore sufficient to substitute Σ (0) in (9). In other words, we observe that Observation 2 In the slow-driving limit, controlling the inverse stiffness matrix T K −1 (t) of the potential is equivalent to directly steering the covariance matrix Σ(t). The irreversible work in the slow-driving limit therefore reads

Generalizations
In the previous section, we have focused on the paradigmatic case of the potential (1) and density distributions (3) centered around x = 0, and a particle-independent friction coefficient. Nevertheless, the obtained results are fully extendable also when removing such assumptions. First, in Appendix B we solve the general case of a quadratic potential with timedependent center z(t), i.e. V (x, t) = 1 2 (x − z(t))K(t)(x − z(t)). As the system is in general driven out of equilibrium, the center of the potential does not necessarily coincide with the average particle position, x ≡ ξ(t) = z(t), and the irreversible work gains an additional contribution (see details in Appendix B). Focusing on the limit of slow driving, it can be expressed as Similarly to (13), the lower bound for W slow irr is in this case Observations 1 and 2 from Section 3 remain valid: if possible, rotations of the covariance matrix thus should be avoided; the state variables in the expression (19) can be substituted by their equilibrium values (ξ, Σ) = (z, T K −1 ). Moreover, in the same limit, ξ(t) − z(t) ∼ O(τ −1 ) while the associated correction to quasistatic ∆E and ∆S is negligible O(τ −2 ) (cf. Appendix B). This implies that moving the trapξ = 0 only contributes to the dissipation (20) and should therefore be avoided when possible, in the same spirit as Observation 1. Second, we comment on the generalization to systems where different particles have different friction coefficients γ i . In such a case, the Langevin equations (2) become, in components, Notice that some of the γ i might refer to different degrees of freedom of the same particle. The white noises η i are mutually uncorrelated. We define y i ≡ √ γ i x i and rewrite the Langevin equations aṡ where the transformed stiffness matrix K ij = K ij √ γ i γ j is still symmetric and positive definite. At the same time, the covariance matrix of the y variable, Σ = yy , satisfies Σ ij = √ γ i Σ ij √ γ i . Finally the energy of the system is given by and similarlyQ = 1 2 Tr[Σ K ] andẆ = 1 2 Tr[Σ K ]. Given the formal equivalence between Eqs. (2,6,7) and (22,23) above, we see that the problem is equivalently mapped to the case with fixed γ i = γ, ∀i.
Finally, throughout the paper we assumed a fixed temperature T . At the same time, the expressions for energy (6), heat and work (7), as well as ∆S = 1 2 ∆ log det Σ do not intrinsically depend on the temperature T . We can therefore relax such assumption and admit a time-dependent temperature T (t) [64,65]. In such case the definition of irreversible work becomes S irr being the irreversible entropy production. From the derivation in Section 3 we get the same expression W irr = 1 2 τ 0 dt g Σ (Σ,Σ), as well as the validity of all the above observations and generalizations.

Applications
Here, we present two examples of application of the formalism, results and observations introduced above. stiffnesses k x = k y , and an additional harmonic interaction, k int , between them. A transformation is performed to modify the local traps stiffnesses. (b) Plots show the minimum values of the dissipated work W irr for transformations from k x (0) = k y (0) = k 1 to k x (τ ) = k y (τ ) = k 2 as a function of k 2 , and for various values of k 1 . W PC irr is the minimum work dissipation attainable with partial control, i.e. when modifying only k x and k y between the end points. W PC irr is always larger than its corresponding total control counterpart, W TC irr . The latter is attainable when controlling also k int and a during the transformation, while returning them to their initial values, k int (τ ) = k int (0) and a(τ ) = a(0). Both k 1 and k 2 are given in units of k int (0) = 1. The remaining parameters, T , a(0) and τ , are set to 1.

Interacting particles in double trap
First, we show how partial control over a system can substantially increase the amount of dissipation when compared to the optimal geodesics transformation. Consider the case of two particles on a line, Romeo and Juliet, who are constrained to be located at two different places, separated by a distance a. That is, Romeo (Juliet) is at position x (y) and subject to a confining harmonic potential of strength k x (k y ) centered at a 2 (− a 2 ). At the same time, the two particles feel a harmonic attraction of strength k int . The complete system is described by the potential (cf. Fig. 1) For two colloidal particles, such an interaction can be realized using optical tweezers [66] or an effective potential induced by feedback control [67]. Besides, it qualitatively mimics the interaction of trapped active particles studied in Ref. [68] or in a similar model [69].
(a) The commutative strategy (i) is highly dissipative, while the rotation strategy (ii) is comparable to the optimal protocol, dissipating just ∼ 20% more energy than (iii) for most values of ω a /ω b . Now imagine that an experimenter can operate a transformation of the Hamiltonian parameters with the goal to increase the strength of the local traps, but with a minimal energetic cost. That is, the boundaries of the transformation are k x (0) = k y (0) < k x (τ ) = k y (τ ) while a(0) = a(τ ) and k int (0) = k int (τ ). We want to know the minimum dissipation that an experimenter can achieve for such a transformation. In the Appendix C.1, we derive and compare the minimum dissipation protocols in the slow-driving limit for two paradigmatic cases: i) partial control in which the experimenter can only tune the values of k x and k y (while a and k int are both constant), ii) total control in which the experimenter can control k int and a as well. The comparison among the two situations is presented in Fig. 1. As expected, the dissipation under partial control, W PC irr , is always larger than that for total control, W TC irr . In particular we observe that having the possibility of controlling all the parameters of the potential (25) allows, in general, substantial savings of up to 50% of energy dissipation with respect to simply tuning the stiffnesses k x,y .

Rotating a 2-dimensional squeezed potential
As a second example, we consider the rotation of a two-dimensional Gaussian system in the xy plane. Specifically, we consider a Gaussian PDF with a non-isotropic covariance matrix of the position coordinates {x, y} which is squeezed with the major axis and the x-axis forming an angle θ in the xy plane.
Denoting the eigenvalues of Σ as ω a and ω b , it can be written in the form Our goal is to find a minimum-dissipation protocol that would overall rotate the system in the xy plane, by an angle δθ = π 4 . We thus consider protocols starting at Σ(0) = Σ 0,ωa,ω b and ending at Σ(τ ) = Σ π 4 ,ωa,ω b (cf. Fig. 2), and compare three strategies for accomplishing this task: (i) Uniform rotation trajectory: the experimenter simply rotates the system, i.e., fixes the variances ω a,b in Σ and increases θ; (ii) Commutative trajectory: the experimenter tunes ω a,b to the same value at an intermediate time, thus making the distribution isotropic. Afterwards, they re-stretch the distribution in the desired direction. Such protocol is commutative in the sense that [Σ(t), Σ(t)] = 0 ∀t (notice however [Σ(0), Σ(τ )] = 0). The intermediate point can be chosen optimally to minimize the dissipation; (iii) Optimal protocol: the experimenter is able to control the system such that it follows the BW-geodesicΣ(t) in (14) between Σ 0,ωa,ω b and Σ π 4 ωa,ω b . Details of the calculations for each of the trajectories are given in Appendix C.2 and the results are depicted in Fig. 2. We find that the commutative strategy is strongly non-optimal, and dissipates at least twice as much as the geodesics trajectory. Notice that this is not in contradiction to Observation 1, as non-commuting boundary condition induce, in general, non-commuting optimal trajectories. At the same time we see that the uniform rotation of the system Σ is close to the optimal (geodesic) trajectory in terms of dissipation, while being much simpler to implement. Finally, no timescale approximation was used in this case, but we notice that in the slow-driving limit the above strategies are equivalently translated to the stiffness matrix K = T Σ −1 (cf. Observation 2).

Discussion
In this paper, we have studied the work, heat exchange, and irreversible work dissipation of overdamped multidimensional classical systems. These may have an arbitrary number of degrees of freedom and are confined by harmonic potentials whose parameters can be partially or totally controlled. Such systems are described by multidimensional Gaussian probability distributions [53]. For uniform friction and non-moving trap centers, we have derived a general analytic expression (9) for the irreversible work (∝ entropy production). This expression is valid for any response trajectory, and allows geometric optimisation based on the Bures-Wasserstein metric for positive matrices. We also discussed straightforward generalizations of these results to non-uniform friction values and nontrivial trap center dynamics. Given that in the slow-driving limit there is a oneto-one mapping between the set of reachable states Σ and the set of reachable controls K, this further allows optimization of control protocols that incorporate experimental constraints, i.e. partial control. Finally, we described general design principles for optimal parameter protocols that minimise dissipation and illustrated them for two examples, increasing local confinement of two interacting particles and rotating a squeezed potential.
The obtained results point towards the manageable optimization of control protocols in experimental systems with many degrees of freedom [67,68,[70][71][72], and they can be directly applied to optical tweezers setups and electric circuits [66,73,74], that wish to minimise dissipation by choosing optimised control parameter protocols. Moreover, our findings can be readily applied to the case of engines and refrigerators described in the low-dissipation regime, characterized by the 1/τ scaling of dissipation [15,30,[75][76][77][78]. Further extensions include the analysis of underdamped classical systems, as well as that of quantum Gaussian systems. When applied to complex positive matrices of unit trace (that is, states in the field of quantum information), such metric represents a fundamental quantifier in problems of quantum statistical inference and metrology [55,[57][58][59]. At the same time it has its own relevance as a distance quantifier between positive real matrices or multivariate distributions, in the context of optimal transport theory [56,60]. The integrated geodesics length between two points Σ 1 and Σ 2 , if no constraints are imposed on the trace of the matrices, reads and the corresponding geodesics is where the square root is the only matrix R satisfying R 2 = Σ 1 Σ 2 and having a positive spectrum (cf. [56]).

Appendix B. General case and slow-driving solution
Here we consider the case in which also the first moment of the quadratic potential can be driven. We assume [53] that the state is Gaussian with covariance matrix Σ(t) and first moments ξ(t) (we avoid expliciting time when possible) while the potential is with ξ = z in general. The energy of the system is therefore where Σ z is the covariance matrix centered in z, that is The Langevin equation (2) becomes accordinglyẋ = −K(x − z) + √ 2T η in natural units, which is translated on the Gaussian moments aṡ where the partial derivative in time is due to the fact that Σ z depends as well on z, i.e.
The work and heat can be computed by simply taking the derivative w.r.t. the driving parameters K, z (for the work), and the dynamical parameters Σ, ξ (for the heat), i.e.
Using the same steps as in the main text (integration by parts and Eq. (B.6)), this expression translates to which can be rewritten as with the BW metrics (10) g. Notice that in general 1 2 ∆ det Σ z = 1 2 ∆ det Σ = ∆S and therefore the expression above cannot be identified as the irreversible work. At the same time, minimizing dissipation requires using finite time protocols in which the system ends in equilibrium with the thermal bath, so that no dissipation follows the end of the protocol. This is automatically satisfied in the case of slow-protocols (see below). For general transformations, it is sufficient to add a final quench of the controls, (K(τ ), z(τ )) = (T Σ −1 z (τ ), ξ(τ )). The condition ξ(τ ) = z(τ ) is sufficient to rewrite (B.10) as Appendix B.1. The slow case In the slow-driving regime a first order expansion is performed around 1 τ 0 [63]. For example in the quasistatic limit of τ → ∞ the solution for the dynamics (B.5,B.6) is clearly ξ (0) = z and Σ (B.12) (B.14) Using the above expressions, the irreversible work reads which can be rewritten in matrix form as To use Eq. (19), we rewrite the potential in the "canonical form" is the effective center of the potential. The scalar − 1 2 a (K + K int )a + 1 2 a Ka is just a global shift in energy that does not depend on the dynamics of the system and vanishes for cyclic protocols.
Confining the particles -Irreversibility parameter. We compute the irreversible work using the slow-driving approximation Eq. (19), in which the center of the distribution can be substituted by the center of the potential, and the covariance matrix can be substituted by the inverse stiffness matrix (cf. Appendix B), leading to Suppose the experimenter wants increase the strength of the local traps to increase the confinement of the two particles. The endpoint of the transformation will therefore be a(0) = a(τ ) , k int (0) = k int (τ ) , k x (0) = k y (0) < k x (τ ) = k y (τ ) . (C.7) We notice that due to the symmetry of the potential at the boundary, the eigenvectors of K + K int are always (1, 1) and (1, −1), independently of the values of k int and k x = k y . That is, [(K +K int )(0), (K +K int )(τ )] = 0 and we can therefore assume that it commutes with itself at all times (cf. Observation1). In such case, the contribution of g Σ (Σ,Σ) to W irr simplifies to (cf. Eq. (16) where ω i are the eigenvalues of Σ = T (K + K int ) −1 , which are easily computed. In particular given k x = k y ≡ k we have (C.12) The associated dissipation of the transformation can be computed from the expressions above for any slow transformation. We now consider the partial control (PC) case in which the distance a and interaction strength k int is fixed, and the experimenter can only control the local stiffnesses k x = k y ≡ k. Substituting (C.10) and (C.12) witḣ a =k int = 0, the irreversible work (C.6) then specifies tȯ For fixed boundary values of k, it can be proven using the Cauchy-Schwarz inequality that the dissipation with partial control (C.13) is lower-bounded by (C.14)