Sustaining a temperature difference

We derive an expression for the minimal rate of entropy that sustains two reservoirs at different temperatures $T_0$ and $T_\ell$. The law displays an intuitive $\ell^{-1}$ dependency on the relative distance and a characterisic $\log^2 (T_\ell/T_0)$ dependency on the boundary temperatures. First we give a back-of-envelope argument based on the Fourier Law (FL) of conduction, showing that the least-dissipation profile is exponential. Then we revisit a model of a chain of oscillators, each coupled to a heat reservoir. In the limit of large damping we reobtain the exponential and squared-log behaviors, providing a self-consistent derivation of the FL. For small damping"equipartition frustration"leads to a well-known balistic behaviour, whose incompatibility with the FL posed a long-time challenge.


Introduction
Temperature differences and gradients are a common motif in the physics of systems out of equilibrium. Primarily they serve as fixed boundary conditions for studying how energy flows within a system [1][2][3] and couples to other currents, e.g. electric or matter ones [4][5][6]. Here instead we will be interested in a dual question: what is the least amount of energy that has to be dissipated outside the system by some apparatus whose only task is to sustain a temperature difference?
More precisely, in this work we study lower bounds to the entropy production rate (EPR) σ of a conductor with respect to the profile of temperatures in the bulk, given those at the boundary. In linear systems held at temperatures T 0 and T at the extremities we find the simple expression σ * ∝ −1 log 2 (T /T 0 ). We first provide a simple heuristic argument based on the Fourier Law (FL) of conduction, and then rederive our results in a stochastic model of a linear chain of harmonic oscillators coupled to heat reservoirs, analyzed in the light of the First and Second laws of (stochastic) thermodynamics [7]. As a main new technical result we obtain explicit second-order expressions for the stationary distribution and the EPR, that may be applied in the optimization of more complex networks of interacting nodes at different temperatures (e.g. power grids [8][9][10]).
We also approach some foundational issues from a new angle. Microscopic derivations of the FL have for long been considered a challenge to the theorist [2,[11][12][13]. The problem here is to reconcile the balistic behaviour in the bulk, that is supposed to be adiabatically isolated from the environment, and the diffusive character of heat conduction. In our approach, along similar lines as in Refs. [14][15][16][17], we open the bulk to interactions with the environment. For example, if we think of a refrigerator with T 0 and T respectively the temperatures inside and outide the cold room, in the "challenge" T (x), x ∈ (0, ) would represent the temperature profile within the adiabatic walls, while in our approach it rather describes the temperature of the refrigerating liquid in the cooling coil (see Fig. 1). In the so-called overdamping limit, the FL emerges self-consistently (but not constructively), while the low-noise limit leads to the known balistic behaviour, and to a phenomenon of "frustration" whereby the system's bulk temperatures differ from their environment's. Furthermore, notice that for fixed T , letting T 0 → 0 the EPR diverges: reaching zero temperature requires ever-increasing dissipated power, providing a self-consistent formulation of the Third Law of thermodynamics -that thus is not independent of the First and Second.

Heuristics
We consider an extended system whose degrees of freedom are localized by position x ∈ X in space, and for which it makes sense to talk about a local temperature T (x), constrained by some boundary values T (∂X). The existence of a meaningful temperature is usually called local equilibrium or local detailed balance. We point out that it is not the system itself to be at equilibrium, but rather a continuous set of thermometers whose local degrees of freedom are labeled by x, and that interact via the system.
To a temperature gradient is associated a thermal force F = ∇T −1 , which generates a heat current J. If we provisionally take the FL for granted, the stationary heat current is given by J = −κ ∇T , with κ the heat conductivity. The macroscopic stationary EPR is given by the scalar product We rewrote the expression in a couple of equivalent useful ways. where n is the unit vector Figure 1: Illustration of constructive vs. self-consistent approaches to heat conduction in temperature gradients. In the former heat flows from the hot reservoir to the cold one through the oscillators. It former draws its motivation in the derivation of the physics of open systems from that of isolated systems (micro-to-meso/scopic); its dissipative observables are the temperatures of the bulk oscillators. In the latter heat flows to and from each reservoir. This approach focuses on the temperature of the baths, and aims to connect to broad-scale energy consumption (meso-to-macro/scopic).
orthogonal to the boundary surface element. Notice that the second expression in Eq. (A.2) manifests an invariance under the transformation T → 1/T . We search for extremals of σ by taking the functional derivative δ/δT (y). After some standard integration by parts and application of Dirac deltas (see details in the Appendix ) we obtain the Laplace equation ∆ log T * = 0, where the asterisk stands for the extremal. Plugging into the EPR we find which is easily checked to be a minimum. Notice the squared-log dependency on boundary temperatures.
We now reduce to one dimension by assuming that the temperature only varies in one extended direction x ∈ [0, ], while at fixed x it is uniform in a perpendicular area of fixed size α. Subjecting the equation to the boundary constraints T (0) = T 0 and T ( ) = T , we obtain as the unique minimum profile The least EPR is then given by Notice that the shorter the distance, the steeper the gradient, the higher the EPR.

Model
Let us now re-derive the above results in a well-known microscopic model, employed e.g. in Refs. [16,17] in attempts at derivations of the FL by self-consistent reservoirs, and in Ref. [15] to discuss energy equipartition of normal modes. We consider a homogeneous chain of n harmonic oscillators of unit mass placed at regularly spaced positions x k = (k − 1) /(n − 1), with k = 1, . . . , n, between boundary x = 0 and x = .
We also set Boltzmann's constant to unity k B = 1. We denote q k the amplitude of oscillation of the k-th oscillator and p k its momentum, and collect z = (q k , p k ) n k=1 . These amplitues are the effective degrees of freedom that describe the local interaction of the system with the baths. The total energy is where q 0 = q n+1 = p n+1 = 0 stand for the non-oscillating endpoints where the chain is anchored and ω is the angular frequency, which for the moment we also set to unity ω = 1 to resume it later in the discussion. Each harmonic oscillator is in contact with a heat bath at temperature T k , which is a source of stochastic white noise (the reduced effect of the bath degrees of freedom). The dynamics is described by the Langevin equatioṅ where ζ k (t) are the formal time derivatives of independent Brownian motions, and γ is the damping coefficient. Letting D be the diagonal positive-definite diffusion matrix D = diag {T k } n k=1 , and A the symmetric tridiagonal matrix and further defining the 2n × 2n matrices with I the identity, the equations of motion take the form of a multidimensional Ornstein-Uhlenbeck (OU) processż where ζ(t) are 2n independent delta-correlated white noises. Letting ρ(z, t)dz be the probability of z(t) being in a neighborhood of z, its density satisfies the Kramers diffusion equation In the second identity in Eq. (10), { · , · } denotes the Poisson bracket and the dissipative current is given by [18,20] j : showing that the underdamped dynamics clearly separates into a balistic term and a diffusive one.
The stationary probability ρ ∞ (z) := lim t→∞ ρ(z, t) is found by assuming as ansatz a centered multinormal distribution is the stationary correlation matrix, with C qq and C pp symmetric, and denoting transposition.

Stationary distribution
For the sake of generalization we momentarily allow for a non-symmetric A. Plugging the stationary distribution into the continuity equation ∂ p · j ∞ = 0 one finds that the covariance matrix is uniquely determined by the (first-order) Lyapunov equation [21] CM + MC = 2D. (13) In view of Eq. (12), one finds that C pq = −C pq is antisymmetric, and furthermore: The first defines C qq in terms of C pq , the second C qq in terms of C pp , and the third C pq in terms of C pp . Combining these formulas together we find the second order Lyapunov equation for the displacements' covariance C := C qq In the overdamping limit γ → ∞ this expression reduces to the well-known AC + CA = 2D [21,22]; furthermore we have C pp → D (local equipartition) and C pq ∼ 1 γ (AC − D). Defining the antisymmetric "curvature" tensor R := D −1 A − A D −1 = 0 [19], the condition of (global) detailed balance states that R vanishes, if and only if C pq vanishes [23]. Proof. If detailed balance holds then C = DA −1 solves the second order Lyapunov equation and C pq = 0. Vice versa, if C pq = 0 then the first term in Eq. (15) vanishes, one obtains D = AC = CA and therefore R = 0.
In our case detailed balance is achieved when D ∝ I, i.e. for all equal temperatures.

Entropy production rate
Before turning to the solution of the second-order Lyapunov equation, we introduce thermodynamic quantities. We now take the time derivative of the mean energy where the average is over realizations of the Brownian motions, and we used the Itô Lemma and the martingale property. The First Law d H + kd Q k = 0 suggests to define the heat flow from the k-th reservoir asd Notice that it vanishes when equipartition holds. The Clausius formula defines the EPR, where S := − dz ρ log ρ is the Gibbs-Shannon entropy. The second identity, providing the Second Law σ ≥ 0, is proven in the SM. At the stationary state, plugging the definition of the heat flux we find σ ∞ = γ tr D −1 (C pp − D) = −tr (C pq R)/2, and finally where we employed the Lyapunov equation for C pq and its antisymmetry. Both this latter expression for the EPR and the Lyapunov equation are clearly independent of the vector basis chosen to represent matrices. We will now work with normal modes, i.e. orthonormal eigenvectors a α of A. Then A = α λ α a α ⊗a α , where λ α are the real eigenvalues and ⊗ denotes the outer product. In this basis the diffusion and covariance matrices have entries respectively D αβ = a α · Da β and C αβ = a α · Ca β . Finally we can express Eq. (15) in this basis to find where we finally resumed the angular frequency ω simply by rescaling all eigevalues λ α → ω 2 λ α . In our problem, the tridiagonal matrix A has real eigenvalues λ α = 2 + 2 cos(απ/(n + 1)) and orthonormal eigenvectors'entries a k α = 2/(n + 1) sin(αkπ/(n + 1)), for α, k = 1, . . . , n, ranging from slower to faster modes (see SM and Ref. [24]). We can thus compute D αβ ; notice that by orthonormality, (D −1 ) αβ can be obtained from D αβ by the duality transformation In the overdamping limit first neighbours are strongly favoured, while in the underdamped limit spatial proximity in the bulk does not play a role, and only the endpoint oscillators are slightly favoured.
where in the first expression Matrix ∆ is symmetric, its diagonal entries are positive, and by orthonormality k a k α a k β = δ α,β one finds that k ∆ kh = 0. Therefore ∆ is a proper discretized Laplacian (see also Refs. [25,26] in the quantum case), and Eq. (23) is a discretized equivalent of the second expression in Eq. (A.2), whose first expression is instead reminiscent of the second in Eq. (23) with F k,h := 1/T h − 1/T k as the thermodynamic force due to the competition of the k-th and the h-th reservoir and J k,h := γ∆ kh (T k − T h ). Then the main difference with the heuristic case is that, despite the fact that the oscillators only interact with nearest neighbors, they create all-to-all nonequilibrium currents (see also Ref. [27]). However, Fig. 2 shows that for large γ first-neighbour contributions indeed dominate.
Finally, we minimize the EPR with respect to the temperatures of the bulk oscillators, subject to constrained values of the temperatures of the first and last oscillators. For n = 3 the free oscillator's temperature is easily found to be T * 2 = √ T 1 T 3 (see SM), yielding an analytical expression of the minimum EPR whose most interesting feature is that as a function of γ it vanishes for γ → 0, ∞ and has a maximum in between, whose physical significance is still an open question.
For n > 1 minimization could only be achieved computationally. Fig. 3 shows that in the overdamping limit in the oscillator model the optimal temperature profile is consistent with that from the FL, while in the underdamped limit the bulk oscillators's optimal temperatures flatten, reproducing the behaviour observed in Ref. [14]. With crosses we plotted the mean squared momentum of the oscillators, calculated as While in the overdamping limit equipartition is reached, in the low-damping limit a phenomenon of temperature frustration occurs, whereby the bulk's internal temperatures are slightly off the environment's. Finally, in the overdamping limit the overlapping of all curves in Fig. 4 proves the squaredlog behaviour and the fact that, for given temperature difference well beyond the linear regime, the minimal EPR scales like n −1 , reproducing the −1 dependence in the FL.

Possible developments
To conclude, in this paper we collected theoretical evidence that, for systems open to the interaction with the environment, on the assumption that thermometers can be defined locally, the minimal entropic cost of maintaining a temperature gradient at the two extremities of a linearly extended body scales with the squared logarithm of the temperature ratio, and inversely with the distance. A main novelty of this work is the idea of viewing a temperature gradient "from the outside" instead of "from the inside". This may be of practical interest in the assessment of the industrial scaling of the energetic demand of technologies running at very low temperatures or in the optimization of power grids. Our law may provide an indirect testing ground for the thermodynamics of open systems based on Markov processes, whose experimental verification is still intertwined with the identification of the kind systems to which it applies. Finally, notice that experimental verification of the squared-log behaviour would entail a form of "minimum entropy production" principle [28]. Figure 4: In the overamped case γ = 100, for chains of n = 10, 20, 30 oscillators, n times the minimum EPR σ * as a function of the temperature of the right-most oscillator ranging from 1 to 100, for given T 0 = 1. The three data sets overlap and are fitted by a squared-log curve as in Eq. (4).

A Minimal EPR
We want to explicitly show that the temperature profile T * that minimizes the macroscopic stationary EPR with our notion of local temperature follows the Laplace equation By considering (κ = 1) the entropy σ = σ[T ] is a functional of the temperature field T (x). By taking the functional derivative of σ respect to T where we applied Dirac's deltas coming from the functional derivative. The first term is and the second one gives which is true when ∆ log T * = 0 holds. The minimal EPR is obtained by integrating (A.2) by parts: The second contribution vanishes because of Eq. (A.1), then When restricting our study to one dimension the Laplace equation gives the minimal exponential profile since ∂ x T ((1/T )∂ x T ) = 0 when ∂ x T = const T which is a first-order differential equation whose solution is a 2-parameters exponential. The parameters are obtained with the two conditions T (0) = T 0 and T ( ) = T .

B Underdamped EPR
Consider the second expression in Eq. (14). We have where in the second term we integrated by parts and recovered the definition of the current. Now, employing the continuity equation Eq. (8) we obtain for the second term where S(t) = −k B zρ(z, t) log ρ(z, t) is the Gibbs-Shannon entropy, and we used the Liouville theorem and probability conservation. Also, notice that dz p · D −1 p ρ = tr [D −1 C pp (t)] leading to Finally this last term can be integrated by parts yielding −nγk B , thus recovering the first expression in Eq. (14). where we used that sin(x) = 1 2i (e ix − e −ix ) and that the partial sum of the geometric series is given by m k=0 x k = (1 − x m+1 )/(1 − x). The result is obtained by noticing that our summations runs from k = 1 and that e 2izπ = cos(2zπ) = 1, ∀z ∈ Z.
With similar arguments we can prove that {a α } is also orthogonal: in fact since for α = β we fall in the previous case to show normality of a α , and for α = β we can say that when α + β is even (odd) also α − β is even (odd), and in both case the expression above vanishes because cos(mπ) = 1 for even m and cos(mπ) = −1 for odd m.

D Exact solution for 3 oscillators
Let us consider a system of 3 interacting harmonic oscillators where the endpoints (oscillators 1 and 3) are coupled with thermal reservoirs with given temperatures T 1 = T c and T 3 = T h .