
==== Front
J Stat Phys
J Stat Phys
Journal of Statistical Physics
0022-4715
1572-9613
Springer US New York

3320
10.1007/s10955-024-03320-w
Article
Optimal Control of Underdamped Systems: An Analytic Approach
http://orcid.org/0009-0001-3690-1682
Sanders Julia 1
http://orcid.org/0000-0003-3559-4032
Baldovin Marco marco.baldovin@cnr.it

2
https://orcid.org/0000-0003-0241-6619
Muratore-Ginanneschi Paolo 1
1 https://ror.org/040af2s02 grid.7737.4 0000 0004 0410 2071 Department of Mathematics and Statistics, University of Helsinki, 00014 Helsinki, Finland
2 https://ror.org/05rcgef49 grid.472642.1 Institute for Complex Systems, CNR, 00185 Rome, Italy
Communicated by Udo Seifert.

17 9 2024
17 9 2024
2024
191 9 11721 3 2024
4 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Optimal control theory deals with finding protocols to steer a system between assigned initial and final states, such that a trajectory-dependent cost function is minimized. The application of optimal control to stochastic systems is an open and challenging research frontier, with a spectrum of applications ranging from stochastic thermodynamics to biophysics and data science. Among these, the design of nanoscale electronic components motivates the study of underdamped dynamics, leading to practical and conceptual difficulties. In this work, we develop analytic techniques to determine protocols steering finite time transitions at a minimum thermodynamic cost for stochastic underdamped dynamics. As cost functions, we consider two paradigmatic thermodynamic indicators. The first is the Kullback–Leibler divergence between the probability measure of the controlled process and that of a reference process. The corresponding optimization problem is the underdamped version of the Schrödinger diffusion problem that has been widely studied in the overdamped regime. The second is the mean entropy production during the transition, corresponding to the second law of modern stochastic thermodynamics. For transitions between Gaussian states, we show that optimal protocols satisfy a Lyapunov equation, a central tool in stability analysis of dynamical systems. For transitions between states described by general Maxwell-Boltzmann distributions, we introduce an infinite-dimensional version of the Poincaré-Lindstedt multiscale perturbation theory around the overdamped limit. This technique fundamentally improves the standard multiscale expansion. Indeed, it enables the explicit computation of momentum cumulants, whose variation in time is a distinctive trait of underdamped dynamics and is directly accessible to experimental observation. Our results allow us to numerically study cost asymmetries in expansion and compression processes and make predictions for inertial corrections to optimal protocols in the Landauer erasure problem at the nanoscale.

Keywords

Optimal control
Multiscale analysis
Underdamped dynamics
http://dx.doi.org/10.13039/501100000781 European Research Council 785932 Baldovin Marco Centre of Excellence in Randomness and Structures of the Academy of Finlandissue-copyright-statement© Springer Science+Business Media, LLC, part of Springer Nature 2024
==== Body
pmcIntroduction

In his remarkable paper [1] (English translation in [2]), Schrödinger addresses the problem of statistical reversibility of a physical system in contact with an environment. In doing so, he puts forward the idea of using entropic indicators to quantify deviations from thermodynamic equilibrium and, therefore, dissipation. Schrödinger identifies what is now commonly known as the Kullback–Leibler divergence or relative entropy [3] as a quantifier between the joint probability distribution of the system’s end states and those of a free diffusion.

In the last decades of the 20th century, Schrödinger’s trailblazing idea was reformulated into the language of stochastic optimal control [4–6], refining Schrödinger’s original “static bridge problem” [1] into a “dynamic Schrödinger bridge", where the relative entropy is computed between probability measures over the systems’ pathspace [7, 8]. Schrödinger bridges have active research interest because they allow computational optimal transport methods to be applied to dynamical models. This enables efficient computation in fields such as neuroscience [9, 10]; data science and machine learning [11]; and generative modeling, sampling, and dataset imputation [12, 13].

Technological advances in the last two decades have paved the way for the observation and manufacturing of nanomachines. At nanoscale, random fluctuations of thermal and topological origin may swamp out any mechanical behavior [14]. A fundamental question is, therefore, how natural or artificial nanosystems can efficiently harness randomness in order to generate controlled motion or perform thermodynamic work on larger scales. Schrödinger bridges find an optimal control protocol to rectify a system obeying stochastic dynamics, thus making it possible to devise systematic methods characterizing the efficiency of nanomachines [15].

In addition, the discovery of fluctuation relations (see Chapter 4 of [16] for a thorough conceptual and historical account) introduces a substantial development with respect to [1]. For Markov stochastic processes, fluctuation relations stem from considering the relative entropy of probability measures connected by a time reversal [17–19]: the corresponding generalized Schrödinger bridges minimize the thermodynamic cost associated with the transition. This means that we can consider thermodynamic cost functionals other than the Kullback–Leibler divergence [20–27]. Because the Kulllback-Leibler divergence is non-negative, the minimum mean entropy production in finite time transitions between assigned probability distributions is strictly larger than zero. Remarkably, the overdamped dynamics minimizer [28, 29] turns out to be the solution of a system of Monge–Ampère–Kantorovich optimal mass transport equations [30]. Two results of [28, 29] stand out. First, minimizers can be determined by efficient numerical algorithms even in the multi-dimensional case [31]. Second, the minimum entropy production is proportional to the squared Wasserstein distance between the probability distributions of the end states divided by the duration of the control horizon. This relation between mean entropy production and squared Wasserstein distance continues to hold as scaling limits for Markov jump processes [32] and underdamped dynamics [33, 34], see also [35, 36].

A detailed description of optimal protocols in the underdamped regime is urgent for several reasons. Optically levitated nanoparticles have become a common tool to study transitions in stochastic thermodynamics. Stable confinement and manipulation of nanoparticles within optical traps requires an account of the momentum dynamics. For instance, particle-environment energy exchanges during isoentropic (isochoric) transition within Brownian Carnot (Stirling) engines occurs through the momentum degrees of freedom [37, 38]. Understanding how to simultaneously control particles’ position and momentum is required to devise robust shortcuts to equilibration protocols [39–43].

A further motivation comes from the design of nanoscale electronic components [44–46]. Increasing the efficiency of such operations toward the bound prescribed by Landauer is a non-trivial task, with potentially relevant consequences for the design of information and computation technology [47]. The presence of inertia has been shown to lower the energetic cost needed to perform logic operations on bits [44, 48, 49]. This has sparked interest in the control of underdamped stochastic systems, with particular emphasis on the non-linear case, needed for the description of information bits [12, 50–52]. Ad-hoc experimental solutions have been found to realize controlled protocols for stochastic dynamics with inertia, confirming that inertial effects allow for fast and precise bit operations [49, 53–55].

With these motivations in mind, we introduce a systematic analytical derivation of optimal protocols for the underdamped dynamics. We consider two paradigmatic cases of running costs:The underdamped version of the Schrödinger dynamic bridge problem, referred to as KL. The cost functional in this case is the Kullback–Leibler divergence between the probability measure of the controlled process and that of an assigned reference process. In [1], Schrödinger motivates this cost functional as a quantifier of the likelihood of non-equilibrium fluctuations, i.e., a large deviation functional in the current literature. Since then, Schrödinger dynamic bridge problems have emerged as relevant efficiency measure for diffusion-mediated transport processes with applications ranging from cybernetics [10] to molecular-scale engines [15]. More recently, it has been realized that when the reference process is a diffusion subject to inertia in the absence of a confining potential, Schrödinger dynamic bridge problems provide a viscous regularization of optimal mass transport [56–58]. This regularization has applications in machine learning [11]. Finally, when the reference process describes motion in a confining potential equal to that in the Maxwell-Boltzmann distribution of the final state, we obtain a model of an optimally controlled shortcut to adiabaticity.

The minimization of the mean entropy production, referred to as EP. This is the cost functional characterizing the second law of thermodynamics and Landauer’s principle [47]. We study this problem in the most general formulation compatible with detailed balance, and which can be self-consistently derived from Hamiltonian mechanics of a system coupled to an infinite bath described by harmonic oscillators. In this generalized formulation, the underdamped entropy production explicitly depends upon the control via a non-dimensional parameter g. Physically, g describes the intensity of momentum coupling between the system and bath. Interactions of these type have been recently observed in Josephson junctions [59, 60].

Besides the cost functional, the specification of an optimal control problem requires a definition of the functional space over which to carry out the optimization. This functional space is called the class of admissible controls. Our focus is on the class of admissible controls described by functions that are sufficiently regular to be differentiable in space and continuous time. Physically, this corresponds to the requirement that the control be slow with respect to the fastest time scale in the problem set by the Wiener process modeling the interaction to the bath. Under these hypotheses, we derive the stationary equations for the cost functionals by taking variations over the class of admissible controls specified by confining mechanical smooth potentials. In such a case [34], extremals of the cost solve a set of integro-differential equations, with features reminiscent of the Vlasov-Poisson-Fokker-Planck problem [61].

We obtain the following main results: I We show that the cumulants of the probability measure describing transitions between Gaussian states are amenable to the solution of a Lyapunov system of equations [62] in any number of dimensions. This immediately yields a body of rigorous results concerning existence, uniqueness and, when applicable, positivity of solutions (Sect. 4).

II For transitions between states described by Maxwell-Boltzmann distributions in phase space, we introduce an infinite dimensional extension of Poincaré–Lindstedt multiscale perturbation theory [63] around the overdamped limit. This method allows us to treat all cumulants of the system probability measure on the same footing in the renormalization group fashion [64]. We hence obtain explicit predictions for the behavior in time of all phase space cumulants within second order accuracy. The method builds on ideas introduced in [65, 66] for dissipative and [67] for conservative dynamics. Although we restrict our analysis to a two-dimensional phase space, the analysis of the Gaussian case shows that extension to higher dimensional phase spaces is possible, albeit cumbersome (Sect. 6).

III In the case of mean entropy production by an underdamped dynamics with purely mechanical coupling, our results support tightness of the lower bound provided by the overdamped dynamics [68, 69]. For more general couplings, both the mean entropy production and the cost of the dynamic Schrödinger bridge receive strictly positive corrections in the presence of inertia (Sect. 6.1.4).

IV The cost of expansion is higher than that of compression when the initial states are thermodynamically equidistant (Sect. 8.1.1). This result is a manifestation of intrinsic asymmetries in thermal kinematics, recently pointed out in [70, 71].

The structure of the paper is as follows. In Sect. 2, we introduce the model of underdamped dynamics of a nanosystem weakly coupled to an environment by both mechanical and momentum dissipation interaction. When the intensity of the momentum coupling g vanishes, the model recovers the most widely applied underdamped dynamics. Next, we consider two thermodynamic cost functionals, KL and EP, and motivate their broad interest for applications to physics and other applied sciences. Our goal is to minimize these functionals towards the mechanical potential Ut governing the underdamped dynamics conditioned on the system’s initial and final probability distributions. For this reason, we present a brief overview of the mathematical results leading to known bounds for the cost functionals KL and EP in the second half of the section.

In Sect. 3, we introduce the Pontryagin-Bismut functional and derive its stationary equations. The Pontryagin-Bismut functional provides a description of optimal control dual to Bellman’s principle.

Section 4 focuses on the Gaussian case in a phase space of arbitrary dimension, and we derive our first main result here.

In Sect. 5, we set the stage for multiscale perturbation theory presented in Sect. 6. As usual, the idea is to use slow scales to cancel secular terms. Our main goal is to obtain a detailed analytical description of experimentally measurable indicators. We therefore summarize the logic of the derivation and the results before proofs. Readers only interested in our results may thus skip the second part of Sect. 6.

In Sect. 7, we briefly return to the Gaussian case and provide the analytic expression of the solution of the cell problem of the multiscale expansion [72]. The solution of the cell problem allows us to determine all cumulants within second order accuracy in the overdamped expansion.

Section 8 applies the results with some numerical computations. We have emphasized the Gaussian case for two reasons. Firstly, methods for accurate numeric integration of the exact optimal control equations are immediately available, meaning we can compare the perturbative approach with exact numeric predictions in the case of Gaussian boundary conditions. Secondly, transitions between Gaussian states are well adapted to model Brownian engines [70, 71, 73–75]. We therefore also study the cost of optimal protocols driving isothermal expansions and compressions of a system to an equilibrium state, which are modelled by a dynamic Schrödinger bridge. Additionally, we solve the cell problem in the case of Landauer’s erasure problem numerically and thus find inertial corrections to the erasure protocol, as well as predictions for the system’s probability measure cumulants.

The final section is devoted to conclusions and outlook. We defer further supplementary material to the Appendices.

Underdamped Control Model

We consider the dynamics of a nanosystem with mass m, whose position qt and momentum pt obey the Langevin–Kramers stochastic differential equations in R2d1a dqt=ptm-gτm(∂Ut)(qt)dt+2gτmβdwt(1)

1b dpt=-ptτ+(∂Ut)(qt)dt+2mτβdwt(2).

In Eqs. (1), wt(1) and wt(2) denote two d-dimensional independent Wiener processes. The Stokes time τ is a constant parameter specifying the characteristic time scale of dissipation.

In (1a), a non-dimensional constant g couples the mechanical force ∂U and the fluctuating environment modeled by the Wiener process wt(1) to the nanosystem position dynamics. For any g≥0, Eq. (1) guarantees convergence towards a Maxwell-Boltzmann equilibrium whenever the potential U is time independent, confining and sufficiently regular. Setting g to zero recovers the standard Langevin–Klein–Kramers model [76].

We emphasize that the dynamics described by (1) are consistent with the general analysis [77] of the conditions guaranteeing the self-consistency of the harmonic environment hypothesis. In fact, (1) can be obtained from a microscopic Hamiltonian dynamics, in which the system interacts with a “bipartite harmonic” environment [78]. By bipartite harmonic environment, we mean an environment modeled by two kinds of oscillators: one type interacts with the system via the commonly assumed position-coupling [79], and the other via a linear momentum coupling [59, 60]. Linear momentum coupling models momentum dissipation observed e.g. in a single Josephson junction interacting with the blackbody electromagnetic field.

As for the force in (1), we only assume that it is the negative gradient of a confining and sufficiently regular mechanical potential, i.e. a potential depending only on the system position. We suppose that potentials of this type give rise to an open set of controls. Within this set, the controls ensure that at every instant of time t in a given time horizon [tι,tf] the probability density of the systemPrx≤qtpt<x+d2dx=ft(x)d2dx

is well defined and satisfies the Fokker–Planck equation.

At an initial time t=tι, we posit that the state of the nanosystem is statistically described by an assigned Maxwell-Boltzmann distribution at inverse temperature β:2 ftι(q,p)=Zι-1exp-β‖p‖22m-βUι(q)

Furthermore, we require that at the end of the control horizon t=tf the probability density of the system satisfies the boundary condition3 ftf(q,p)=Zf-1exp-β‖p‖22m-βUf(q).

These assumptions on the probability distributions of the system end states are not necessary for the considerations that follow. They have the merit, however, to be both physically admissible and to lead to simplifications in the multiscale analysis of Sect. 6.Fig. 1 Stylized representation of a Schrödinger bridge modeling Landauer’s erasure of one bit of memory at minimum dissipation

The set of confining potentials Ut that give rise to phase space diffusions with probability marginals (2), (3) define the class of admissible controls of (1).

Our aim is to determine the optimal mechanical potentials Ut among the admissible ones that minimizes the thermodynamic cost functionals defined below conditioned on the initial and final probability distributions (2) and (3).

Thermodynamic Cost Functionals

We focus our attention on two physically relevant cases, hereafter referred to as KL and EP.

KL: Underdamped dynamic Schrödinger bridge [1]. The thermodynamic cost functional to minimize is the Kullback–Leibler divergence of the measure P=Pιf generated by (1) subject to (2) and (3), from the measure Q=Qι generated by (1) when the mechanical force is ∂U⋆ and only the initial density (2) is assigned. The cost functional reads (see Appendix A)4 K(P‖Q):=EPlndPdQ=βτ(1+g)4mEP∫tιfdt‖(∂Ut)(qt)-(∂U⋆)(qt)‖2.

The notation EP emphasizes that the expectation value over the diffusion path is with respect to the measure P, and dP/dQ denotes the Radon-Nikodym derivative between P and Q.

In the mathematics literature, the minimization of (4) atU⋆=0

is referred to as entropic interpolation [7] or entropic transportation cost [80]. This terminology is due to the discovery [56] that the minimization of the overdamped counterpart of (4) yields a viscous regularization of the Monge–Ampère–Kantorovich optimal transport problem (see [30]). Finally, [15] supports the use of the cost of a Schrödinger bridge as a natural efficiency measure for nano-engines in highly fluctuating environments; see [10, 81] for a wider class of applications. In Sect. 8.1.1 we show how the optimization of (4) provides a plausible model of shortcut to equilibration.

EP: Mean entropy production. In stochastic thermodynamics, the average entropy production is identified with the Kullback–Leibler divergence of the forward measure P from a measure PR obtained by a combined time-reversal and path-reversal operation (see e.g. [18, 19, 29] and Appendix A for further details):5 E=EPlndPdPR=EPlnptι(qtι,ptι)ptf(qtf,ptf)+EP∫tιfdtβ‖pt‖2mτ-dτ+βgτmEP∫tιfdt(∂U)(qt)2-(∂2U)(qt)β.

Some observations are in order. To start with, the identification of (5) as mean entropy production of during a thermodynamic transition is legitimate in consequence of its relation with the heat release by the system evolving according to (1). This is a consequence of the general theory expounded in [19]. In Appendix B for reader convenience, we reproduce the calculation that justifies the identification.

For any g, the entropy production vanishes for a system in a Maxwell-Boltzmann equilibrium. For a bridge process, equilibrium means that the boundary conditions (2), (3) are specified by the same Maxwell-Boltzmann distribution. The corresponding optimal control problem becomes trivial. In any non-trivial case, the Gibbs-Shannon entropy difference appearing in the first row of (5) does not play a role in the optimization as it is fully specified by the boundary conditions.

Finally, the entropy production is non-coercive, i.e it is not a convex functional of the control at g equal zero. As a practical consequence, none of the infinitely many time-dependent protocols that connect the selected end states can be said to minimize the entropy production. Precise treatments of the optimal control problem in such a case are possible either by regularizing the problem [34], or in special cases [75], by considering non-purely mechanical controls [33]. Studying, as we propose here, the mean entropy production at finite g has the advantage of making the cost functional coercive with respect to the mechanical force.

At this point, it is worth commenting on our working hypotheses. The cost functionals in both case KL and EP are readily convex in the mechanical potential. We surmise the existence of an open set of admissible potentials that allows us to look for a minimum in the form of a regular extremal of a variational problem [82]. To justify this assumption we recall that Hörmander’s theorem (see e.g [83]) ensures that any potential Ut (1) that is sufficiently regular, bounded from below, and growing sufficiently fast at infinity results in a smooth density.

Bounds of the Thermodynamic Cost Functionals

In practice, the cost functionals (4) and (5) are the limit of Riemann sums on ratios of transition probability densities evaluated over increasingly small time increments. This construction is recalled in Appendix A. The construction immediately implies that (4) is bounded from below by the Kullback–Leibler divergence of the joint probability distribution of the system state at the end-times of the control horizon.

The measure theoretic analysis in Sect. 3 of [6] permits drawing more precise qualitative conclusions without making direct reference to the details of the dynamics. To summarize them, let us denote by S the state space of dimension dS, where the stochastic process {xt,t∈[tι,tf]} with probability measure P takes values. We also denote by Pyx (Qyx) the probability measures subject to the bridge conditionsxtι=y&xtf=x.

Under technical hypotheses guaranteeing that the optimization problem is well-posed, the main takeaways of [6] are the following. First, the Kullback–Leibler divergence is always amenable to the decomposition [4]6 K(P‖Q)=K(℘‖℘⋆)+∫S2ddSxddSy℘(x,y)K(Pyx‖Qyx).

The first addend is the quantity originally considered by Schrödinger in [1], namely the static Kullback-Leibler divergence7 K(℘‖℘⋆)=∫S2ddSxddSy℘(x,y)ln℘(x,y)℘⋆(x,y)

of the joint probability density ℘ of xtι and xtf from the two point probability℘⋆(x,y)=Ttf,tι(Q)(x∣y)ftι(y),

which is uniquely defined by the transition probability density Ttf,tι(Q)(·∣·) of the reference process and the probability distribution of xtι.

Both addends in (6) are positive. Furthermore, the static divergence (7) vanishes if and only if℘(x,y)=℘⋆(x,y).

Many possible P are compatible with the same ℘. Once ℘ is fixed, K(P‖Q) attains an infimum, in fact a minimum, for the P that makes the second term of (6) vanish. A necessary condition [4, 6] enforcing this requirement is that for any tι≤s≤t≤tf8 P(xt=x∣xs=y,xtι=z)=ht,tι(x,z)Tt,s(Q)(x∣y)hs,tι(y,z)

where the function h is defined byht,tι(x,y)=∫SddSzhtf,tι(z,y)Ttf,t(Q)(x∣z).

Once (8) holds true, the control of the abstract optimization problem is htf,tι. Correspondingly, (6) reduces to (7), which in turn we can couch into the formKopt(P‖Q)=Kopt(℘‖℘⋆)=∫S2ddSxddSyht,tι(x,y)Ttf,tι(Q)(x∣y)ftι(y)lnhtf,tι(x,y).

The general form (8) of the necessary condition for the reduction to a static problem does not require the optimal process to enjoy the Markov property; the transition probability (8) may carry memory of the value taken by xtι. The results of [6] ensure that (8), under further regularity assumptions, reduces for any s≤t∈tι,tf to a Markov transition probability densityTt,s(Pιf)(x∣y)=ht(x)Tt,s(Q)(x∣y)hs(y).

From the physics point of view, the assumptions leading to Markov transition probability densities immediately include an overdamped dynamics [5] or an underdamped dynamics driven by force field depending on both the position and momentum of the system [68, 84], and thus distinct from (1).

Our discussion so far refers to case KL. The connection to case EP stems from the Talagrand-Otto-Villani inequalities [85, 86]. These inequalities show that the static Kullback–Leibler divergence between probability densities is bounded from below by the squared Wasserstein distance between the densities multiplied by a proportionality factor. For the overdamped dynamics considered in [5], Mikami [56] (see also [57, 80, 87]) later proved that the bound becomes tight in a suitable scaling limit and the proportionality factor reduces to the inverse of the duration of the control horizon. More explicitly, the entropic transport cost (U⋆=0) multiplied by the viscosity becomes equal to the cost of a Monge–Ampère–Kantorovich optimal mass transport problem [30] in the limit of vanishing viscosity.

The connection to problem EP consists in the proof [28, 29] and [58] that the minimization of the mean entropy production by bridge processes obeying the overdamped dynamics can be exactly mapped into a Monge–Ampère–Kantorovich optimal mass transport. The reason is that the optimal control problem admits an equivalent reformulation, in which the current velocity of the admissible processes [88] play the role of control instead of the drift.

In the underdamped case, the presence of inertial effects complicates the picture. The mean entropy production cannot be written as the square of the current velocity. This prevents a direct application of the Benamou–Brenier inequality [89] (see also Appendix C). The Benamou–Brenier inequality allows one to couch the minimum mean overdamped entropy production into the squared Wasserstein distance between the densities at the end of the control horizon. It is, however, possible to show [33, 34] that the underdamped mean entropy production admits its overdamped counterpart as a lower bound. In particular, for (5), the following inequality9 Etιf≥mβ(1+g)τEPqtf-qtι2tf-tι

holds true. Bounds of the type (9) for the mean entropy production appeared in [33, 34] and later in [90]. The proof of (9) presented in [69] is motivated by [91]. For reader convenience, we reproduce the proof in Appendix C.

The above considerations suggest that for both the over- and underdamped dynamics (1), the inequality10 K(P‖Q)≥CTOVmβ(1+g)τEPqtf-qtι2tf-tι

should also hold true with CTOV a positive constant in agreement with the Talagrand-Otto-Villani theory [85, 86]. We refer to [68] for a mathematical proof of the bound for the underdamped dynamics, and for the overdamped dynamics [92] (see also [58, 80, 93]), including an explicit prediction of the constant CTOV.

Optimal Control Formulation

Optimal control problems can be turned into variational problems by coupling the dynamics to the cost functional through Lagrange multipliers. In hydrodynamics, such an approach is referred to as the adjoint equation method, and has a long history going back to [94, 95]. By themselves alone, solutions of the variational equations only provide a necessary condition for the existence of regular extremals of an optimal control problem: extremals continuously satisfying (partial) differential equations. Optimality follows from convexity of the cost, as in our case, or from the study of the second variation. A mathematical formulation of the adjoint method in the context of stochastic optimal control is due to Bismut see e.g. [96]. Bismut theory may be regarded as extending the Pontryagin principle of deterministic optimal control [97] see also [82, 98]. Based on these considerations, we construct from KL and EP variational functionals by imposing that the mean forward derivative [88] of a Lagrange multiplier is an exact differential when the density ft(x) with respect to the Lebesgue measure satisfies a Fokker–Planck equation associated with probability preserving boundary conditions. As in [33, 34] we refer to these functionals as Pontryagin-Bismut. In such a setup, we look for variational extremals of11 A[f,U,V]=∫R2dd2dx(Vtι(x)fι(x)-Vtf(x)ff(x))+∫tιfdt∫R2dd2dxft(x)(Ct(Ut)(x)+(∂t+Lx)Vt(x)).

Here we collectively denote phase space coordinates asx=(q,p)

and define the running cost functional asCt(Ut)(x)=βτ(1+g)4m‖(∂Ut)(q)-(∂U⋆)(q)‖2[KL]β‖p‖2mτ-dτ+βgτm(∂Ut)(q)2-(∂2Ut)(q)β.[EP]

In writing (11) we conceptualize the fields f, V, and U as unknown variational fields. The existence of the functional requires integrability with respect to ft which we assume to be a probability density taking the values fι and ff at the start and end of the control horizon respectively, fixed by (2), (3). The field V becomes the value function of Bellman’s formulation of optimal control theory [97]. In (11), it plays the role of a Lagrange multiplier enforcing the dynamics.

Accordingly, we denote by Lx the differential generator of the dynamics determined by (1):12 Lx=p-τg(∂Ut)(q)m·∂q-pτ+(∂Ut)(q)·∂p+gτmβ∂q2+mτβ∂p2.

Thus, if ft is the instantaneous density of (1), then(DV)t(x)=(∂t+Lx)Vt(x)

is the mean forward derivative of V along the paths of (1) and by definitionEPVtf(xtf)-Vtι(xtι)=∫tιfdtEP(DV)t(xt).

This observation justifies the introduction of the value function as a Lagrange multiplier.

Our definition of the value function in EP omits the contribution from the variation of the Gibbs-Shannon entropy to the mean entropy production from the Pontryagin-Bismut functional. This is because the Gibbs-Shannon entropy in (5) is fully specified by the assigned boundary conditions and therefore does not enter the determination of the optimal control.

Variational Equations

We determine the optimal control equations by a stationary variation of (11). As expected, the variation with respect to the value function yields the Fokker–Planck equation for the probability density13 (∂t-Lx†)ft(x)=0.

The variation with respect to the probability density yields the dynamic programming equation [97]14 (∂t+Lx)Vt(x)+Ct(Ut)(x)=0.

In the overdamped case [5, 28, 58], and in the case when the control is a function of both position and momentum [68, 84, 99], the variation with respect to the potential yields a local, exactly integrable condition for the optimal control. In stark contrast, we find that the optimal control potential in the underdamped case must solve an integral equation coupled to the Fokker–Planck and dynamic programming equations [34]:15 ∂q·∫Rdddpft(q,p)gτm∂qVt(q,p)+∂pVt(q,p)=∂q·(f~t(q)bt(q))

with f~t(q) as the position marginal of ft(q,p) (see Eq. (126) ) andbt(q)=βτ(1+g)2m((∂Ut)(q)-(∂U⋆)(q))[KL]βgτm2(∂Ut)(q)+(∂lnf~t)(q)β.[EP]

Finding regular extremals amounts to finding the simultaneous solutions of Eqs. (13), (14) and (15). The integro-differential stationary condition (15) is hard to approach due to its non-local nature (in momentum space). These issues are to some extent reminiscent of the Vlasov–Poisson–Fokker–Planck (see e.g. [61]) and the McKean–Vlasov (see e.g. [100]) equations. The condition somewhat simplifies when the configuration space is one dimensional. We can write16 ∫Rdpft(q,p)f~t(q)gτm∂qVt(q,p)+∂pVt(q,p)=βτ(1+g)2m((∂Ut)(q)-(∂U⋆)(q))[KL]βgτm2(∂Ut)(q)+(∂lnf~t)(q)β.[EP]

Dual Expression of the Optimal Cost

When the dynamic programming equation (14) holds, the Pontryagin-Bismut functional (11) reduces to17 A[f,U,V]|d.p.=∫R2dd2dx(Vtι(x)fι(x)-Vtf(x)ff(x)).

The optimum value of the cost hence coincides with the minimum, or at least infimum of (17), taken over all value functions satisfying the dynamic programming equation. This observation is the basis for the aforementioned duality relation used in [92], and later in [58, 68, 80, 93]. In what follows, we use (17) to compute the expression of minimum costs predicted by multiscale perturbation theory.

Gaussian Case

In view of the complexity of the optimal control condition (15), it is instructive to analyze the case of Gaussian boundary conditions. A similar analysis was performed in [34] for a case closely related to EP, but only in one dimensional configuration space.

Gaussian boundary conditions lead to major simplifications. The structure of the Fokker–Planck and dynamic programming equations preserve the space of Gaussian probability densities and second order polynomials in phase space for any at most quadratic controlUt(q)=ut+ut·q+12q⊤Utq

and referenceU⋆(q)=u⋆+u⋆·q+12q⊤U⋆q

potentials. In the above expressions ut, u⋆ are vectors in Rd and Ut, U⋆ are d×d real symmetric matrices. Thus, the probability density is fully specified by the set of first and second order cumulantsQt=EP(qt⊗qt)-EPqt⊗EPqtCt=EP(qt⊗pt)-EPqt⊗EPptPt=EP(pt⊗pt)-EPpt⊗EPpt.

Here and below, we use ⊗ to denote the outer product of vectors in Rd. Correspondingly, a value function of the form18 Vt(q,p)=vt+vt(q)·q+vt(p)·p+q⊤Vt(q,q)q+p⊤Vt(p,p)p+p⊤Vt(p,q)q+q⊤Vt(q,p)p2

satisfies the dynamic programming equation (14). In (18), Vt(q,q) and Vt(p,p) are d×d symmetric matrices, andVt(q,p)=Vt(p,q)⊤.

The Fokker–Planck equation (13) reduces to a closed system of differential equations of first order for the second order cumulants of the Gaussian statistics19 ddtQt=Ct+Ct⊤m-gτmUtQt+QtUt+2gτmβ1ddtCt=-1τCt-gτmUtCt-UtQt+1mPtddtPt=-2τPt-UtCt-Ct⊤Ut+2mβτ1

and to a system of differential equations of first order for the first order cumulants sustained by the solution of second order ones:20 ddtEPqt=EPptm-gτmut+UtEPqtdd tEPpt=-EPptτ+ut+UtEPqt.

The full system of cumulant equations (19)-(20) is complemented by boundary conditions at both ends of the control horizon:Q0=Qι&Qtf=QfC0=Ctf=0P0=Ptf=mβ1d

andEPqtι=qι&EPqtf=qfEPptι=EPptf=0.

The boundary conditions can be satisfied because the potential couples the cumulant equations to a first-order differential system of equal size for the coefficients of the value function in (18).

Analysis of Case KL

For case KL, we get 21a ddtVt(q,q)=(UtVt(p,q)+Vt(q,p)Ut)+gτmUtVt(q,q)+Vt(q,q)Ut-βτ(1+g)2mUt-U⋆Ut-U⋆

21b ddtVt(q,p)=Vt(q,p)τ+UtVt(p,p)+gτmUtVt(q,p)-Vt(q,q)m

21c ddtVt(p,p)=2τVt(p,p)-Vt(p,q)+Vt(p,q)⊤m

and 22a ddtvt(q)=Vt(p,q)⊤ut+Utvt(p)+gτmUtvt(q)+Vt(q,q)u-βτ(1+g)2mUt-U⋆ut-u⋆

22b ddtvt(p)=vt(p)τ-vt(q)m+Vt(pp)ut+gτmVt(pq)u.

Finally, we find23 ddtvt=vt(p)·ut+gτmvt(q)·ut-TrmβτV(p,p)+gτmβV(q,q)-(1+g)βτ4mut-u⋆2.

The structure of (21)–(22) is analogous to that of the cumulant equations. The coefficients of second degree monomials in (18) satisfy a closed system whose solution sustains the equation for the coefficients of first order monomial.

We now turn to the solution of (15). A straightforward exercize in Gaussian integration yields the explicit expression of the “osmotic force” [88] or “score function” [11] of the position marginal24 (∂lnf~t)(q)=-Qt-1(q-EPqt)

as well as an explicit expression for the conditional expectation25 EP(pt∣qt=q)=EPpt+CtQt-1(q-EPqt).

Upon inserting (24), (25) into (15) and matching the coefficients of monomials of same degree in q, we arrive at the equations 26a Qt-1Ut+UtQt-1=Qt-1Mt+Mt⊤Qt-1

26b TrMt-Ut=0

withMt=U⋆+2m(1+g)τβVt(p,p)CtQt-1+Vt(p,q)+2gτ(1+g)βVt(q,p)CtQt-1+Vt(q,q)

and the dependent conditions27 ut=u⋆+(U⋆-Ut)EPqt+2m(1+g)βτvt(p)+Vt(p,q)EPqt+Vt(p,p)EPpt+2g(1+g)βvt(q)+Vt(q,q)EPqt+Vt(q,p)EPpt.

Clearly, the conditions (27) are always satisfied if (26) is solvable. In fact we recognize that equation (26a) is in fact a Lyapunov equation. Uniqueness, symmetry and positivity of the solution are very well understood [62]. In particular, for every t∈[tι,tf] we can write the solution of (26a) as28 Ut=∫0∞dse-Qt-1sQt-1Mt+Mt⊤Qt-1e-Qt-1s.

The solution is well defined because by definition Qt is a positive matrix. Finally, taking the trace of both sides of (28) readily recovers (26b) thus completing the proof that the Gaussian case is solvable.

Analysis of Case EP

The equations that change are (21a)ddtVt(q,q)=(UtVt(p,q)+Vt(q,p)Ut)+gτmUtVt(q,q)+Vt(q,q)Ut-2βτgmUtUt,

Eq. (22a) which is replaced byddtvt(q)=Vt(p,q)⊤ut+Utvt(p)+gτmUtvt(q)+Vt(q,q)u-2gτβmUtut

and, finally, Eq. (23) which for the mean entropy production readsddtvt=dτ+vt(p)·ut+gτmvt(q)·ut-TrmβτV(p,p)+gτmβV(q,q)-gβτmut2-1βTrUt.

A qualitative difference with case KL occurs for vanishing g when the mean entropy production does not explicitly depend upon the control potential. This is most evident when inspecting (15). We get 29a gQt-1Ut+gUtQt-1=Qt-1M~t+M~t⊤Qt-1

29b TrM~t-gUt=0

where nowM~t=g2Qt-1+m2τβVt(p,p)CtQt-1+Vt(p,q)+gτ2βVt(q.p)CtQt-1+Vt(q,q)

and30 gut=-gUtQt-1EPqt+m2βτvt(p)+Vt(p,q)EPqt+Vt(p,p)EPpt+g2βvt(q)+Vt(q,q)EPqt+Vt(q,p)EPpt.

Whereas for g>0 the optimal potential is uniquely determined by the solution of the Lyapunov equation (29a), the limit g↓0 is singular. The Lyapunov equation becomes a constraint imposed on the coefficients of the value function. The upshot is that for vanishing g it is not possible to satisfy boundary conditions imposed on all phase space cumulants. In other words, the problem is not generically solvable for a generic assignment of Gaussian probability densities (2), (3). The problem admits, however, a solution if boundary data are just the position marginals. A detailed slow manifold analysis performed in the one-dimensional case in the supplementary material of [34] shows that the equation for g equal zero coincides with the slow manifold equations (see e.g. [101]) of the limit g↓0 optimal control equations. This gives a precise mathematical meaning to the idea of δ-Dirac optimal control upheld in [21]. It also shows that even if the optimal control does not exist for g equal zero, the strictly positive lower bound on the mean entropy production is always in agreement with (9).

General Case in One Dimension

A distinctive trait of the underdamped extremal equations (13), (14) and (15) is the integral term in Eq. (15), which introduces a non-local condition in the momentum variable. This is in stark contrast with the overdamped counterpart of (15). Indeed the latter is exactly integrable and thus reduces the extremal conditions to a pair of local hydrodynamics equations [5, 28]. In this section we construct a systematic multiscale expansion of (13) - (15) around the overdamped limit. By proceeding in this way we manage to reabsorb the non-locality in phase space into effective parameters of local equations—the cell problem—in configuration space. We perform our analysis in two-dimensional phase space. Extension to higher dimensional phase space is possible at the price of dealing with far more cumbersome algebra.

The approach we follow is inspired by [65, 66]. The first step is to project the momentum dependence in Eqs. (13)–(15) onto the basis of Hermite polynomials orthonormal with respect to the Maxwell thermal equilibrium distribution. We obtain an kinetic-theory-type hierarchy of coupled equations that do not depend on the momentum. Despite the additional complication of dealing with an infinite number of equations, this description turns out to be the ideal starting point for a multiscale expansion approach (in time).

Non-dimensional Variables

In order to neaten our notation, it is expedient to preliminary introduce non-dimensional variables:t=tτ,q=qℓ,p=βmp

where ℓ is the typical length-scale set by the mechanical potentials in the boundary conditions.

Next, we introduce the non-dimensional counterparts of the phase space density, value function and mechanical control potential:ft(q,p)=ℓmβft(q,p)Vt(q,p)=Vt(q,p)Ut(q)=βUt(q).

In non-dimensional variables, the generator of the phase space process  (12) becomes31 Lx=-(p-∂p)∂p+εp∂q-ε(∂qUt)∂p-ε2g(∂qUt)-∂q∂q=τLx

where now the order parameter of the overdamped expansion32 ε=τ2βℓ2m

explicitly appears. Equipped with these definitions, we rewrite the Fokker–Planck 33a ∂t-Lx†ft=0,

the dynamic programming33b ∂t+LxVt=-ε2(1+g)∂qUt-∂qU⋆24[KL]1-p2ε2g((∂qUt)2-∂q2Ut)[EP],

and the stationary condition equations33c ∫Rdpft(q,p)ft(0)(q)∂p+εg∂qVt(q,p)=ε(1+g)2(∂qUt(q)-∂qU⋆(q))[KL]εg∂q(2Ut(q)+lnft(0)(q))[EP],

whereft(0)(q)=∫Rdpft(q,p)

denotes the position marginal of the probability density function.

Expansion in Hermite Polynomials

Calling Hn the n-th Hermite polynomial (see Appendix D for details), we expand the probability density and the value function as34 ft(q,p)=e-p222π∑n=0∞ft(n)(q)Hn(p)

and35 Vt(q,p)=∑n=0∞vt(n)(q)Hn(p).

The expansion coefficients ft(n) and Vt(n) are scalar functions of the position and time. At equilibrium, the expansion for the probability density consists of the term n=0 only. The remaining contributions are non-zero only in out-of-equilibrium conditions. In particular, all ft(n) and Vt(n) for n>0 vanish at the beginning and at the end of the control horizon because of the boundary conditions (2) and (3). The expansion in Hermite polynomials turns the extremal equations (33) into an infinite hierarchy of equations whose n-th elements are: 36a ∂t+nft(n)+ε(∂q-∂qUt)ft(n-1)+ε(n+1)∂qft(n+1)=ε2g∂q(∂qUtft(n))+∂q2ft(n)

36b ∂t-nvt(n)+ε(n+1)(∂q-(∂qUt))vt(n+1)+ε∂qvt(n-1)=-δn,0ε2(1+g)4∂qUt-∂qU⋆2[KL]-δn,0gε2(∂qUt)2-∂q2Ut-δn,2[EP]

36c ∑n=0∞n!(n+1)ft(n)vt(n+1)+εg∂qvt(n)=ε(1+g)2ft(0)∂qUt-∂qU⋆[KL]εgft(0)∂q(2Ut+lnft(0)).[EP]

More detail on the derivation of the above equations is given in Appendix D. The hierarchy is complemented by equilibrium boundary conditions on the probability density, that, in the non-dimensional variables, read: 37a ftι(0)(q)=exp(-Uι(q))∫Rdyexp(-Uι(y))

37b ftf(0)(q)=exp(-Uf(q))∫Rdyexp(-Uf(y))

37c ftι(n)(q)=ftf(n)(q)=0n≥1.

Multiscale Perturbation Theory

The hierachy (36) is equivalent to the original extremal equations (33), and holds for any value of ε. We are interested in cases where ε≪1 in order to solve (36) with a perturbative strategy. The limit of vanishing ε is, however, singular and cannot be handled by regular perturbation theory. We therefore resort to multiscale perturbation theory. In doing so, we need to take into account an essential difference with respect to the multiscale treatment of the overdamped limit of the underdamped dynamics [65, 66]. The difference is that the mechanical potential is not assigned but must be determined by solving the stationary conditions (36c). In addition, we are dealing with a time-boundary value problem rather than with an initial data problem. To overcome these difficulties, we formulate the multiscale expansion drawing from the Poincaré–Lindstedt technique [63] and renormalization group ideas that in recent years have been successfully applied to the resummation of perturbative series arising from Hamiltonian and dissipative dynamical systems [67]. Our strategy is based on the following considerations.We suppose that the time variation of all functions in the hierarchy (36) occurs through effective time variables 38 tj=εjt,j≥0.

occasioned by the overdamped order parameter ε. As a consequence, the partial derivative with respect to t breaks down into a differential operator 39 ∂t=∂t0+ε∂t1+ε2∂t2+⋯

thus introducing a new dynamical variable at each order of the overdamped expansion.

We assume that the mechanical potential has a finite limit when ε tends to zero. This assumption [33, 34] is central in order to recover the overdamped dynamics [5, 28]. The a priori justification of the assumption is that momentum marginals of the boundary conditions already describe a Maxwell thermal equilibrium. The need for a controlled dynamics only arises in consequence of the boundary conditions imposed on the position process. In the generator (31), the mechanical potential is coupled to the dynamics by the overdamped expansion order parameter ε. This fact leads to the inference that the control potential should admit a regular expansion in ε as a function of the position variable, varying in time on scales set by ε.

The Poincaré–Lindstedt method is usually formulated for initial value problems. In such a context, the dynamics of the slow times tj with j>0 is fixed by canceling secular terms (equivalently: resonances), i.e. polynomial terms in the time variable which as times increases would lead to a breakdown of perturbation theory. Such a secular term subtraction scheme is equivalent to a renormalization group type partial resummation of the perturbative expansion [67]. We need to adapt the subtraction scheme to a boundary problem. At ε equal zero the Fokker–Planck hierarchy (36a) is decoupled from the dynamic programming one (36b). As a consequence, the boundary conditions (37) cannot be satisfied at zero order of the regular perturbative expansion. We therefore use the boundary conditions to determine the slow time dependence of the f(n).

The value function expansion coefficients f(n) are not subject to anything other than satisfying the dynamics (36a). The logical basis for the resonance subtraction scheme is the duality relation (17). We reason that a cost can only be generated on the same time scales over which the mechanical control potential varies. We thus require that the dependence of f(0) must be constant with respect to the fastest time t0. We also observe that, although physically motivated, a non-uniqueness is intrinsic in any secular term cancellation or finite renormalization scheme exactly because these techniques involve a partial and not a complete resummation of the perturbative expansion [64]. Consistent alternative schemes may differ by higher order terms in the regular perturbative expansion.

The introduction of the slow time variables (38) is justified under a sufficiently wide scale separation. In principle, the perturbative expansion only holds for ε≪1 and tf≫1 (i.e. tf≫τ). Yet, we hope that extrapolating the results for finite control horizons will give sufficiently accurate results if the resummation scheme correctly captures the “turnpike behavior” of the exact solution of the optimal control. Turnpike behavior means the tendency of optimal controls to approximate the solution of the adiabatic limit, corresponding to a vanishing cost, as much as possible. We refer to [68] for further discussion and references on this point.

In summary, our aim is to look for a solution of (36) in the form of multiscale power series 40a ft(n)(q)=∑i=0∞εift0(n:i)(q)

40b vt(n)(q)=∑i=0∞εivt0(n:i)(q)

40c Ut(q)=∑i=0∞εiUt0(i)(q),

where each addend of the above series depends, a priori, on all time scales41 tj=(tj,tj+1,tj+2,...).

Results

We report the main results of the overdamped multiscale expansion, while deferring their derivation to Sect. 6.2. Without loss of generality, we settι=0

to neaten the notation. Within second order in ε in the multiscale expansion, the solution of the Fokker–Planck equation takes the form42 ft(q,p)=ft0,t2(0:0)(q)+εpft0,t2(1:1)(q)e-p222π+ε2ft0,t2(0:2)(q)+pft0,t2(1:2)(q)+(p2-1)ft0,t2(2:2)(q)e-p222π+e-p222πO(ε3).

We emphasize that t0=t and t2=ε2t and that Eq. (42) is independent of t1=εt. We also neglect slower time scales tj, j>2, as they only provide higher order corrections. Hence, for all the results presented in this Subsection, we drop the explicit dependence on t1 and t3.

Cell Problem Equations

As customary [72], we refer to the secular term subtraction conditions emerging at order O(ε2) in regular perturbation theory as the cell problem. Secular term subtraction fixes the functional dependence upon the slow time t2. As a consequence, we find it expedient to denote the unknown quantities of the cell problem as43 ρt2(q)=ft0,t2(0:0)(q)

and as an auxiliary field σt2(q) related to vtf,t2(0:0)(q) and ft0,t2(0:0)(q) by equation (101) in Sect. 6.2 below. The formulation of the cell problem in terms of the pair ρt2, σt2 exactly recovers the optimal control equations governing the overdamped limit in KL and EP: 44a ∂t2ρt2=∂q(ρt2∂qσt2+α∂qρt2)

44b ∂t2σt2=12(∂qσt2)2-α∂q2σt2+α2∂q2U⋆-(∂qU⋆)22-α2χt2q

44c χt2=BA∫Rdqρt2(q)∂q∂q2U⋆(q)-(∂qU⋆(q))22.

We fully specify the cell problem by complementing (44) with the exact boundary conditions imposed by the position marginals of (2), (3) and written in non-dimensional variables as in (37)45 ρt2(q)|t2=0=e-Uι(q)∫Rdye-Uι(y)ρt2(q)|t2=ε2tf=e-Uf(q)∫Rdye-Uf(y).

In (44) we implyU⋆(q)=0

for EP. Depending upon the problem under consideration, the constant α takes the values46 α=(1+g)A[KL][0.3cm]0[EP],

while the constants A and B are given by (see Eq. (98) below)47 A1+g=1-(ω2-4)tanhωtf2tanhtfωtfωtanhωtf2-2tanhtfB1+g=-tanhωtf2tfωtanhtf-2tanhωtf2ωtanhωtf2-2tanhtf

with48 ω=1[KL]1+gg.[EP]

A and B always admit a finite limit as g tends to zero: by (48) the limit of vanishing g entails ω tending to infinity in case EP. Furthermore, they depend upon the size of the control horizon tf so that49 limtf↗∞A-1-g=limtf↗∞B=0.

When U⋆=0, the cell problem reduces to a coupled system of a Fokker–Planck and Burgers’ equation 50a ∂t2ρt2=∂q(ρt2∂qσt2+α∂qρt2)

50b ∂t2σt2=12(∂qσt2)2-α∂q2σt2.

By (49) and (46) in the limit of infinite control horizon (tf↗∞), we recover the result of [57] that in the overdamped limit optimal entropic transport is a viscous regularization of the minimization of the mean entropy production. As the optimal control of the latter problem [28] is equivalent to optimal mass transport, we also recover Mikami’s result [56]. In fact, for tf finite but sufficiently large to justify the scale separation required by the multiscale approach, the cell problem (50) allows us to extract information about corrections to the overdamped limit.

Cumulants and Marginal Distribution

Solving the cell problem allows us to evaluate the leading order corrections to the overdamped limit of all phase space cumulants within order O(ε2). Namely, all cumulants turn out to be linear combinations of functionals of the pair ρt2 and σt2, weighed by pure functions of the fast time t0.

We denote the non-dimensional counterparts of the position and momentum processes with a tildeq~t=qτtℓ,p~t=βmpτt.

Unlike the cell problem, the cumulants are functions of both the fast t0 and slow t2 time variables. To neaten the expressions, we denote the moments of the position process with respect to the probability density specified by the cell problem as51 μt2(n)=∫Rdqρt2(q)qn,n=1,2.

Correspondingly, we also writeμ˙t2(n)=∂t2∫Rdqρt2(q)qn.

In particular,μ˙t2(1)=-∫Rdqρt2(q)(∂σt2)(q)

is a constant, i.e.μ¨t2(1)=0,

for both the optimal entropic transport (KL with zero reference potential) and minimum mean entropy production EP problems. We justify this claim in Appendix F.

Momentum mean. Recalling (42), the expectation value of the momentum process conditioned on the position process is52 EPp~t|q~t=qρt(q)=∫Rdppft(q,p)=εft0,t2(1:1)(q)+O(ε2).

In section 6.2, we show how to compute ft0,t2(1:1) from the solution of the cell problem. We obtain53 ft0,t2(1:1)=-at0ρt2∂(σt2+αlnρt2)A+ρt2Bat0-Abt0A(A-B)μ˙t2(1).

Here at0 and bt0 are pure functions of the fast time t0: 54a at0=1+sinh(ωt0)tanhωtf2-cosh(ωt0)+bt0

54b bt0=ωe-tfcosh(ωt0)-e2t0ωcoshtf-2sinhtfcothωtf2+ωe-tfe2tf-cosh(ωtf)sinh(ωt0)ωcoshtf-2sinhtfcothωtf2sinh(ωtf).

It is straightforward to verify thata0=atf=b0=btf=0

hence enforcing the boundary conditions imposed on (53). The derivation of an explicit expression of this quantity requires the solution of the O(ε4) cell problem in the same way as (44a) specifies f(1:1).

Equipped with the above definitions, by integrating (52) in q we arrive at55 EPp~t=εat0-bt0A-Bμ˙t2(1)+O(ε2).

In order to interpret this result, we recall a standard result of multiscale analysis (see e.g. § 2.5.1 of [72]) ensuring that∫0fdtEPp~t≈ε∫0fdt0at0-bt0A-B∫0ε2tfdt2μ˙t2(1)

when the separation of scales is sufficiently large: tf≫O(1) with ε2tf=O(1). The relation (98) between the integral over the functions (54) and the constants A and B (47) implies that as the duration of the control horizon grows, the momentum expectation tends to∫0fdtEPp~t≈tf≫1εtfμε2tf(1)-μ0(1)1+g.

Once recast into dimensional quantities, the identity reads56 ∫0fdtEPpt≈tf≫ττtfβμε2tf(1)-μ0(1)ℓ(1+g).

Correction to the position marginal distribution. Upon integrating out the momentum variable in (42), we get57 ∫Rdpft(q,p)=ρt2(q)+ε2ft0,t2(0:2)(q)+O(ε3).

In Sect. 6.2 we show that58 ft0,t2(0:2)=-g∂qft0,t2(1:1)+(1+g)t0tf∫0fds-∫00ds∂qfs,t2(1:1).

Inspection of (58) reveals that the marginal (57) exactly satisfies the non-perturbative boundary conditions (45) and preserves normalization within accuracy. We avail ourselves of (57) to evaluate the remaining linear and second order cumulants.

Position mean. We readily obtain59 EPq~t=μt2(1)+ε2gat0-bt0A-Bμ˙t2(1)-ε2(1+g)μ˙t2(1)A-Bt0tf∫0fds-∫00ds(as-bs)+O(ε3)

As expected, at the boundaries of the control horizon the cumulant are fully specified by the boundary conditions and so are independent of ε.

Position-momentum cross correlation. After straightforward algebra we find60 EPq~tp~t-EPq~tEPp~t=εat0ς˙t22A+O(ε2)

whereςt2=μt2(2)-(μt2(1))2

and61 ς˙t2=-2∫Rdqρt2q∂qσt2+αlnρt2-2μ˙t2(1)μt2(1)=∂t2μt2(2)-(μt2(1))2≡∂t2ςt2.

Position variance. We obtain its expression by evaluating the difference betweenEPq~t2=μt2(2)+2ε2g∫Rdqqft0,t2(1:1)(q)-2ε2(1+g)∫Rdqqt0tf∫0fds-∫00dsfs,t2(1:1)(q)

and the squared mean valueEPq~t2=(μt2(1))2+2ε2gat0-bt0A-Bμt2(1)μ˙t2(1)-2ε2(1+g)μt2(1)μ˙t2(1)A-Bt0tf∫0fds-∫00ds(as-bs)+O(ε3).

After some algebra, we find that the expression of the variance reduces to62 EPq~t2-EPq~t2=ςt2+ε2gat0Aς˙t2-ε21+gAς˙t2t0tf∫0fds-∫00dsas+O(ε3)

with ςt2 and ς˙t2 defined by (61).

Momentum variance. The expectation value of the squared momentum conditioned on the position process isEp~t2|q~t=qρt2(q)=∫Rdpp2ft(q,p)=ρt2(q)+ε2ft0,t2(0:2)(q)+2ft0,t2(2:2)(q)+O(ε3).

with63 ft0,t2(2:2)=ft0,t2(1:1)22ρt2-ρt2∫00dse-2(t0-s)∂qft0,t2(1:1)ρt2.

After some tedious algebra we arrive at64 Ep~t2-(Ep~t)2=1-ε2at02A2(μ˙t2(1))2+2ε2∫00dse-2(t0-s)asA∫Rdqρt2∂q2(σt2+αlnρt2)+ε2at02A2∫Rdqρt2(∂q(σt2+αlnρt2))2+O(ε3).

We notice that the variance satisfies the boundary conditions in consequence of the identity∫0fdse-2(t0-s)as=0

which follows from (92).

Optimal Control Potential

Similarly, the cell problem yields the leading order expression for the gradient of the optimal control potential65 ∂qUt(q)=-∂qlnρt2(q)-a˙t0+at0A∂qct2(q)-(Ba˙t0-Ab˙t0)+(Bat0-Abt0)A(A-B)μ˙t2(1)+O(ε)

witha˙t0=∂t0at0

and66 ct2(q)=-σt2(q)-αlnρt2(q)

the current potential of the cell problem. In other words, the gradient of (66) is the current velocity [88] which allows us to represent (44a) as a mass conservation equation for any strictly positive α.

By definition, the current velocity vanishes when the system is in a Maxwell-Boltzmann equilibrium state. Hence, finite time transitions at minimum cost are not between Maxwell-Boltzmann equilibrium states, as we see from the explicit expression of the drift at the end times∂qU0(q)=-∂qlnρ0(q)-a˙0A∂qc0(q)-Ba˙0-Ab˙0A(A-B)μ˙0(1)+O(ε)

and∂qUε2tf(q)=-∂qlnρε2tf(q)-a˙ε2tfA∂qcε2tf(q)-Ba˙ε2tf-Ab˙ε2tfA(A-B)μ˙ε2tf(1)+O(ε).

From the physics point of view, this means transitions minimizing thermodynamic cost functionals have non-vanishing current velocity at the start and end of the protocol. Mathematically, this is unsurprising because the boundary conditions associated to the optimal control problem do not impose any conditions on the terminal values of the control potentials.

For all practical purposes, the shape of potential corresponding to the boundary equilibrium states can be matched at zero cost, through an instantaneous change of the control.

Minimum Cost

We evaluate the expression for the minimum cost using the duality relation (17). KL : The projection onto Hermite polynomials couches (17) into the form K(P‖Q)=∫Rdqf0(0)v0(0)-ftf(0)vtf(0).

At leading order, multiscale perturbation theory yields the approximation K(P‖Q)=∫Rdqf0,0(0:0)vtf,0(0:0)-f0,ε2tf(0:0)vtf,ε2tf(0:0)+O(ε2).

This is because the non-perturbative boundary conditions only allow contributions that are proportional to f0,0(0:n). In addition, we subtract secular terms in the value function expansion by requiring vtf,t2(0:2)=v0,t2(0:2).

Thus, in our multiscale framework, the value of vtf,ε2tf(0:2) can be only determined by higher order cell problems. To gain insight into the predicted features of the minimum, we couch the optimum value of the divergence into the form K(P‖Q)=-∫0ε2tfdt2∫Rdq∂t2f0,t2(0:0)vtf,t2(0:0)+O(ε2).

The above representation allows us to express the divergence in terms of the cell problem density (43) and the identity vtf,t2(0:0)=σt2+(α-A)lnρt22A-U⋆2-Bqμ˙t2(1)2A(A-B)-B4A(A-B)∫02ds(μ˙s(1))2

stemming from (89) and (101) in Sect. 6.2. Indeed, straightforward algebra yields 67 K(P‖Q)=∫0ε2tfds∫Rdqρs(∂q(σt2-αU⋆))24A+A-α2A∫Rdqρε2tflnρε2tfρ⋆-ρ0lnρ0ρ⋆+B4A(A-B)∫0ε2tfds(μ˙s(1))2+O(ε2)

where lnρ⋆=-U⋆-ln∫Rdqe-U⋆

and A-B=(1+g)1-2tanhtf2tf

which is positive definite when tf>2. In (67), all terms but the first vanish in the limit of infinite scale separation tf tending to infinity. Further elementary considerations shed more light on the sign of the corrections. Recalling (66) and the properties of the current velocity, we obtain the identity 68 ∫Rdqρt2(∂q(σt2-αU⋆))2=∫Rdqρt2(∂qct2)2+α2∫Rdqρt2∂qlnρt2ρ⋆2+2α∂t∫Rdqρt2lnρt2ρ⋆.

We then re-write (67) as K(P‖Q)=14α∫0ε2tfdt2∫Rdqρt2(∂q(σt2-αU⋆))2-B4A(A-B)∫0ε2tfdt2∫Rdqρt2(∂qct2)2-(μ˙t2(1))2+α-(A-B)4α(A-B)∫0ε2tfdt2∫Rdqρt2(∂qct2)2+αα-A4A∫0ε2tfdt2∫Rdqρt2∂qlnρt2ρ⋆2+O(ε2).

The identity μ˙(1)=∫Rdqρt2(∂qct2)

and the Cauchy–Schwarz inequality then ensure that all corrections are positive for tf>2. Thus for any tf sufficiently large to ensure a separation of time scales, we arrive at the inequality K(P‖Q)≥14α∫0ε2tfdt2∫Rdqρt2(∂q(σt2-αU⋆))2

whence we read the multiscale perturbation theory prediction of the Talagrand-Otto-Villani constant CTOV in (10). To do so, we focus on entropic transport and set U⋆ to zero. Next, we recall that any solution of the cell problem (50a)–(50b) enjoys the lower bound [102] 69 ∫0ε2tfdt2∫Rdqρt2(∂qσt2)2≥EP~q~ε2tf-q~02ε2tf

where P~ is the measure generated by the overdamped Schrödinger bridge in [0,ε2tf] associated to the stochastic differential equation 70 dq~t2=-(∂σt2)(q~t2)dt2+2αdwt2.

The inequality (69) is a consequence of the law of iterated expectation (see e.g. [76] pag. 310). Indeed, it ensures that ∫0ε2tfdt2∫Rdqρt2(∂qσt2)2≡∫0ε2tfdt2EP~((∂σt2)(q~t2))2=∫0ε2tfdt2EP~EP~((∂σt2)(q~t2))2|q~0≥∫0ε2tfdt2EP~(EP~(∂σt2)(q~t2)|q~0)2.

We now invert the order of integration and apply the Benamou–Brenier argument [89] to the stochastic paths generated by (70) and find ∫0ε2tfdt2∫Rdqρt2(∂qσt2)2≥EP~(EP~q~ε2tf-q~0-2αw~ε2tf|q~0)2ε2tf.

The inequality now follows because w~ is the Wiener process with respect to the measure P~ and as such has zero conditional expectation with respect to q~0. In Appendix E we present a path integral derivation of the same result. The upshot is that for entropic transport we get a Talagrand-Otto-Villani type inequality K(P‖Q)≥14αEP~q~ε2tf-q~02ε2tf.

In dimensional units, the same result reads K(P‖Q)≥14αβmEP~qtf-q02τtf.

EP : Upon contrasting (5) with (11), the exact expression of the minimum mean entropy production reads E=∫Rdqf0(0)(v0(0)+lnf0(0))-ftf(0)(vtf(0)+lnftf(0)).

The multiscale approximation then is E=-∫0ε2tfdt2∫Rdq∂t2f0,t2(0:0)(vtf,t2(0:0)+lnf0,t2(0:0))+O(ε2)

where, by (89) and (101), the identity vtf,t2(0:0)+lnf0,t2(0:0)=2σt2A+2Bqμ˙(1)A(A-B)+B(μ˙(1))2t2A(A-B)

with μ˙(1)=με2tf(1)-μ0(1)ε2tf

holds true. After some algebra, we arrive at 71 E=11+g∫0ε2tfdt2∫Rdqρt2(∂qσt2)2+1+g-A(1+g)A∫0ε2tfdt2∫Rdqρt2(∂qσt2)2-(μ˙(1))2+1+g-(A-B)4(1+g)(A-B)(μ˙(1))2ε2tf+O(ε2)

with A-B=(1+g)1-2tanhωtf2ωtf.

All addends in (71) are positive. Furthermore, the last two vanish both in the limit of infinite scale separation and upon recalling the definition (48) of ω when the coupling constant g is vanishing limg↘0E=∫0ε2tfdt2∫Rdqρt2(∂qσt2)2.

We also emphasize for case EP, the field σ satisfies the compressible Euler equation. As a consequence, we can directly apply the Benamou-Brenier inequality [89] to (71) and straightforwardly recover the bound (9).

Accuracy of the Multiscale Approximation

Infinite hierarchies of equations such as (36) appear in the study of Liouville’s and Boltzmann equations [65, 103]. Many numerical methods resort to a phenomenological truncation of the hierarchy. The multiscale method provides a controlled truncation at the level of second order equations. In fact, all cumulants up to second order can be reconstructed from an effective first order system embodied by the cell problem.

In Fig. 2, we summarize how the secular term cancellation (or, equivalently, solvability) conditions allow us to re-order contributions of the regular perturbative expansion within the hierachy. The upshot is that the predictions for cumulants and total cost obtained from the solution of the cell problem have different accuracies in ε.

Order-by-Order Solution

In this Section, we solve the hierarchy of equations (36) in a multiscale perturbative series in powers of ε. To this goal, we insert Eqs. (40a)–(40c) into Eq. (36), and identify equations of distinct order in the power expansion, taking into account the time differentiation, which acts on the multiscale dependence of the probability density and value function according to (39). The derivation of the results is briefly outlined in words below.

At order zero in ε, the equations for the density and value function give rise to two decoupled infinite systems of first order differential equations in the fast time t0. These systems are trivially integrable with respect to the fast time t0, implicitly keeping all information about the boundary condition in the unresolved dependence of the integration constants upon the slow times.

Remarkably, at order ε1 the boundary (37) and stationary conditions reduce the non-trivial contribution of the two infinite hierarchies of equations to a system of two first order differential equations in the fast time for f(1:1) and v(1:1). Dependence upon higher order coefficients of the expansion in Hermite polynomials enters these equations in the form of functions of the slow time t2 that must be determined at order ε2 in the regular perturbative expansion. As no secular term appears at this order we can assume within accuracy independence of the solution of the extremal equations from t1.

At order ε2, we can determine all unknown quantities inherited from lower orders in the regular perturbative expansion by imposing the cancellation of secular terms. This fixes the dynamical dependence upon the slow time t2 in the form of a cell problem. We enforce the correct boundary conditions in terms of f(0:2), f(2:2) and v(0:2), v(2:2). Finally, if we set all the f(n:0), f(n:1) that are not sustained by the drift and all the v(n:0), v(n:1) that are not needed to control the non-vanishing contributions to the density to zero, it is self-consistent to setf(1:2)=v(1:2)=0.

Figure 2 is a stylized summary of the procedure. Additional details are provided in Appendix F.

In principle, it is possible to extend the analysis to orders higher than ε2, as done in [65]. The appearance of spatial derivatives of higher order than the second may, however, call for the introduction of appropriate variables to perform partial resummations [103]. We return to this point in Sect. 9Fig. 2 Scheme of the multiscale approach presented in the text. The logical order of the calculation is represented by the black solid arrows, going through the sequence of solutions of the differential systems at the different orders. The quantities computed with this strategy are reported in the coloured boxes, where different colours correspond to different steps of the order-by-order multiscale calculation. Dashed arrows show the functional dependencies of the computed quantities. ⋆The calculation actually shows that there is no dependence on the t1-time scale

Boundary Conditions

The boundary conditions (2), (3) are by hypothesis independent of the Stokes time and therefore remain the same once expressed in non-dimensional units. Consequently, all ft0(n:i)’s with n≥1 vanish at the boundaries, so that 72a f0,t1(n:i)(q)=f0(0)(q)δn,0δi,0

72b ftf,t1(n:i)(q)=ftf(0)(q)δn,0δi,0,

where δi,j is a Kronecker delta. Without loss of generality, we set tι=0.

The non-perturbative boundary behavior is not assigned a priori but is determined by that of the probability density. However, in multiscale perturbation theory, we have the freedom to choose how partial resummations to cancel secular terms are performed [63]. We have reasoned that contributions to the cost can only come from the same time scale as those where the control varies, which gives the following resonance subtraction condition73 vtf,t1(0:i)(q)-v0,t1(0:i)(q)=vtf(0)(q)-v0(0)(q)δi,0.

Solution of the Problem at Order Zero

The calculation starts at order zero of the ε-expansion. Equation (36a) can be written at order zero in ε as:∂t0+nft0,t1(n:0)=0,

which implies74 ft0(n:0)=ct1e-nt0.

Here ct1 is fixed by imposing the initial condition at time t0=0:ft0(n:0)|t0=0=δn,0f0,t1(0:0),

following from Eq. (72). This observation leads to75 ft0(n:0)=δn,0f0,t1(0:0).

Solving the value function equation (36b) at order zero in ε gives76 vt0(n:0)-vtf,t1(n:0)et0-tf=0[KL]δn,2(1-e2(t0-tf))/2.[EP]

This time we have no boundary conditions to impose. However, it follows from Eq. (36c) that77 vt0(1:0)=0,

hence78 vt0(n:0)-(1-δn,1)vtf,t1(n:0)et0-tf=0[KL]δn,2(1-e2(t0-tf))/2.[EP]

The t0-dependence of the probability density and the value function is completely determined at order zero. We have no way to enforce the boundary condition at t=tf for ft0(n:0) at this stage: we will need to impose it on a slower time scale, in this way exploiting the additional freedom provided by the multiscale approach.

Solution of the Problem at Order One

By expanding Eq. (36a) at order one in ε, one gets∂t0ft0(n:1)+∂t1ft0(n:0)+nft0(n:1)+n+1∂qft0(n+1:0)+∂q+∂qUt0(0)ft0(n-1:0)=0.

The boundary conditions for the probability density force all terms of order higher than zero in ε vanish at t=0 and t=tf. This is a consequence of our assumption that the protocol starts and ends in equilibrium states, which cannot depend on the relaxation time scale ε. They must coincide with the stationary states of the overdamped limit ε→0. The n=0 case of Eq. (6.2.3), by recalling Eq. (75), implies therefore that ft0(0:0) is independent of t1, hence79 f0(0:0)=f0,t1(0:0)=f0;t2(0:0),

where we have introduced the notation80 ft0;t2≡ft0,0,t2.

Similarly, the equations with n≥2 lead to81 f0(n:1)=0,n≥2.

The case n=1 is less trivial and brings about the relation82 ∂t0ft0(1:1)+ft0(1:1)+∂qf0;t2(0:0)+f0;t2(0:0)∂qU0(0)=0.

Similarly, the value function equation (36b) at order one, for the case n=1, gives83 ∂t0-1v0(1:1)=-∂qv0(0:0)-2∂q+(∂qUt0(0))v0(2:0).

Once complemented with a condition for the drift ∂qU0(0), Eqs. (82) and (83) form a closed system of differential equations. The missing relation can be obtained from the stationarity condition (36c), which at order one in ε reads84 gf0;t2(0:0)∂qv0(0:0)+f0;t2(0:0)v0(1:1)+2ft0(1:1)v0(2:0)=1+g2ft0;t2(0:0)∂qU0(0)[KL]2gft0;t2(0:0)∂qU0(0)+g∂qft0,t2(0:0).[EP]

Eq. (84) provides an expression for the drift, which can be inserted into Eq. (82) to obtain a relation for v0(1:1) (see Eq. (140) in Appendix F). The system is then solved by differentiating the resulting equation with respect to t0, and eliminating ∂t0v0(1:1) through (83) and v0(1:1) through (140). A second-order ODE for ft0(1:1) is found:85 ∂t02ft0(1:1)-ω2ft0(1:1)=Ft0,

with ω as defined in (48). The dependence of Ft0(q) on t0 is known; for its explicit expression, see Eq. (141). The equation can be solved by recalling that the Green function for the second order differential equation (85) is86 Gt,s=Jt,s+Js,t

withJt,s=-θ(t-s)sinhω(tf-t)sinhωsωsinhωtf,

with θ(·) being the Heaviside step-function. By introducing the notation87 Gt(k)=∫0fdsGt,se-k(tf-s),

one obtains for ft0(1:1) the relation88 ft0;t2(1:1)=ω2G0(0)f0;t2(0:0)∂qζt2+41+g∂qvtf;t2(2:0)f0;t2(0:0)G0(2)[KL]1g∂qvtf;t2(2:0)f0;t2(0:0)-f0;t2(0:0)2G0(2)[EP]

where89 ζt2(q)=2vtf;t2(0:0)(q)+U⋆(q)+lnf0;t2(0:0)(q)[KL]12vtf;t2(0:0)(q)+lnf0;t2(0:0)(q),[EP]

and a notation analogous to (80) is adopted also for the value function. The last step of the solution of the order ε1 consists in writing the optimal control potential as a function of (88). From Eq. (82) we find90 ∂qUt0(0)(q)=-∂qf0;t2(0,0)(q)+∂t0ft0;t2(1:1)(q)+ft0;t2(1:1)(q)f0;t2(0,0)(q).

Since there are no equations for ∂t1ft0(1:1) nor ∂t1v0(1:1), i.e. no secular terms are found on the time scale t1, we can assume the solution to be independent of t1. Once ft0(1:1) is known, an explicit expression for v0(1:1) can also be found from (83) (see Eq. (140) in Appendix F). From the above equation it is possible to derive the expression (65) for the optimal drift, by using the expressions of ft0;t2(0:0) and ft0;t2(1:1) that will be found in the next subsection.

Solution of the Problem at Order Two

The order-two expansion of (36a) provides the following relations 91a ∂t0ft0(0:2)+∂t2f0;t2(0:0)=-∂qft0;t2(1:1)+g∂q2f0;t2(0:0)+g∂qft0;t2(0:0)∂qUt0;t2(0)

91b ∂t0ft0(1:2)+ft0(1:2)=-∂qft0;t2(0:1)-ft0;t2(0:0)∂qUt0(1)

91c ∂t0ft0(2:2)+2ft0(2:2)=-∂qft0;t2(1:1)-ft0;t2(1:1)∂qUt0;t2(0)

91d ∂t0ft0(n:2)+nft0(n:2)=0,n>2.

The last equation (91d) ensures that all terms ft0(n:2) with n>2 vanish once equilibrium boundary conditions are taken into account. Equation (91b) provides a relation for ft0(1:2) that requires knowledge of ∂qU0(1). If we expand the stationary condition (36c) to second order in ε and assume that all vt0(n:0), vt0(n:1) that are not needed to control the non-vanishing ft0(n:2)’s can be set to zero, we get∂qU0(1)=2g+1vt0(1:2)+2ft0(1:2)vt0(2:0)ft0(0:0)[KL]ω2-12vt0(1:2)+2ft0(1:2)vt0(2:0)ft0(0:0).[EP]

We insert this result into (91b) and the corresponding equation for vt0(1:2), and after straightforward, albeit tedious, algebra we arrive at∂t0ft0(1:2)-ω2ft0(1:2)=0.

Taking into account the boundary conditions, we getf0,t1(1:2)=ftf,t1(1:2)=0

and we conclude that for any t0,ft0,t1(1:2)=0.

The same applies to vt0,t1(1:2).

Let us focus first on Eq. (91c). Integrating over t0, one hasf0(2:2)=-∫00dse-2(t0-s)∂qfs;t2(1:1)+fs;t2(1:1)∂qUs;t2(0).

By substituting the expression of the drift obtained from Eq. (90) and integrating the term proportional to ∂sfs;t2(1:1) by parts, we findft0;t2(2:2)=ft0;t2(1:1)22f0;t2(0:0)-f0;t2(0:0)∫00dse-2(t0-s)∂qfs;t2(1:1)f0;t2(0:0),

which is Eq. (63). This relation implies, recalling the boundary conditions, that92 ∫0fdse-2(tf-s)∂qfs;t2(1:1)f0;t2(0:0)=0.

By substituting Eq. (88), an equation for the term ∂qvtf;t2(2:0)f0;t2(0:0) can be derived (see Eq. (142) in Appendix F). Once plugged back into Eq. (88) itself, it yields93 ft0;t2(1:1)(q)=-f0;t2(0:0)at0(∂ζt2)(q)-bt0κt2.

Here we introduce the functions whose explicit expression we gave in (54)at0=-ω2G0(0)-G(0:2)G(2:2)G0(2)bt0=ω2G(0:2)G(2:2)G0(2).

By (87) the two functions at0 and bt0 are non homogeneous solution of the unstable oscillator equation weighed by constant coefficient also depending upon integrals over the Green function94 G(k:l)=∫0fdse-l(tf-s)Gs(k).

In (93) we also introduce95 κt2=∫Rdqft0;t2(0:0)(q)(∂ζt2)(q),

where we use the function ζt2 defined in Eq. (89). Equation (93) will be crucial in the following, as it allows to write a closed system of differential equations for f0;t2(0:0) and vtf;t2(0:0), which can be reshaped as in Eqs. (44).

Taking into account the boundary conditions and Eq. (90), Eq. (91a) can be integrated over t0 to give96 ∂t2f0;t2(0:0)+g+1tf∫0fds∂qfs;t2(1:1)=0.

If we now substitute Eq. (93) we get97 ∂t2f0;t2(0:0)=A∂qf0;t2(0:0)(∂qζt2)-Bκt2∂qf0;t2(0:0)

where 98a A=g+1tf∫0fdsas=-ω2(1+g)tfG(0:0)-G(0:2)2G(2:2)

98b B=g+1tf∫0fdsbs=ω2(1+g)tfG(0:2)2G(2:2).

The above relations lead to Eq. (47).

We now need to find an equation for ζt2 in order to close the differential system and find the t2-dependence of f0;t2(0:0). To this aim, we consider the case n=0 for the expansion of Eq. (36b) at order two in ε. It reads99 ∂t0vt0;t2(0:2)+∂t2vtf;t2(0:0)+∂q-(∂qUt0;t2(0))vt0;t2(1:1)+g∂qvtf;t2(0:0)=-g+14∂q(Ut0;t2(0)-U⋆)2[KL]-g(∂qUt0;t2(0))2-∂q2Ut0;t2.[EP]

We integrate the above equation over t0. By substituting (90) and making repeated use of Eqs. (85), (93) and (92) (see Appendix F for details), one finds100 ∂t2ζt2=A-B2(∂qζt2)2+B2∂qζt2-κt22+α2AW⋆+∂q2f0;t2(0:0)(f0;t2(0:0))2-12∂qf0;t2(0:0)f0;t2(0:0)2,

which is the closure equation for ζt2. The constant α is defined by Eq. (46).

The differential system for f0;t2(0:0) and ζt2 can be rewritten in a much more convenient form by introducing the auxiliary field101 σt2(q)=Aζt2(q)-αlnf0;t2(0:0)(q)-Bqκt2+A-B2∫02dsκs,t32.

Indeed, taking into account Eq. (149), it is easy to verify that Eqs. (97) and (100) are amenable to the form (44). Let us stress that, in terms of the field σ, Eq. (97) becomes102 κt2=1A-B∫Rdqft0;t2(0:0)(q)(∂σt2)(q).

Upon inserting this identity,  (101), and (51), in (93), we recover the expression (53).

Finally, by plugging Eqs. (96) and  (90) in Eq. (91a) one gets:∂t0ft0(0:2)-1+gtf∫0fds∂qfs;t2(1:1)=-(1+g)∂qft0(1:1)-g∂t0∂qft0(1:1).

Integrating over t0 leads to Eq. (58).

Analytic Results for the Gaussian Case

As discussed in Sect. 6.1 and shown analytically in Sect. 6.2, in order to find the explicit solution of the optimal problem, one first needs to address the differential system (44). In most cases, the solution can only be found numerically: this is discussed in the next Section. However, if the assigned initial and final conditions are Gaussian probability density functions (meaning that the particle is subject to harmonic confinement), the solution can be found analytically.

To do this, we plug a Gaussian ansatz for the density and a parabolic one for σt2, namely 103a ρt2=12πςt2exp-(q-μt2(1))22ςt2

103b σt2=σ2(0)+σ2(1)q+σ2(2)q2

where μ(1) and μ(2) are consistent with Eq. (51), into Eqs. (44).

Next, we solve for the coefficients, taking into account the boundary conditions. The derivation is straightforward and not carried out here. For both cases KL and EP, and U⋆=0, the explicit expressions for the relevant coefficients appearing in Eqs. (103) areμt2(1)=μ0(1)+t2ε2tfμε2tf(1)-μ0(1)ςt2=(t2-ε2tf)2ς0+t22(ε2tf-t2)λε2tf+t2ςε2tfε4tf2σt2(1)=μ0(1)ςε2tf-ε2tfα+λε2tf-με2tf(1)ς0+ε2tfα+λε2tfε4tf2ςt2

whereλε2tf=ε4tf2α2+ς0ςε2tf

and for ς˙t2=∂t2ςt2104 σt2(2)=2α-ς˙t22ςt2

Knowing these coefficients allows us to compute the cumulants discussed in Sect. 6.1.2 for the general case.

It is worth noticing that these results result in a remarkably simple expression for mean entropy production at g=0E=με2tf(1)-μ0(1)2+ςε2tf-ς022ε2tf.

WhenU⋆=U⋆(1)q+12U⋆(2)q2

(104) remains valid, whereas it is possible to close the hierachy with a second order equation for the variance of the position process105 2ςt2ς¨t2-ς˙t22-4α2U⋆(2)ςt2-1=0

The general solution of this equation takes the form106 ςt2=c1e2αU⋆(2)t2+c2e-2αU⋆(2)t2+c3U⋆(2)

with the constants ci, i=1,2,3, related by the algebraic equation4c1c2+1-c32=0

Unfortunately resolving the ci’s in terms of generic boundary conditions leads to somewhat cumbersome expressions. In Sect. 8.1.1 we consider a special case of particular relevance.

Numerically Assisted Applications

In this section, we apply numerical methods to the multiscale expansion to analyze the underdamped dynamics, both in the case of Gaussian boundary conditions and in more complex boundary conditions, in particular, those modelling Landauer’s one bit of memory erasure.

In the Gaussian case, we have a system of differential equations specifying the non-perturbative solution, and we can therefore use numerical integration to solve the associated boundary value problems, from which we can obtain the first and second order phase space cumulants. We then use the perturbative approach to compute the same values, which show good agreement: see Figs. 3 for case KL and  4 for case EP. Additionally, we take a look at expansion and compression with Gaussian boundary conditions in Sect. 8.1.1.

Furthermore, the perturbative approach can be used to make predictions for the cumulants when no analytic solution is available. We demonstrate this using boundary conditions modelling Landauer’s one bit of memory erasure, as illustrated in Fig. 1. This requires numerically solving the cell problem (50), from which we obtain the optimal control protocol and the marginal distribution of the position in the overdamped dynamics. We can then compute leading order corrections to approximate the quantities in the underdamped dynamics.

Gaussian Case

In cases KL and EP, when the boundary conditions are assigned as Gaussian random variables, we have two boundary value problems for the first and second order cumulants. For case KL, we compute approximate solutions to the systems (19) and (20), and, for case EP, we make the amendments as described in Sect. 4.2.

The perturbative approach follows Sect. 7, and instead we have only one boundary value problem. The dependant quantities: momentum mean, momentum variance and the position-momentum cross correlation, as well as the higher order corrections to the position mean and variance can then be computed.

The respective boundary value problems are integrated numerically using the DifferentialEquations.jl [104] library in the Julia programming language. The results of the perturbative and non-perturbative integrations for case KL are in Fig. 3 and for case EP are in Fig. 4. In both case KL and case EP, we see that the perturbative expansion gives a very good approximation of the true solution.Fig. 3 Position mean (a) and variance (b); position-momentum cross correlation (c); and momentum mean (d) and variance (e) for the underdamped problem minimising the Kullback–Leibler divergence (Case KL) from a free diffusion with assigned Gaussian initial and final conditions. We plot the expressions of cumulants up to second order predicted by our perturbative approach, solid dark blue line, and contrast them with the numeric solution of the corresponding non-perturbative equation of Sect. 4.1 shown with a dashed light blue line. To plot the perturbative expressions we use the exact solution of the cell problem given in Sect. 7 and then we determine the full perturbative prediction for each of the cumulants using the expressions given in Sect. 6.1.2. We impose Gaussian boundary conditions through the first and second order cumulants of the position and momentum (tildes denoting non-dimensional coordinates): at initial time t=0, we set the position and momentum variance Q~0=P~0=1, the position-momentum cross correlation C~0=0; and the position and momentum means EPq~tι=EPp~tι=0. At final time t=tf: the momentum variance P~tf=1; the position-momentum cross correlation C~tf=0; the position variance Q~tf=1.7; the momentum mean EPp~tf=0; and the position mean EPq~tf=2. We use tf=5, ε=0.2 and g=0. Numerical integration is performed by a fourth order co-location method in the DifferentialEquations.jl library [104]

Fig. 4 Position mean (a) and variance (b); position-momentum cross correlation (c); and momentum mean (d) and variance (e) for the underdamped problem minimising the entropy production EP from a free diffusion with assigned Gaussian initial and final conditions. We use tf=5, ε=0.2, and compare the values computed by the perturbative approach (solid lines) with the numeric solution of the non-perturbative system at g=0.5 (blue, dashed) and g=0.1 (green, dotted), and ω=(1+g)/g. We impose Gaussian boundary conditions through the first and second order cumulants of the position and momentum. We use the same values for the boundary conditions and the same numerical integration method as in Fig. 3. We compute the perturbative predictions in the same way as explained in the caption of Fig. 3. The equations of the non-perturbative system are those of Sect. 4.2

Asymmetry in Optimal Approaches to Equilibrium

Fig. 5 Thermal kinematics in the underdamped dynamics. We compare a compression (blue) and expansion (orange) process starting from thermodynamically equidistant states from the final state. We fix the final position variance as Qtf=4, the corresponding reference potential U⋆=1/4, and β=1 in all panels. Panels (a)-(c) show the picture where the initial variances are chosen to be close together. We use Q0(c)=4.20 for compression and Q0(e)≈3.8065 for expansion. For (a)-(b), we use tf=5 and ε=0.1. Panel (a) shows the variance of the position, with the solid lines computed non-perturbatively and the overlayed dashed lines computed using the linearized approach outlined in Sect. 8.1.1. Panel (b) shows the static Kullback–Leibler divergence (107), with the inset axes showing the difference between the compression and expansion process. Panel (c) shows the difference in the dynamic Kullback–Leibler divergence as a function of the time horizon tf, with ε=0.5. Panels (d)–(e) illustrate the case when the initial variances are chosen to be further apart. We use Q0(c)≈10.3467 for compression and Q0(e)=1 for expansion. We use tf=5 and ε=0.1. Panel (d) shows the variances and (e) shows the static Kullback–Leibler divergence as functions of the time interval. Panel (f) shows the dynamic Kullback–Leibler divergence as a function of the time horizon tf, with ε=0.5. First order cumulants do not play a role in the analysis and are equal to 0 in all panels. All numerical integration is performed by DifferentialEquations.jl [104], as in Fig (3)

Very recently, [70, 71] highlighted the existence of a cooling versus heating asymmetry in the relaxation to a thermal equilibrium from hotter and colder states that are “thermodynamically equidistant”. Although not strictly a distance, the Kullback–Leibler divergence from the thermal state may be used to identify the dual processes [70]. We show that a similar asymmetry also occurs in optimally controlled isothermal compressions versus expansions of a small system.

To this goal we make the following observations. Choosing a reference potential U⋆ in (4) equal to the potential in the final condition (3) forces the current velocity specified by the optimal protocol to be as small as possible at the end of the control horizon. In this sense, the optimal control problem models a relaxation to a thermal equilibrium in finite time. Well-established laboratory techniques [37, 105, 106] use the fact that the optical potential generated by a laser to trap a colloidal nanoparticle is effectively Gaussian. We combine these two observations to compare the compression versus the expansion of a nanosystem in an isothermal environment when the initial data are thermodynamically equidistant from the final equilibrium state. Mathematically, this means that the position marginals of the boundary conditions (2), (3) are centered Gaussians that differ only in the variance. In such a case, the only non-trivial optimal control equations are (19) and (21). Our aim is to compare a compression and an expansion process starting from “dual” initial states. Duality is with respect the Kullback–Leibler divergence from the end state whose value is initially the same for the two opposite processes. In the notation of Sect. 4 the Kullback–Leibler divergence for d=1 reads107 K(f~t‖f~tf)=12QtQtf-1-lnQtQtf

We fix the terminal conditionQtf(i)=(βU⋆)-1=:Q⋆,i=e,c

and compare the evolution of probability densities specified by initial conditions at tι=0Q0(e)<Q⋆<Q0(c)

such that the initial position marginals have equal Kullback–Leibler divergence from the final stateK(f~0(e)‖f~⋆)=K(f~0(c)‖f~⋆)

The dynamic Schrödinger bridge with boundary conditions (Q0(e),Q⋆) / (Q0(c),Q⋆) provides a model of optimal expansion/compression of the system towards the equilibrium state characterized by Q⋆.

The multiscale prediction for the position variance isQtℓ2=ςε2tτ-ε2ς˙ε2tτAttf∫0tfτds-∫0tτdsas+O(ε3)

where for the sake of simplicity we set g=0. To relate non-dimensional quantities to their dimensional counterparts, we explicitly write the Stokes time τ and the typical length-scale ℓ of the transition. We suppose that the variance of the non-dimensional cell problem at the beginning of the control horizon isς0=vU2,withU2=βU⋆ℓ2

How much the non-dimensional constant v differs from unity controls the thermodynamic distance from the final state. In such a case we find that the coefficients ci’s in (106) arec1=yc2c2=2e2αε2tfU2ye4αε2tfU2+1ye4αε2tfU2-12c3=-1+6e4αε2tfU2y+e8αε2tfU2y2ye4αε2tfU2-12

withy=2cosh2αε2tfU2-3v+1e4αε2tfU2-2e2αε2tfU2+v-22sinhαε2tfU22v+cosh2αε2tfU2-1v+1e4αε2tfU2-2e2αε2tfU2

The above expressions are exact and provide a useful benchmark for exact numerical integration of the cumulant hierarchy (see Fig 5).

For transitions describing small deformations of the position marginal of the system, it is however expedient to resort to simpler approximated expressions. We obtain these by linearizing (105) around the final condition of the transition. In other words, we look for a solution of the formςt2=1U2+ςt2′+⋯

with dots corresponding to higher order terms in the non-linearity. We obtainςt2′=-v-1e-2αt2U2e4αt2U2-e4αε2tfU2U2e4αε2tfU2-1

This expression allows us to analytically compare the behavior of the divergence from a common end state of system undergoing an expansion and a compression. We see that if we choosev(e)=1-η

for η∼O(10-1) then within O(10-4) accuracy the initial data for the dual compression process isv(c)=1+η+2η23+4η39+44η4135+O(η5)

A straightforward calculation then shows that within leading order accuracyK(f~t(c)‖f~⋆)-K(f~t(e)‖f~⋆)≥0,∀0≤t≤tf

The result holds analytically for small deformations of the potential η≪1 and close to the overdamped limit ε≪1.

Another thermodynamic indicator encoding similar information is the cost of the dynamic Schrödinger bridge (4). This quantity is a global indicator of the transition that can be studied versus the duration of the horizon. Consistently with the analytic perturbative result the evaluation of (4) shows that the divergence from equilibrium is larger for compression processes. The difference between compression and expansion tends to zero as the duration of the horizon tends to infinity, thus indicating symmetry restoration for adiabatic processes.

Our findings are summarized in Fig. 5. Our analysis is in line with the findings of [70]. If we interpret the divergence from equilibrium at any fixed time as an indirect quantifier of the speed with which the system ultimately thermalizes, our analytic and numerical results confirm that expansion is faster than compression for Gaussian models.

Landauer’s Erasure Problem

Fig. 6 Predictions for the position mean (a) and variance (b); position-momentum cross correlation (c); momentum variance (d) and mean (e) in the underdamped dynamics minimizing the Kullback–Leibler divergence (Case KL) from a free diffusion computed using the perturbative expansion. The boundary conditions are assigned on the marginal density of the position: the initial state at t=0 is a single peaked distribution centered at xo=1 (see Eq. (108)) and the final is a double peaked distribution with peaks at -1 and 1 (see Eq. (109)). We set tf=5, ε=0.2, g=0, ω=1 and use α=(1+g)A≈0.64, where A is as defined in (98a). The functions ρt2(q) and σt2(q) are approximated numerically for values of t2∈[0,ε2tf] as the solution of Eq. (50) by a forward-backward iteration. We perform a total of 15 forward and backward passes of the iteration. At each step, the factors ϕ and ϕ^ are computed by means of Eq. (112) and normalized for numerical stability. We use a total of 8000 sample points for q in the interval [-6,6] and evolve 50000 independent Monte Carlo trajectories started from each q by an Euler-Maruyama discretization of the SDE (113) with step-size h=0.005. The forward-backward iteration is initialized by ϕ^tf with a vector of ones. We recover the functions ρ and σ from ϕ and ϕ^ using Eqns. (114); these are then smoothed using a convolution with a box filter with window size δ=0.06. The predictions for the underdamped moments are then computed using expressions found in Sect. (6.1.2)

We model the Landauer’s one bit of memory erasure [47] as a Schrodinger bridge problem between an initial state single-peaked distribution and final state as a double-peaked distribution, as illustrated in Fig. 1. We can make predictions for the first and second order cumulants of the position and momentum distributions from the perturbative expansion, by computing the numerical solution to the cell problem (50) and hence the appropriate corrections. We focus only on case KL.

We assign the initial and final state of the position marginal distributionρε2tι(q)=∫Rdpptι(q,p)=:Pι(q)ρε2tf(q)=∫Rdpptf(q,p)=:Pf(q)

where Pι and Pf denote the assigned initial and final distributions, and here take the explicit forms108 Pι(q)=1Zιexp(-(q-xo)4)

109 Pf(q)=1Zfexp(-(q2-xo2)2)

with Zι,Zf normalizing constants. The initial condition is a single peaked distribution centered at xo, and final condition is a double peaked distribution, with peaks at xo and -xo.

We look at the case of U⋆=0. The cell problem (44) can be approximated numerically using a forward-backward iteration. This specifically means computing the numerical solution of two coupled non-linear partial differential equations to obtain the functions ρ and σ of the slow time t2=ε2t.

We adopt the methodology of [84], beginning with the Hopf-Cole transformϕ^t2(q)=ρt2(q)expσt2(q)2α,ϕt2(q)=exp-σt2(q)2α

yielding a pair of Fokker–Planck equations 110a ∂t2ϕt2(q)+α∂q2ϕt2(q)=0

110b ∂t2ϕ^t2(q)-α∂q2ϕ^t2(q)=0

with coupled boundary conditions 111a ϕε2tf(q)=Pf(q)/ϕ^ε2tf(q)

111b ϕ^0(q)=Pι(q)/ϕ0(q)

In this form, the cell problem can be solved using the forward-backward iteration, an adaptation of Algorithm 1 of [84]. We make a slight simplification, in that we perform the numerical integration of equations (110) by a Monte Carlo method, computing 112a ϕ^ε2tf(q)=E(ϕ^0qε2tf|q0=q)

112b ϕ0(q)=E(ϕε2tfq0|qε2tf=q)

using the forward and backward evolution respectively of the underlying auxiliary (Ito) stochastic process113 dqt2=2αdwt2,

where {wt2}t2≥0 denotes a standard Wiener process. The values of qt2 are approximated with discretized trajectories of (113) by the Euler-Maruyama scheme.

The forward-backward iteration goes as follows: We begin by sampling a set of values for q from an interval on which both the initial and final assigned distributions Pι and Pf are compactly supported. We initialize the forward-backward iteration by taking a set of (positive) values for ϕ^ε2tf, which are then used to compute the boundary condition (111a) for equation (110a). We integrate equation (110a) using the expression (112b) to obtain ϕ0 and recompute ϕ^0 using (111b). By integrating (110b) using (112a) up to tf, we once again obtain ϕ^ε2tf. This procedure is then repeated until convergence; we verify that the boundary condition relations (111) are satisfied, and the mean-squared difference between two iterations of ϕtf and ϕ^tf is less than a specified tolerance. We can then recover the values of ρt2 and σt2 by the relations114 ρt2(q)=ϕ^t2(q)ϕt2(q)σt2(q)=-2αlogϕt2(q).

The optimal control protocol in the overdamped case is σ. From here, we use the relevant equations in Sects. 6.1.2 and 6.1.3 to make predictions for the first and second order cumulants of the position and momentum in the underdamped dynamics, which are shown in Fig. 6. The predicted marginal distribution of the position and the gradient of the optimal control protocol is shown in Fig. 7. Figure 8 contrasts the heights of the peaks of the marginal distribution of the position in the underdamped and overdamped dynamics over the time interval.Fig. 7 Predictions for the marginal density of the position (a)-(e) and (k) -(o) and the gradient of the optimal control protocol (f)-(j) and (p) -(t) in the overdamped (orange) and underdamped (blue) dynamics minimizing the Kullback–Leibler divergence (Case KL) from a free diffusion, U⋆=0. We show the distribution ρ and the gradient of the optimal control protocol -∂qσ (114) for the overdamped dynamics, and compute the corrections needed to obtain the corresponding quantities ft and -∂qUt for the underdamped dynamics. The shaded region in panels (a) and (o) show the assigned boundary conditions (108) and (109) respectively: the initial state at t=0 is a single peaked distribution centered at xo=1 and the final is a double peaked distribution with peaks at -1 and 1. We set tf=5, ε=0.2, g=0, ω=1 and α=(1+g)A≈0.64, where A is as defined in (98a). The functions ρt2(q) and σt2(q) are computed as in Fig. (6), and predictions for the underdamped are computed using the expressions in 6.1. All distributions are normalized

Fig. 8 Predictions for the heights of the peaks of the marginal distribution of the position in the overdamped (orange) and underdamped (blue) dynamics. Both peaks in the overdamped remain at roughly equal height, while the underdamped diverge. The distribution in the overdamped and the corrections to obtain the underdamped are computed as in Fig. 7, using tf=5, ε=0.2, and g=0

Conclusions and Outlook

In this paper, we address the problem of finding optimal control protocols analytically for finite time stochastic thermodynamic transitions described by underdamped dynamics. To such end, we introduce a multiscale expansion whose order parameter vanishes in the overdamped limit. Within second order accuracy, we are able to find corrections for the linear and quadratic moments of the process. When the boundary conditions are Gaussian, our results are in excellent agreement with the solutions found by non-perturbative numerical methods.

We expect our theoretical predictions to provide a necessary benchmark for design and interpretation of experiments on nanomachine thermodynamics. In particular, this is the case for statistical indicators of the momentum process, whose dynamical properties are a distinctive trait of the underdamped regime. Our predictions for the momentum variance and the position-momentum cross correlation are in qualitative agreement with the very recent experimental observations in related laboratory setups [40].

We envisage several directions to extend the present work. In our view, the most urgent and possibly relevant for applications is devising efficient numerical algorithms to determine regular extremals for general (non-Gaussian) boundary conditions. The non-local nature of the equations determining the regular extremals hamper the direct application of proximal algorithms [84, 100] and Monte Carlo methods. We address the problem of generalizing these methods to the underdamped case in a companion contribution [107]. Here, we also compute inertial corrections to the numerical solution of the overdamped problem [31] for minimal entropy production in Landauer’s problem [28].

A second main result of the present work is the proof that the optimal control for transitions between Gaussian states solve a Lyapunov equation in any number of dimensions. This is a strong indication of the existence of regular extremals in phase spaces of any number of dimensions: in view of [66], the extension of the multiscale method is very cumbersome, but otherwise conceptually straightforward. A more subtle issue is instead the computation of corrections of orders higher than two, which are prone to instabilities already at third order. Ideas motivated by normal form theory [103] offer a promising way to overcome this difficulty. Yet, the application to optimal control on a finite time horizon is still an open challenge.

From the physics perspective, the multiscale expansion appears best suited to deal with nanoscale dynamics when inertial effects are present, but are small in comparison to thermal fluctuations. A possible alternative approach is the underdamped expansion (see e.g. Chapter 6 of [76]). This technique could be used to extract complementary information to that obtained here.

In terms of applications, our results are relevant for all physical contexts where random fluctuations and inertial effects cannot be disregarded. This is the case, for example, in bit manipulation in electronic devices. Information bits are encoded using bi-stable states governed by double-well potentials. Inertia is required to improve the efficiency of most logic operations [44].

Our results find natural applications also in biophysics. The control of biological systems such as bacteria suspensions and swarms is nowadays accessible to experimentation through several techniques [108–111]. This has generated increasing interest in the theoretical challenge of applying control theory to active matter models, i.e. out-of-equilibrium dynamics showing complex phenomena inspired by biology [26, 112–114]. So far, however, only overdamped dynamics have been considered. While this does describe the behaviour of microscopic biological systems at high Reynolds numbers (e.g., bacteria in liquid suspensions) fairly, it is well known that inertial effects do play a fundamental role in some classes of such systems [115–117]. A meaningful description of the collective behaviour of flocks and swarms requires taking into account inertial effects that allow efficient propagation of information within the system [118–120]. Any approach to the control of these models should therefore be carried out in the underdamped regime: even if our results cannot be straightforwardly applied to collective dynamics, they may provide a promising starting point for the development of control theory in this context.

Appendices

Derivation of the Cost Functionals

The physics-style derivation of (4) and (5) proceeds by constructing finite dimensional approximation on families of time lattices tι≤t0≤⋯≤tN+1=tf with mesh sizeh=tf-tιN+2.

The one-step approximation of the transition probability density of (1) in the pre-point prescription is115 Tti+1,ti(h)(xi+1∣xi)=exp(-Ati+1,ti(h)(xi+1∣xi))Zi

where xi=qi⊕pi and116 Ati+1,ti(h)(xi+1∣xi)=βm4gτhqi+1-qi-pim-gτm(∂Uti)(qi)h2+βτ4mhpi+1-pi-piτ+(∂Uti)(qi)h2,

while Zi is a normalization constant irrelevant for the present considerations. Within accuracy (115) satisfies the Chapman–Kolmogorov equation [76]. Hence, we obtain the transition probability over any finite time interval by means of the limit117 Tt,t~(x∣x~)=limh↓0,Nh=t-t~∏i=1N∫d2dziTsi+1,si(h)(zi+1∣zi)Ts1,s0(h)(z1∣z0)

where we hold fixed in the limit t~=s0≤sN+1=t and z0=x~, zN+1=x. For any admissible potential, (117) satisfies by hypothesis the bridge boundary conditionsftf(x)=∫R2dd2dyTtf,tι(x∣y)ftι(y)

with ftι, ftf respectively assigned by (2) and (3).

Case KL

Proceeding in a similar fashion, the finite dimensional approximation of (4) is by definition118 K(PN||Q⋆N)=∫R2d(N+2)d2dx0d2dxN+1∏j=1Nd2dxjftι(x0)×Ttj+1,tj(h)(xj+1∣xj)ln∏k=1NTtk+1,tk(h)(xk+1∣xk)T⋆tk+1,tk(h)(xk+1∣xk),

where T⋆(h) is defined with respect to the reference potential U⋆. Using the properties of the logarithm and the normalization of the transition probability, the definition reduces to the sum119 K(PN||Q⋆N)=∑i=1N∫R2d(i+2)d2dx0d2dxi+1∏j=1id2dxj×ftι(x0)Ttj+1,tj(h)(xj+1∣xj)lnTtj+1,tj(h)(xj+1∣xj)T⋆tj+1,tj(h)(xj+1∣xj)

Next, we observe thatlnTti+1,ti(h)(xi+1∣xi)T⋆ti+1,ti(h)(xi+1∣xi)=βτ(1+g)h4m((∂U⋆)2(qi)-(∂Uti)2(qi))+β2qi+1-qi-hpim·((∂U⋆)(qi)-(∂Uti)(qi))+βτ2mpi+1-pi+hpiτ·((∂U⋆)(qi)-(∂Uti)(qi))

The outermost integrals in (119) over qi+1,pi+1 are Gaussian and equal to∫d2dxi+1Tti+1,ti(h)(xi+1∣xi)qi+1pi+1=qi+pim-gτm(∂Uti)(qi)h[0.3cm]pi-piτ+(∂Uti)(qi)h

We thus arrive at∫d2dxi+1Tti+1,ti(h)(xi+1∣xi)lnTti+1,ti(h)(xi+1∣xi)Tt1,t0(h)(xi+1∣xi)=βτ(1+g)h4m(∂U⋆)(qi)-(∂Uti)(qi)2

Inserting this result into (119) and passing to the continuum limit recovers (4).

Case EP

The starting point is (118) where we replace T⋆(h) with the transition probability generated by the backward stochastic differential equationsd♭qt=ptm+gτm(∂Ut)(qt)dt+2τgmβd♭wt(1)d♭pt=ptτ-(∂Ut)(qt)dt+2mτβd♭wt(2),

where the label ♭ recalls that the evolution proceeds backwards. The one-step approximation of the transition probability density on the lattice using the adapted post-point prescription yields120 Tti,ti+1♭(h)(xi∣xi+1)=exp(-Ati,ti+1♭(h)(xi∣xi+1))Zi♭

withAti,ti+1♭(h)(xi∣xi+1)=βm4gτhqi-qi+1+pi+1m+gτm(∂Uti+1)(qi+1)h2+βτ4mhpi-pi+1+pi+1τ-(∂Uti+1)(qi+1)h2.

We recover (5) by contrasting ratios of (115) and (120) over the same time intervals and by identifying the sum of two finite dimensional approximations of stochastic integrals over the same integrand but evaluated in the pre-point and post-point prescription as twice the same integral in the Stratonovich prescription [76]. We refer to [18] for the details of the calculation or, e.g., to [69] a derivation directly in the continuum limit using stochastic calculus and Girsanov formula.

Consistency of the Definition of Mean Entropy Production with Stochastic Thermodynamics

The calculation follows the same steps as [33]. LetHt(x)=p22m+Ut(q)

the kinetic plus potential Hamiltonian specified by the control potential in (5). We define the work done on the system during a time interval [0, t] asWt=∫0tds(∂sHs)(xs)≡∫0tds(∂sUs)(xs)

Here, xt ia a realization of (1). The partial derivative only affects the explicit time dependence of the control potential. Correspondingly, we identify the heat released by the system during the same realization with the Stratonovich stochastic integralQt=-∫0tdxs·(∂Hs)(xs),

so that for any path satisfying (1) the identity121 Ht=Wt-Qt

ensures the validity of the first law of thermodynamics. The dynamics (70) describes an open system in contact with an environment at constant temperature β-1. Hence the average heat released by the system in [0, t] also specifies the change of entropy in the environment122 ΔSt(e)=βEQt

The total entropy change is the sum of this quantity and the change of the Gibbs-Shannon entropy of the system [28]123 ΔSt=ΔSt(e)+Elnp0(x0)pt(xt)

In order to justify referring to (5) as mean entropy production, we need to show how it is related to (122). Indeed, from the properties of the Stratonovich integralβEQt=-βE∫0tdsvs(xs)·(∂Hs)(xs)

On the right hand side we introduce the current velocityvt(x)=J·(∂Ht)(x)-mτSg·(∂Ht)(x)+1β(∂lnpt)(x)

which is most conveniently written in terms of the 2d×2d real matricesJ=01d-1d0&Sg=gτ2m21d001d

Probability conservation and antisymmetry of J implyE∫0tds∂slnps(xs)=E∫0tds(∂Hs)(xs)J·(∂lnps)(xs)=0

and thus allow us to arrive atβEQt=Elnpt(xt)p0(x0)+mτE∫0tdsSg·(∂Hs)(xs)+1β(∂lnps)(xs)2

From this expression we readily verify that the thermodynamic mean entropy production (123) is positive definite. Next, we unfold the quadratic form in the integrandΔSt=mβτE∫0tds‖Sg·(∂Hs)(xs)‖2+mτE∫0tds2(∂H)(xs)·(Sg∂lnps)(xs)+1β‖Sg·(∂lnps)(xs)‖2

Under our working hypothesis that confining mechanical potentials produce probability density decreasing at infinity sufficiently fast, an integration by parts yields0≤E‖Sg·(∂lnps)(xs)‖2=-ETrSg(∂xt⊗∂xtlnps)(xs)

We are therefore in the position to apply the chain of identitiesElnpt(xt)p0(x0)=∫0tdsE(∂s+Lxs)lnps(xs)=mτ∫0tdsE-(∂Hs)(xs)·(Sg∂lnps)(xs)+1βTrSg(∂⊗∂lnps)(xs)

to couch (123) into the formΔSt=-Elnpt(xt)p0(x0)+mτE∫0tdsβ‖Sg·(∂Hs)(xs)‖2+(∂Hs)(xs)·(Sg∂lnps)(xs)

The last step is to make use of the kinetic plus potential form of the Hamiltonian. Straightforward algebra and an integration by parts yieldΔSt=E

Hence, the Kullback–Leibler divergence (5) coincides with the thermodynamic entropy production. The above calculation can also be conceptualized as a special case of the general theory expounded in [19].

Proof of the Mean Entropy Production Lower Bound

It is well known (see e.g. [19, 33]) that (5) can be couched into the explicitly positive form124 E=EP∫tιfdtmβτptm+1β∂ptlnft(qt,pt)2+EP∫tιfdtβgτm(∂Ut)(qt)+1β∂qtlnft(qt,pt)2

We add and subtract125 kt(q)=∫Rdddpft(q,p)f~t(q)pm

in the squared norm of the first integrand in (124) andht(q)=(∂Ut)(q)+1β∂qlnf~t(q)

to the second one. In both expressions we introduce the position marginal density126 f~t(q)=∫Rdddpft(q,p)

Upon expanding the norm squared into inner products, and taking advantage of the cancellation of the mixed term, we get127 E≥mβτEP∫tιfdtgτ2m2ht(qt)2+k(qt)2

For any g>0 and a, b arbitrary real numbers, the inequality128 (ga+b)2≤(1+g)(ga2+b2)

holds true. The upshot isE≥mβτ(1+g)EP∫tιfdtv~t(qt)2

The vector field appearing on the right hand-side of the inequality129 v~t(q)=kt(q)-gτmht(q)

is exactly the current velocity transporting the position marginal distribution:130 ∂tf~t(q)+∂q·(v~t(q)f~t(q))=0

We are therefore in the position to apply the Benamou–Brenier inequality [89]131 EP∫tιfdtv~t(qt)2≥EP∫tιfdtv~t(qt)2tf-tι=EPqtf-qtι2tf-tι

whence we finally recover (9).

Details of the Expansion in Hermite Polynomials

The Hermite polynomials are defined as132 Hn(p)=(-1)nep2/2dndpne-p2/2,

so thatH0(p)=1,H1(p)=p,H2(p)=p2-1,…

and so on. They fulfill the following orthonormality condition:133 Hn,Hm=∫Rdpe-p2/22πHn(p)Hm(p)=n!δn,m.

Let us notice that from the definition of Hermite polynomials it follows:134 ∂pHn(p)=nHn-1(p)

and135 pHn(p)=Hn+1(p)+nHn-1(p).

The above identities, together with decomposition (34), can be used to write136 p∂qft=pe-p222π∂q∑nft(n)Hn=e-p222π∑n∂qft(n)(Hn+1+nHn-1)=e-p222π∑n≥1∂qft(n-1)Hn+∑n(n+1)∂qft(n+1)Hn.

Similarly, one has137 ∂pft=-e-p222π∑nft(n)Hn+1

and138 ∂p2ft=e-p222π∑nft(n)Hn+2.

Substituting the above relations into the Fokker–Planck equation (33a) and projecting onto Hn we get139

which can be recast into Eq. (36a).

A similar approach can be followed for the value function. Recalling Eq. (35) one gets∂t-nvt(n)=-ε(n+1)(∂q-(∂qUt))vt(n+1)-ε∂qvt(n-1)-δn,0ε24(∂q(U⋆-Ut))2,

hence Eq. (36b) follows. Equation (36c) comes from an analogous expansion of Eq. (33c).

Path Integral Proof of the Inequality (69)

We start from the discretized stochastic differential equationqi+1-qi=bi(qi)h+2αηi+1h

where the label i runs over bins of the time discretization. We take uniform mesh h. b is a sufficiently regular drift, and the ηi’s are independent identically distributed centered Gaussian random variables with unit variance. We set out to computeE∫0fdtb(ξt)2=limh↓0Nh=tf∫∏k=0N+1dxkZhp(x0)∏j=0NTj(xj+1∣xj)∑k=0Nb(xk)h2

withTi(xi+1∣xi)=exp-(xi+1-xi-b(xi)h)24αh

and Zh a mesh dependent normalization constant. We emphasize the use of the pre-point prescription in our construction of finite dimensional approximations of the path integral. We perform the change of variablesyi=xi-xi-1-b(xi-1)hi≥1x0=y0

in consequence whereof the chain of identities∑i=0Nb(xi)h=∑i=0N(xi+1-xi-yi+1)=xN+1(yN+1,⋯,y0)-y0-∑i=0Nyi+1

holds true. As the change of variables in the pre-point prescription has unit Jacobian, we arrive atE∫0fdtb(ξt)2=limh↓0Nh=tf∫∏k=0N+1dykp(y0)Zh×e-∑i=1N+1yk24αhxN+1(yN+1,⋯,y0)-y0-∑i=0Nyi+12≥limh↓0Nh=tf∫dy0p(y0)(XN+1(y0)-y0)2

by the Cauchy inequality, withXN+1(y0)=∫∏k=1N+1dyke-yk24αhxN+1(yN+1,⋯,y0).

This inequality holds for any drift and therefore also for the one implementing the optimal bridge.

Further Details on the Order-by-Order Multiscale Expansion

In this Appendix, we present additional details about the order-by-order solution of the multiscale problem presented in Sect. 6.2. While it is not essential to follow the logic of our method, these intermediate steps may be a useful reference for the reader interested in the detailed derivation of the results.

In Sect. 6.2.3, we describe how to obtain a second order differential equation in t0 for ft0(1:1), namely Eq. (85). The first step is to find a relation for v0(1:1) by plugging Eq. (84) in Eq. (82). We get140 v0(1:1)=-g∂qvtf,t1(0:0)-1+g2∂qU⋆+lnf0;t2(0:0)+∂t0ft0;t2(1:1)f0;t2(0:0)+1+4vt0,t1(2:0)1+gft0;t2(1:1)f0;t2(0:0)[KL]-g∂qlnf0;t2(0:0)+vtf,t1(0:0)-2gf0;t2(0:0)∂t0ft0;t2(1:1)+1+vt0,t1(2:0)gft0;t2(1:1).[EP]

By differentiating in t0 and eliminating vt0,t1(1:1) and its time derivative through (140) itself and (83), one gets Eq. (85). The explicit expression of its right hand side reads141 Ft0=f0;t2(0:0)∂q2vtf,t1(0:0)+U⋆+lnf0;t2(0:0)+41+g∂qvt0,t1(2:0)f0;t2(0:0)[KL][0.3cm]ω22f0;t2(0:0)∂qvtf,t1(0:0)+lnf0;t2(0:0)+1g∂qvt0,t1(2:0)f0;t2(0:0)-∂qf0;t2(0:0)2g.[EP]

The dependence of the value function on t1 can be actually dropped, since no resonant equation holds for v0(0:0) and v0(2:0) on that time scale. This result allows us to find (88) through the use of the Green function (86):ft0;t2(1:1)=∫RdsGt0,sFs;t2.

Once an explicit expression for ft0;t2(1:1) is known (i.e., Eq. (88)), it can be substituted into Eq. (92) to get the relationG(0:2)ω2G(2:2)f0;t2(0:0)∂qζt2=4(1+g)∂qvtf;t2(2:0)f0;t2(0:0)+C1f0;t2(0:0)[KL][0.3cm]1g∂qvtf;t2(2:0)f0;t2(0:0)-12g∂qf0;t2(0:0)+C2f0;t2(0:0),[EP]

where C1 and C2 are constants that can be evaluated considering the limits for q→±∞. One then has142 ∂qvtf;t2(2:0)f0;t2(0:0)=-(1+g)4G(0:2)ω2G(2:2)f0;t2(0:0)∂qζt2-κt2[KL]-gG(0:2)ω2G(2:2)f0;t2(0:0)∂qζt2-κt2+f0;t2(0:0)2.[EP]

In Sect. 6.2.4 we derive the differential equation (44b), which allows closing the system, providing the t2-dependence of f0;t2(0:0). To this aim, one first needs to substitute the expression for vt0;t2(1:1) given by Eq. (99), obtaining143 ∂t0vt0;t2(0:2)+∂t2vtf;t2(0:0)-2∂q-(∂qU(0))vt0;t2(2:0)ft0;t2(1:1)f0;t2(0:0)=1+g2W⋆-W(0)[KL]g-∂qlnf0;t2(0:0)-W(0).[EP]

where we took into account Eq. (90) and we introducedW(0)=∂q2U(0)-12(∂qU(0))2

andW⋆=∂q2U⋆-12(∂qU⋆)2.

Because of the boundary conditions (73) one has144 ∫0fdt0vt0;t2(0:2)=0.

Besides, using (36b) at order 0,2vt0;t2(2:0)=∂t0vt0;t2(2:0)[KL]∂t0vt0;t2(2:0)+1[EP],

hence145 ∫0fdt0∂t0ft0;t2(1:1)+ft0;t2(1:1)2vt0;t2(2:0)ft0;t2(1:1)=0[KL]∫0fdt0ft0;t2(1:1)2[EP].

We also notice, recalling Eq. (85), that146 ∫0fdt0∂t0ft0;t2(1:1)2=-∫0fdt0ft0;t2(1:1)∂t02ft0;t2(1:1)=∫0fdt0ft0;t2(1:1)ω2ft0;t2(1:1)-Ft0

where Ft0 obeys Eq. (141). Taking into account Eqs. (90), (92), (144), (145) and (146), one obtains from Eq. (143) [KL]∂t2vtf;t2(0:0)=1+g2W⋆+∂q2f0;t2(0:0)(f0;t2(0:0))2-12∂qf0;t2(0:0)f0;t2(0:0)2+1+g4tf∫0fdt02∂qft0;t2(1:1)-ft0;t2(1:1)∂qζt2f0;t2(0:0)+∫0fdt0ft0;t2(1:1)e-2(tf-t0)tf(f0;t2(0:0))2∂qf0;t2(0:0)vtf;t2(0:0)[EP]∂t2vtf;t2(0:0)=1+gtf∫0fdt0∂qft0;t2(1:1)-ft0;t2(1:1)∂qζt2f0;t2(0:0)+∫0fdt0ft0;t2(1:1)e-2(tf-t0)tf(f0;t2(0:0))2∂qf0;t2(0:0)vtf;t2(0:0)-12.

Now we can substitute the explicit expressions for ft0;t2(1:1) and ∂qf0;t2(0:0)vtf;t2(0:0) provided by Eqs. (88) and (142). The integrals are evaluated by making repeated use of Eq. (94). We arrive at [KL]∂t2vtf;t2(0:0)-1+g2W⋆+∂q2f0;t2(0:0)(f0;t2(0:0))2-12∂qf0;t2(0:0)f0;t2(0:0)2=A2∂q2ζt2+A2∂qζt2lnf0;t2(0:0)-B2κt2lnf0;t2(0:0)+A4∂qζt22-B4κt22∂qζt2-κt2[EP]∂t2vtf;t2(0:0)=A∂q2ζt2+A∂qζt2lnf0;t2(0:0)-Bκt2lnf0;t2(0:0)-A∂qζt22-Bκt22∂qζt2-κt2.

By recalling Eq. (97), we obtain Eq. (100). It is now possible to compute the time derivative of κt2, by making use of (97) and (100):149 ∂t2κt2=∫Rdq∂t2f0;t2(0:0)∂qζt2+f0;t2(0:0)∂q∂t2ζt2=α2A∫Rdq∂q2U⋆f0;t2(0:0)∂qU⋆-∂qf0;t2(0:0).

The right hand side vanishes for case EP, and also for case KL when either U⋆ is a linear function of q (including the physically relevant case U⋆=0), or U⋆ is symmetric with symmetric boundary conditions. By taking into account this result and the definition (101), Eq. (100) is straightforwardly recast into Eq. (44b).

We finally observe that, recalling (51) and (102), one hasμ˙t2(1)=-(A-B)κt2.

Therefore,μ˙t2(1)=0

when the right hand side of (149) vanishes.

Acknowledgements

The authors are pleased to acknowledge discussions with Luca Peliti and Paolo Erdman. JS was supported by the Centre of Excellence in Randomness and Structures of the Academy of Finland and by a University of Helsinki funded doctoral researcher position, Doctoral Programme in Mathematics and Statistics. MB was supported by ERC Advanced Grant RG.BIO (Contract No. 785932).

Funding

Open access funding provided by Consiglio Nazionale Delle Ricerche (CNR) within the CRUI-CARE Agreement.

Data Availibility

Data sets generated during the current study are available from the corresponding author on reasonable request.

Declarations

Conflict of interest

The authors have no relevant financial or non-financial interests to disclose.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Julia Sanders, Marco Baldovin and Paolo Muratore-Ginanneschi have contributed equally to this work.
==== Refs
References

1. Schrödinger E Über die Umkehrung der Naturgesetze Sitzungsberichte der preussischen Akademie der Wissenschaften, physikalische mathematische Klasse 1931 8 9 144 153
Schrödinger, E.: Über die Umkehrung der Naturgesetze. Sitzungsberichte der preussischen Akademie der Wissenschaften, physikalische mathematische Klasse 8(9), 144–153 (1931). 10.1002/ange.19310443014
2. Chetrite R Muratore-Ginanneschi P Schwieger K E. Schrödinger’s 1931 paper “On the Reversal of the Laws of Nature” [“Über die Umkehrung der Naturgesetze”, Sitzungsberichte der preussischen Akademie der Wissenschaften, physikalisch-mathematische Klasse, 8 N9 144-153] Eur. Phys. J. 2021 46 1 1 29
Chetrite, R., Muratore-Ginanneschi, P., Schwieger, K.: E. Schrödinger’s 1931 paper “On the Reversal of the Laws of Nature’’ [“Über die Umkehrung der Naturgesetze’’, Sitzungsberichte der preussischen Akademie der Wissenschaften, physikalisch-mathematische Klasse, 8 N9 144-153]. Eur. Phys. J. 46(1), 1–29 (2021). 10.1140/epjh/s13129-021-00032-7
3. Cover, T.M., Thomas, J.A.: Elements of Information Theory, 2nd edn. Telecommunications and Signal Processing, p. 776. Wiley-Blackwell, Hoboken (2006). 10.1002/0471200611
4. Föllmer H Hennequin PL École d’Étè de Probabilitès de Saint-Flour XV-XVII Random Fields and Diffusion Processes 1988 New York Springer
Föllmer, H.: École d’Étè de Probabilitès de Saint-Flour XV-XVII. In: Hennequin, P.L. (ed.) Random Fields and Diffusion Processes. Springer, New York (1988). 10.1007/BFb0086180
5. Dai Pra P A stochastic control approach to reciprocal diffusion processes Appl. Math. Optim. 1991 23 1 313 329
Dai Pra, P.: A stochastic control approach to reciprocal diffusion processes. Appl. Math. Optim. 23(1), 313–329 (1991). 10.1007/BF01442404
6. Föllmer H Gantert N Entropy minimization and Schrödinger processes in infinite dimensions Ann. Probab. 1997 25 2 901 926
Föllmer, H., Gantert, N.: Entropy minimization and Schrödinger processes in infinite dimensions. Ann. Probab. 25(2), 901–926 (1997)
7. Léonard C A survey of the Schrödinger problem and some of its connections with optimal transport Discrete Contin. Dyn.l Syst. Ser. A 2014 34 4 1533 1574
Léonard, C.: A survey of the Schrödinger problem and some of its connections with optimal transport. Discrete Contin. Dyn.l Syst. Ser. A 34(4), 1533–1574 (2014). 10.3934/dcds.2014.34.1533. arXiv:1308.0215 [math.PR]
8. Chen Y Georgiou TT Pavon M Stochastic control liaisons: Richard Sinkhorn meets gaspard monge on a Schrödinger bridge SIAM Rev. 2021 2 249 313
Chen, Y., Georgiou, T.T., Pavon, M.: Stochastic control liaisons: Richard Sinkhorn meets gaspard monge on a Schrödinger bridge. SIAM Rev. 2, 249–313 (2021). 10.1137/20m1339982
9. Todorov E Optimality principles in sensorimotor control Nat. Neurosci. 2004 7 9 907 915 15332089
Todorov, E.: Optimality principles in sensorimotor control. Nat. Neurosci. 7(9), 907–915 (2004). 10.1038/nn130915332089
10. Todorov E Efficient computation of optimal actions Proc. Nat. Acad. Sci. 2009 106 28 11478 11483 19574462
Todorov, E.: Efficient computation of optimal actions. Proc. Nat. Acad. Sci. 106(28), 11478–11483 (2009). 10.1073/pnas.071074310619574462
11. Peyré G Cuturi M Computational optimal transport: with applications to data science Found. Trends Mach. Learn. 2019 5–6 355 607
Peyré, G., Cuturi, M.: Computational optimal transport: with applications to data science. Found. Trends Mach. Learn. 5–6, 355–607 (2019). 10.1561/2200000073. arXiv: 1803.00567
12. De Bortoli, V., Thornton, J., Heng, J., Doucet, A.: Diffusion schrödinger bridge with applications to score-based generative modeling. NeurIPS 2021 (spotlight) and arXiv: 2106.01357 (2021)
13. Vargas, F., Ovsianas, A., Fernandes, D., Girolami, M., Lawrence, N.D., Nüsken, N.: Bayesian Learning via Neural Schrödinger-Föllmer Flows. arXiv:2111.10510 (2021)
14. Kay ER Leigh DA Zerbetto F Synthetic molecular motors and mechanical machines Angew. Chem. Int. Ed. 2007 46 1–2 72 191
Kay, E.R., Leigh, D.A., Zerbetto, F.: Synthetic molecular motors and mechanical machines. Angew. Chem. Int. Ed. 46(1–2), 72–191 (2007). 10.1002/anie.200504313
15. Filliger R Hongler M-O Relative entropy and efficiency measure for diffusion-mediated transport processes J. Phys. A 2005 38 1247 1255
Filliger, R., Hongler, M.-O.: Relative entropy and efficiency measure for diffusion-mediated transport processes. J. Phys. A 38, 1247–1255 (2005). 10.1088/0305-4470/38/6/005
16. Peliti L Pigolotti S Stochastic Thermodynamics 2020 Princeton Princeton University Press
Peliti, L., Pigolotti, S.: Stochastic Thermodynamics. Princeton University Press, Princeton (2020)
17. Maes C The fluctuation theorem as a Gibbs property J. Stat. Phys. 1999 95 367 392
Maes, C.: The fluctuation theorem as a Gibbs property. J. Stat. Phys. 95, 367–392 (1999). 10.1023/A:1004541830999. arXiv: math-ph/9812015
18. Maes C Redig F Moffaert AV On the definition of entropy production, via examples J. Stat. Phys. 2000 41 3 1528 1554
Maes, C., Redig, F., Moffaert, A.V.: On the definition of entropy production, via examples. J. Stat. Phys. 41(3), 1528–1554 (2000). 10.1063/1.533195
19. Chetrite R Gawȩdzki K Fluctuation relations for diffusion processes Commun. Math. Phys. 2008 282 2 469 518
Chetrite, R., Gawȩdzki, K.: Fluctuation relations for diffusion processes. Commun. Math. Phys. 282(2), 469–518 (2008). 10.1007/s00220-008-0502-9. arXiv:0707.2725 [math-ph]
20. Schmiedl T Seifert U Optimal finite-time processes in stochastic thermodynamics Phys. Rev. Lett. 2007 98 108301 17358574
Schmiedl, T., Seifert, U.: Optimal finite-time processes in stochastic thermodynamics. Phys. Rev. Lett. 98, 108301 (2007). 10.1103/PhysRevLett.98.10830117358574
21. Gomez-Marin A Schmiedl T Seifert U Optimal protocols for minimal work processes in underdamped stochastic thermodynamics J. Chem. Phys. 2008 129 2 024114 18624523
Gomez-Marin, A., Schmiedl, T., Seifert, U.: Optimal protocols for minimal work processes in underdamped stochastic thermodynamics. J. Chem. Phys. 129(2), 024114 (2008). 10.1063/1.2948948. arXiv:0803.0269 [cond-mat.stat-mech]18624523
22. Esposito M Broeck CV Second law and Landauer principle far from equilibrium Europhys. Lett. 2011 95 4 40004
Esposito, M., Broeck, C.V.: Second law and Landauer principle far from equilibrium. Europhys. Lett. 95(4), 40004 (2011). 10.1209/0295-5075/95/40004. arXiv: 1104.5165
23. Aurell E Mejía-Monasterio C Muratore-Ginanneschi P Optimal protocols and optimal transport in stochastic thermodynamics Phys. Rev. Lett. 2011 106 25 250601 21770620
Aurell, E., Mejía-Monasterio, C., Muratore-Ginanneschi, P.: Optimal protocols and optimal transport in stochastic thermodynamics. Phys. Rev. Lett. 106(25), 250601 (2011). 10.1103/PhysRevLett.106.250601. arXiv: 1012.2037 21770620
24. Sivak DA Crooks GE Thermodynamic metrics and optimal paths Phys. Rev. Lett. 2012 108 19 190602 23003019
Sivak, D.A., Crooks, G.E.: Thermodynamic metrics and optimal paths. Phys. Rev. Lett. 108(19), 190602 (2012). 10.1103/PhysRevLett.108.190602. arXiv:1201.4166 [cond-mat.stat-mech]23003019
25. Rotskoff GM Crooks GE Vanden-Eijnden E Geometric approach to optimal nonequilibrium control: minimizing dissipation in nanomagnetic spin systems Phys. Rev. E 2017 95 1 012148 28208424
Rotskoff, G.M., Crooks, G.E., Vanden-Eijnden, E.: Geometric approach to optimal nonequilibrium control: minimizing dissipation in nanomagnetic spin systems. Phys. Rev. E 95(1), 012148 (2017). 10.1103/PhysRevE.95.012148. arXiv:1607.07425 [cond-mat.stat-mech]28208424
26. Baldovin M Guéry-Odelin D Trizac E Control of active Brownian particles: an exact solution Phys. Rev. Lett. 2023 131 118302 37774311
Baldovin, M., Guéry-Odelin, D., Trizac, E.: Control of active Brownian particles: an exact solution. Phys. Rev. Lett. 131, 118302 (2023). 10.1103/PhysRevLett.131.11830237774311
27. Chennakesavalu S Rotskoff GM Unified, geometric framework for nonequilibrium protocol optimization Phys. Rev. Lett. 2023 10 130
Chennakesavalu, S., Rotskoff, G.M.: Unified, geometric framework for nonequilibrium protocol optimization. Phys. Rev. Lett. 10, 130 (2023). 10.1103/physrevlett.130.107101
28. Aurell E Gawȩdzki K Mejía-Monasterio C Mohayaee R Muratore-Ginanneschi P Refined second law of thermodynamics for fast random processes J. Stat. Phys. 2012 147 3 487 505
Aurell, E., Gawȩdzki, K., Mejía-Monasterio, C., Mohayaee, R., Muratore-Ginanneschi, P.: Refined second law of thermodynamics for fast random processes. J. Stat. Phys. 147(3), 487–505 (2012). 10.1007/s10955-012-0478-x. arXiv:1201.3207 [cond-mat.stat-mech]
29. Gawȩdzki, K.: Fluctuation Relations in Stochastic Thermodynamics. Lecture notes, arXiv.org:1308.1518 (2013)
30. Villani C Optimal Transport: Old and New 2009 Berlin Springer 973
Villani, C.: Optimal Transport: Old and New. Grundlehren der mathematischen Wissenschaften, vol. 338, p. 973. Springer, Berlin (2009)
31. Brenier Y Frisch U Hénon M Loeper G Matarrese S Mohayaee R Sobolevskiǐ A Reconstruction of the early Universe as a convex optimization problem Monthly Not. R. Astronom. Soc. 2003 346 2 501 524
Brenier, Y., Frisch, U., Hénon, M., Loeper, G., Matarrese, S., Mohayaee, R., Sobolevskiǐ, A.: Reconstruction of the early Universe as a convex optimization problem. Monthly Not. R. Astronom. Soc. 346(2), 501–524 (2003). 10.1046/j.1365-2966.2003.07106.x. arXiv:astro-ph/0304214 [astro-ph]
32. Muratore-Ginanneschi P Mejía-Monasterio C Peliti L Heat release by controlled continuous-time Markov jump processes J. Stat. Phys. 2013 150 1 181 203
Muratore-Ginanneschi, P., Mejía-Monasterio, C., Peliti, L.: Heat release by controlled continuous-time Markov jump processes. J. Stat. Phys. 150(1), 181–203 (2013). 10.1007/s10955-012-0676-6. arXiv:1203.4062 [cond-mat.stat-mech]
33. Muratore-Ginanneschi P (2014) On extremals of the entropy production by “Langevin–Kramers” dynamics J. Stat. Mech. 2014 5 05013
Muratore-Ginanneschi, P.: (2014) On extremals of the entropy production by “Langevin–Kramers’’ dynamics. J. Stat. Mech. 5, 05013 (2014). 10.1088/1742-5468/2014/05/p05013. arXiv:1401.3394 [cond-mat.stat-mech]
34. Muratore-Ginanneschi P Schwieger K How nanomechanical systems can minimize dissipation Phys. Rev. E 2014 90 6 060102
Muratore-Ginanneschi, P., Schwieger, K.: How nanomechanical systems can minimize dissipation. Phys. Rev. E 90(6), 060102 (2014). 10.1103/PhysRevE.90.060102. arXiv:1408.5298 [cond-mat.stat-mech]
35. Shiraishi N Funo K Saito K Speed limit for classical stochastic processes Phys. Rev. Lett. 2018 121 7 070601 30169075
Shiraishi, N., Funo, K., Saito, K.: Speed limit for classical stochastic processes. Phys. Rev. Lett. 121(7), 070601 (2018). 10.1103/PhysRevLett.121.070601. arXiv: 1802.06554 30169075
36. Remlein B Seifert U Optimality of nonconservative driving for finite-time processes with discrete states Phys. Rev. E 2021 103 5
Remlein, B., Seifert, U.: Optimality of nonconservative driving for finite-time processes with discrete states. Phys. Rev. E 103, 5 (2021). 10.1103/physreve.103.l050105
37. Martínez IA Roldán E Dinis L Petrov D Parrondo JMR Rica RA Brownian carnot engine Nat. Phys. 2016 12 67 70 27330541
Martínez, I.A., Roldán, E., Dinis, L., Petrov, D., Parrondo, J.M.R., Rica, R.A.: Brownian carnot engine. Nat. Phys. 12, 67–70 (2016). 10.1038/nphys3518. arXiv:1412.1282 [cond-mat.stat-mech]27330541
38. Dinis, L., Martínez, I.A., Roldán, E., Parrondo, J.M.R., Rica, R.A.: Thermodynamics at the microscale: from effective heating to the Brownian Carnot engine. J. Stat. Mech. , 2016 (2016). 10.1088/1742-5468/2016/05/054003
39. Guéry-Odelin D Ruschhaupt A Kiely A Torrontegui E Martínez-Garaot S Muga JG Shortcuts to adiabaticity: concepts, methods, and applications Rev. Modern Phys. 2019 91 045001
Guéry-Odelin, D., Ruschhaupt, A., Kiely, A., Torrontegui, E., Martínez-Garaot, S., Muga, J.G.: Shortcuts to adiabaticity: concepts, methods, and applications. Rev. Modern Phys. 91, 045001 (2019). 10.1103/RevModPhys.91.045001
40. Raynal D Guillebon T Guéry-Odelin D Trizac E Lauret J-S Rondin L Shortcuts to equilibrium with a levitated particle in the underdamped regime Phys. Rev. Lett. 2023 131 087101 37683149
Raynal, D., Guillebon, T., Guéry-Odelin, D., Trizac, E., Lauret, J.-S., Rondin, L.: Shortcuts to equilibrium with a levitated particle in the underdamped regime. Phys. Rev. Lett. 131, 087101 (2023). 10.1103/PhysRevLett.131.087101. arXiv: 2303.09542 37683149
41. Plata CA Prados A Trizac E Guéry-Odelin D Taming the time evolution in overdamped systems: shortcuts elaborated from fast-forward and time-reversed protocols Phys. Rev. Lett. 2021 127 190605 34797129
Plata, C.A., Prados, A., Trizac, E., Guéry-Odelin, D.: Taming the time evolution in overdamped systems: shortcuts elaborated from fast-forward and time-reversed protocols. Phys. Rev. Lett. 127, 190605 (2021). 10.1103/PhysRevLett.127.19060534797129
42. Baldovin M Guéry-Odelin D Trizac E Shortcuts to adiabaticity for Lévy processes in harmonic traps Phys. Rev. E 2022 106 054122 36559466
Baldovin, M., Guéry-Odelin, D., Trizac, E.: Shortcuts to adiabaticity for Lévy processes in harmonic traps. Phys. Rev. E 106, 054122 (2022). 10.1103/PhysRevE.106.05412236559466
43. Guéry-Odelin D Jarzynski C Plata CA Prados A Trizac E Driving rapidly while remaining in control: classical shortcuts from Hamiltonian to stochastic dynamics Rep. Progress Phys. 2023 86 3 035902
Guéry-Odelin, D., Jarzynski, C., Plata, C.A., Prados, A., Trizac, E.: Driving rapidly while remaining in control: classical shortcuts from Hamiltonian to stochastic dynamics. Rep. Progress Phys. 86(3), 035902 (2023)
44. López-Suárez M Neri I Gammaitoni L Sub-k bt micro-electromechanical irreversible logic gate Nat. Commun. 2016 7 1 12068 27350333
López-Suárez, M., Neri, I., Gammaitoni, L.: Sub-k bt micro-electromechanical irreversible logic gate. Nat. Commun. 7(1), 12068 (2016). 10.1038/ncomms1206827350333
45. Deshpande A Gopalkrishnan M Ouldridge TE Jones NS Designing the optimal bit: balancing energetic cost, speed and reliability Proc. R. Soc. A 2017 473 2204 20170117 28878557
Deshpande, A., Gopalkrishnan, M., Ouldridge, T.E., Jones, N.S.: Designing the optimal bit: balancing energetic cost, speed and reliability. Proc. R. Soc. A 473(2204), 20170117 (2017). 10.1098/rspa.2017.011728878557
46. Ciampini, M.A., Wenzl, T., Konopik, M., Thalhammer, G., Aspelmeyer, M., Lutz, E., Kiesel, N.: Experimental nonequilibrium memory erasure beyond Landauer’s bound (2021)
47. Lent CS Anderson NG Sagawa T Porod W Ciliberto S Lutz E Orlov AO Hänninen IK Campos-Aguillón CO Celis-Cordova R McConnell MS Szakmany GP Thorpe CC Appleton BT Boechler GP Snider GL Energy Limits in Computation 2019 Cham Springer
Lent, C.S., Anderson, N.G., Sagawa, T., Porod, W., Ciliberto, S., Lutz, E., Orlov, A.O., Hänninen, I.K., Campos-Aguillón, C.O., Celis-Cordova, R., McConnell, M.S., Szakmany, G.P., Thorpe, C.C., Appleton, B.T., Boechler, G.P., Snider, G.L.: Energy Limits in Computation. Springer, Cham (2019). 10.1007/978-3-319-93458-7
48. Ray KJ Boyd AB Wimsatt GW Crutchfield JP Non-markovian momentum computing: Thermodynamically efficient and computation universal Phys. Rev. Res. 2021 3 023164
Ray, K.J., Boyd, A.B., Wimsatt, G.W., Crutchfield, J.P.: Non-markovian momentum computing: Thermodynamically efficient and computation universal. Phys. Rev. Res. 3, 023164 (2021). 10.1103/PhysRevResearch.3.023164
49. Dago S Pereda J Barros N Ciliberto S Bellon L Information and thermodynamics: fast and precise approach to Landauer’s bound in an underdamped micromechanical oscillator Phys. Rev. Lett. 2021 126 170601 33988419
Dago, S., Pereda, J., Barros, N., Ciliberto, S., Bellon, L.: Information and thermodynamics: fast and precise approach to Landauer’s bound in an underdamped micromechanical oscillator. Phys. Rev. Lett. 126, 170601 (2021). 10.1103/PhysRevLett.126.17060133988419
50. Proesmans K Ehrich J Bechhoefer J Finite-time Landauer principle Phys. Rev. Lett. 2020 125 10 100602 32955336
Proesmans, K., Ehrich, J., Bechhoefer, J.: Finite-time Landauer principle. Phys. Rev. Lett. 125(10), 100602 (2020). 10.1103/physrevlett.125.100602. arXiv: 2006.03242 32955336
51. Proesmans K Ehrich J Bechhoefer J Optimal finite-time bit erasure under full control Phys. Rev. E 2020 3 102
Proesmans, K., Ehrich, J., Bechhoefer, J.: Optimal finite-time bit erasure under full control. Phys. Rev. E 3, 102 (2020). 10.1103/physreve.102.032105
52. Zhen YZ Egloff D Modi K Dahlsten O Universal bound on energy cost of bit reset in finite time Phys. Rev. Lett. 2021 19 127
Zhen, Y.Z., Egloff, D., Modi, K., Dahlsten, O.: Universal bound on energy cost of bit reset in finite time. Phys. Rev. Lett. 19, 127 (2021). 10.1103/physrevlett.127.190602
53. Gonzalez-Ballestero C Aspelmeyer M Novotny L Quidant R Romero-Isart O Levitodynamics: Levitation and control of microscopic objects in vacuum Science 2021 6564 374
Gonzalez-Ballestero, C., Aspelmeyer, M., Novotny, L., Quidant, R., Romero-Isart, O.: Levitodynamics: Levitation and control of microscopic objects in vacuum. Science 6564, 374 (2021). 10.1126/science.abg3027
54. Dago S Bellon L Dynamics of information erasure and extension of Landauer’s bound to fast processes Phys. Rev. Lett. 2022 128 070604 35244423
Dago, S., Bellon, L.: Dynamics of information erasure and extension of Landauer’s bound to fast processes. Phys. Rev. Lett. 128, 070604 (2022). 10.1103/PhysRevLett.128.07060435244423
55. Dago S Pereda J Ciliberto S Bellon L Virtual double-well potential for an underdamped oscillator created by a feedback loop J. Stat. Mech. 2023 5 202
Dago, S., Pereda, J., Ciliberto, S., Bellon, L.: Virtual double-well potential for an underdamped oscillator created by a feedback loop. J. Stat. Mech. 5, 202 (2023)
56. Mikami T Monge’s problem with a quadratic cost by the zero-noise limit of h -path processes Probab. Theory Relat. Fields 2004 129 2 245 260
Mikami, T.: Monge’s problem with a quadratic cost by the zero-noise limit of h -path processes. Probab. Theory Relat. Fields 129(2), 245–260 (2004). 10.1007/s00440-004-0340-4
57. Muratore-Ginanneschi P On the use of stochastic differential geometry for non-equilibrium thermodynamics modeling and control J. Phys. A 2013 46 27 275002
Muratore-Ginanneschi, P.: On the use of stochastic differential geometry for non-equilibrium thermodynamics modeling and control. J. Phys. A 46(27), 275002 (2013). 10.1088/1751-8113/46/27/275002
58. Chen Y Georgiou TT Pavon M On the relation between optimal transport and Schrödinger bridges: a stochastic control viewpoint J. Optim. Theory Appl. 2016 169 2 671 691
Chen, Y., Georgiou, T.T., Pavon, M.: On the relation between optimal transport and Schrödinger bridges: a stochastic control viewpoint. J. Optim. Theory Appl. 169(2), 671–691 (2016). 10.1007/s10957-015-0803-z
59. Cuccoli A Fubini A Tognetti V Vaia R Quantum thermodynamics of systems with anomalous dissipative coupling Phys. Rev. E 2001 6 64
Cuccoli, A., Fubini, A., Tognetti, V., Vaia, R.: Quantum thermodynamics of systems with anomalous dissipative coupling. Phys. Rev. E 6, 64 (2001). 10.1103/physreve.64.066124
60. Ankerhold J Pollak E Dissipation can enhance quantum effects Phys. Rev. E 2007 4 75
Ankerhold, J., Pollak, E.: Dissipation can enhance quantum effects. Phys. Rev. E 4, 75 (2007). 10.1103/physreve.75.041103
61. Bonilla LL Carrillo JA Soler JS Asymptotic behavior of an initial-boundary value problem for the Vlasov–Poisson–Fokker–Planck System SIAM J. Appl. Math. 1997 57 5 1343 1372
Bonilla, L.L., Carrillo, J.A., Soler, J.S.: Asymptotic behavior of an initial-boundary value problem for the Vlasov–Poisson–Fokker–Planck System. SIAM J. Appl. Math. 57(5), 1343–1372 (1997). 10.1137/S0036139995291544
62. Bhatia, R., Elsner, L.: Positive linear maps and the lyapunov equation. In: Linear Operators and Matrices, pp. 107–120. Birkhäuser, Basel (2002). 10.1007/978-3-0348-8181-4_9
63. Verhulst F Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics 2005 New York Springer
Verhulst, F.: Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics. Texts in Applied Mathematics, Springer, New York (2005). 10.1007/0-387-28313-7
64. Amit, D.J., Martin-Mayor, V.: Field Theory, the Renormalization Group, and Critical Phenomena, 3rd edn. In: International series in pure and applied physics, p. 568. World Scientific Publishing, Singapore (2005). 10.1142/5715
65. Wycoff D Bálazs NL Multiple time scales analysis for the Kramers–Chandrasekhar equation Physica A 1987 146 1–2 175 200
Wycoff, D., Bálazs, N.L.: Multiple time scales analysis for the Kramers–Chandrasekhar equation. Physica A 146(1–2), 175–200 (1987). 10.1016/0378-4371(87)90227-5
66. Wycoff D Bálazs NL Multiple time scales analysis for the Kramers–Chandrasekhar equation with a weak magnetic field Physica A 1987 146 1–2 201 218
Wycoff, D., Bálazs, N.L.: Multiple time scales analysis for the Kramers–Chandrasekhar equation with a weak magnetic field. Physica A 146(1–2), 201–218 (1987). 10.1016/0378-4371(87)90228-7
67. Gentile G Quasi-periodic motions in dynamical systems. Review of a Renormalisation group approach J. Math. Phys. 2010 51 015207
Gentile, G.: Quasi-periodic motions in dynamical systems. Review of a Renormalisation group approach. J. Math. Phys. 51, 015207 (2010). 10.1063/1.3271653. arXiv:0910.0755 [math.DS]
68. Chiarini, A., Conforti, G., Greco, G., Ren, Z.: Entropic turnpike estimates for the kinetic Schrödinger problem. eprint arXiv:2108.09161 [math.PR] (2021)
69. Muratore-Ginanneschi P Peliti L Classical uncertainty relations and entropy production in non-equilibrium statistical mechanics J. Stat. Mech. 2023 8 083202083202
Muratore-Ginanneschi, P., Peliti, L.: Classical uncertainty relations and entropy production in non-equilibrium statistical mechanics. J. Stat. Mech. 8, 083202083202 (2023). 10.1088/1742-5468/ace3b3
70. Lapolla A Godec A Faster uphill relaxation in thermodynamically equidistant temperature quenches Phys. Rev. Lett. 2020 11 125
Lapolla, A., Godec, A.: Faster uphill relaxation in thermodynamically equidistant temperature quenches. Phys. Rev. Lett. 11, 125 (2020). 10.1103/physrevlett.125.110602
71. Ibáñez M Dieball C Lasanta A Godec A Rica RA Heating and cooling are fundamentally asymmetric and evolve along distinct pathways Nat. Phys. 2024 20 1 135 141
Ibáñez, M., Dieball, C., Lasanta, A., Godec, A., Rica, R.A.: Heating and cooling are fundamentally asymmetric and evolve along distinct pathways. Nat. Phys. 20(1), 135–141 (2024). 10.1038/s41567-023-02269-z. arXiv:2302.09061
72. Pavliotis GA Stuart AM Multiscale Methods: Averaging and Homogenization 2008 New York Springer 307
Pavliotis, G.A., Stuart, A.M.: Multiscale Methods: Averaging and Homogenization. Texts in Applied Mathematics, p. 307. Springer, New York (2008)
73. Schmiedl T Seifert U Efficiency of molecular motors at maximum power EPL (Europhys. Lett.) 2008 83 3 30005
Schmiedl, T., Seifert, U.: Efficiency of molecular motors at maximum power. EPL (Europhys. Lett.) 83(3), 30005 (2008). 10.1209/0295-5075/83/30005. arXiv:0801.3743 [cond-mat.stat-mech]
74. Muratore-Ginanneschi P Schwieger K Efficient protocols for Stirling heat engines at the micro-scale EPL (Europhys. Lett.) 2015 112 20002
Muratore-Ginanneschi, P., Schwieger, K.: Efficient protocols for Stirling heat engines at the micro-scale. EPL (Europhys. Lett.) 112, 20002 (2015). 10.1209/0295-5075/112/20002. arXiv:1503.05788 [cond-mat.stat-mech]
75. Dechant A Kiesel N Lutz E Underdamped stochastic heat engine at maximum efficiency EPL (Europhys. Lett.) 2017 119 50003
Dechant, A., Kiesel, N., Lutz, E.: Underdamped stochastic heat engine at maximum efficiency. EPL (Europhys. Lett.) 119, 50003 (2017). 10.1209/0295-5075/119/50003
76. Pavliotis GA Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations 2014 New York Springer
Pavliotis, G.A.: Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations. Springer, New York (2014). 10.1007/978-1-4939-1323-7
77. Caldeira AO Leggett AJ Quantum tunnelling in a dissipative system Ann. Phys. 1983 149 2 374 456
Caldeira, A.O., Leggett, A.J.: Quantum tunnelling in a dissipative system. Ann. Phys. 149(2), 374–456 (1983). 10.1016/0003-4916(83)90202-6
78. Maile D Andergassen S Rastelli G Effects of a dissipative coupling to the momentum of a particle in a double well potential Phys. Rev. Res. 2020 1 2
Maile, D., Andergassen, S., Rastelli, G.: Effects of a dissipative coupling to the momentum of a particle in a double well potential. Phys. Rev. Res. 1, 2 (2020). 10.1103/physrevresearch.2.013226
79. Zwanzig R Nonlinear generalized Langevin equations J. Stat. Phys. 1973 9 3 215 220
Zwanzig, R.: Nonlinear generalized Langevin equations. J. Stat. Phys. 9(3), 215–220 (1973). 10.1007/bf01008729
80. Conforti G Ripani L Around the entropic Talagrand inequality Bernoulli 2020 26 1431 1452
Conforti, G., Ripani, L.: Around the entropic Talagrand inequality. Bernoulli 26, 1431–1452 (2020). 10.48550/ARXIV.1809.02062
81. Theodorou, E.A., Todorov, E.: Relative entropy and free energy dualities: Connections to Path Integral and Kullback–Leibler control. In: Annual Conference on Decision and Control (CDC), 2012 IEEE 51st, pp. 1466–1473 (2012). 10.1109/CDC.2012.6426381
82. Boscain U Sigalotti M Sugny D Introduction to the pontryagin maximum principle for quantum optimal control PRX Quantum 2021 3 2
Boscain, U., Sigalotti, M., Sugny, D.: Introduction to the pontryagin maximum principle for quantum optimal control. PRX Quantum 3, 2 (2021). 10.1103/prxquantum.2.030203
83. Rey-Bellet L Ergodic Properties of Markov Processes. In: Quantum Open Systems II. The Markovian Approach 2006 Berlin Springer 1 39
Rey-Bellet, L.: Ergodic Properties of Markov Processes. In: Quantum Open Systems II. The Markovian Approach. Lecture Notes in Mathematics, pp. 1–39. Springer, Berlin (2006)
84. Caluya KF Halder A Wasserstein proximal algorithms for the Schrödinger bridge problem: density control with nonlinear drift IEEE Trans. Automatic Control 2022 3 67
Caluya, K.F., Halder, A.: Wasserstein proximal algorithms for the Schrödinger bridge problem: density control with nonlinear drift. IEEE Trans. Automatic Control 3, 67 (2022). 10.1109/TAC.2021.3060704
85. Talagrand M Transportation cost for Gaussian and other product measures Geometric Funct. Anal. 1996 6 3 587 600
Talagrand, M.: Transportation cost for Gaussian and other product measures. Geometric Funct. Anal. 6(3), 587–600 (1996). 10.1007/bf02249265
86. Otto F Villani C Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality J. Funct. Anal. 2000 173 2 361 400
Otto, F., Villani, C.: Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality. J. Funct. Anal. 173(2), 361–400 (2000). 10.1006/jfan.1999.3557
87. Léonard C From the Schrödinger problem to the Monge–Kantorovich problem J. Funct. Anal. 2012 262 4 1879 1920
Léonard, C.: From the Schrödinger problem to the Monge–Kantorovich problem. J. Funct. Anal. 262(4), 1879–1920 (2012). 10.1016/j.jfa.2011.11.026. arXiv:1011.2564 [math.OC]
88. Nelson E Dynamical Theories of Brownian Motion 2001 2 Princeton Princeton University Press 148
Nelson, E.: Dynamical Theories of Brownian Motion, 2nd edn., p. 148. Princeton University Press, Princeton (2001). 10.2307/j.ctv15r57jg
89. Benamou J-D Brenier Y A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem Numerische Mathematik 2000 84 3 375 393
Benamou, J.-D., Brenier, Y.: A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numerische Mathematik 84(3), 375–393 (2000). 10.1007/s002110050002
90. Dechant A Sasa SI Entropic bounds on currents in Langevin systems Phys. Rev. 2018 6 97
Dechant, A., Sasa, S.I.: Entropic bounds on currents in Langevin systems. Phys. Rev. 6, 97 (2018). 10.1103/physreve.97.062101
91. Gawedzki, K.: Improved 2nd Law of Stochastic Thermodynamics for underdamped Langevin process. Note written for Salambô Dago in July 2020, and communicated by the author to Erik Aurell (2021)
92. Mikami T Thieullen M Duality theorem for stochastic optimal control problem Stoch. Process. Appl. 2006 116 1815 1835
Mikami, T., Thieullen, M.: Duality theorem for stochastic optimal control problem. Stoch. Process. Appl. 116, 1815–1835 (2006)
93. Gentil I Léonard C Ripani L About the analogy between optimal transport and minimal entropy Ann. Facult. Sci. Toulouse 2017 26 3 569 600
Gentil, I., Léonard, C., Ripani, L.: About the analogy between optimal transport and minimal entropy. Ann. Facult. Sci. Toulouse 26(3), 569–600 (2017). 10.5802/afst.1546
94. Serrin, J.: Mathematical principles of classical fluid mechanics. In: Fluid Dynamics I / Strömungsmechanik I. Encyclopedia of Physics / Handbuch der Physik, vol. 3 / 8 / 1, pp. 125–263. Springer, Berlin, Heidelberg (1959). 10.1007/978-3-642-45914-6_2
95. Seliger RL Whitham GB Variational principles in continuum mechanics Proc. R. Soc. A 1968 305 1480 1 25
Seliger, R.L., Whitham, G.B.: Variational principles in continuum mechanics. Proc. R. Soc. A 305(1480), 1–25 (1968). 10.1098/rspa.1968.0103
96. Bismut JM Kohlmann M Vogel W An introduction to duality in random mechanics Stochastic Control Theory and Stochastic Differential Systems 1979 42 60
Bismut, J.M.: An introduction to duality in random mechanics. In: Kohlmann, M., Vogel, W. (eds.) Stochastic Control Theory and Stochastic Differential Systems. Lecture Notes in Control and Information Sciences, pp. 42–60. (1979). 10.1007/BFb0009375
97. Bechhoefer J Control Theory for Physicists 2021 Cambridge Cambridge University Press
Bechhoefer, J.: Control Theory for Physicists. Cambridge University Press, Cambridge (2021). 10.1017/9780511734809
98. Liberzon D Calculus of Variations and Optimal Control Theory. A Concise Introduction 2012 Princeton Princeton University Press
Liberzon, D.: Calculus of Variations and Optimal Control Theory. A Concise Introduction. Princeton University Press, Princeton (2012)
99. Chen Y Georgiou TT Pavon M Fast cooling for a system of stochastic oscillators J. Math. Phys. 2015 11 56
Chen, Y., Georgiou, T.T., Pavon, M.: Fast cooling for a system of stochastic oscillators. J. Math. Phys. 11, 56 (2015). 10.1063/1.4935435
100. Caluya KF Halder A Gradient flow algorithms for density propagation in stochastic systems IEEE Trans. Autom. Control 2020 65 10 3991 4004
Caluya, K.F., Halder, A.: Gradient flow algorithms for density propagation in stochastic systems. IEEE Trans. Autom. Control 65(10), 3991–4004 (2020). 10.1109/tac.2019.2951348
101. Hayes M Kaper TJ Kopell N Ono K On the application of geometric singular perturbation theory to some classical two point boundary value problems Int. J. Bifurcat. Chaos 1998 08 02 189 209
Hayes, M., Kaper, T.J., Kopell, N., Ono, K.: On the application of geometric singular perturbation theory to some classical two point boundary value problems. Int. J. Bifurcat. Chaos 08(02), 189–209 (1998). 10.1142/s0218127498000140
102. Lehec J Representation formula for the entropy and functional inequalities Annales de l’Institut Henri Poincaré, Probabilités et Statistiques 2013 3 49
Lehec, J.: Representation formula for the entropy and functional inequalities. Annales de l’Institut Henri Poincaré, Probabilités et Statistiques 3, 49 (2013). 10.1214/11-aihp464
103. Bobylev AV Instabilities in the Chapman–Enskog expansion and hyperbolic Burnett equations J. Stat. Phys. 2006 124 2–4 371 399
Bobylev, A.V.: Instabilities in the Chapman–Enskog expansion and hyperbolic Burnett equations. J. Stat. Phys. 124(2–4), 371–399 (2006). 10.1007/s10955-005-8087-6
104. Rackauckas C Nie O Differentialequations.jl-: a performant and feature-rich ecosystem for solving differential equations in julia J. Open Res. Softw. 2017 1 15
Rackauckas, C., Nie, O.: Differentialequations.jl-: a performant and feature-rich ecosystem for solving differential equations in julia. J. Open Res. Softw. 1, 15 (2017)
105. Bérut A Arakelyan A Petrosyan A Ciliberto S Dillenschneider R Lutz E Experimental verification of Landauer’s principle linking information and thermodynamics Nature 2012 483 187 189 22398556
Bérut, A., Arakelyan, A., Petrosyan, A., Ciliberto, S., Dillenschneider, R., Lutz, E.: Experimental verification of Landauer’s principle linking information and thermodynamics. Nature 483, 187–189 (2012). 10.1038/nature1087222398556
106. Martínez IA Petrosyan A Guéry-Odelin D Trizac E Ciliberto S Engineered swift equilibration of a Brownian particle Nat. Phys. 2016 12 843 846 27610190
Martínez, I.A., Petrosyan, A., Guéry-Odelin, D., Trizac, E., Ciliberto, S.: Engineered swift equilibration of a Brownian particle. Nat. Phys. 12, 843–846 (2016). 10.1038/nphys375827610190
107. Sanders, J., Baldovin, M., Muratore-Ginanneschi, P.: Minimal-work protocols for inertial particles in non-harmonic traps. Eprint arXiv:2407.15678 (2024) 10.48550/ARXIV.2407.15678
108. Sipos O Nagy K Di Leonardo R Galajda P Hydrodynamic trapping of swimming bacteria by convex walls Phys. Rev. Lett. 2015 114 258104 26197146
Sipos, O., Nagy, K., Di Leonardo, R., Galajda, P.: Hydrodynamic trapping of swimming bacteria by convex walls. Phys. Rev. Lett. 114, 258104 (2015). 10.1103/PhysRevLett.114.25810426197146
109. Peng C Turiv T Guo Y Wei Q-H Lavrentovich OD Command of active matter by topological defects and patterns Science 2016 354 6314 882 885 27856907
Peng, C., Turiv, T., Guo, Y., Wei, Q.-H., Lavrentovich, O.D.: Command of active matter by topological defects and patterns. Science 354(6314), 882–885 (2016). 10.1126/science.aah693627856907
110. Cavagna A Giardina I Gucciardino MA Iacomelli G Lombardi M Melillo S Monacchia G Parisi L Peirce MJ Spaccapelo R Characterization of lab-based swarms of anopheles Gambiae mosquitoes using 3d-video tracking Sci. Rep. 2023 13 1 8745 37253765
Cavagna, A., Giardina, I., Gucciardino, M.A., Iacomelli, G., Lombardi, M., Melillo, S., Monacchia, G., Parisi, L., Peirce, M.J., Spaccapelo, R.: Characterization of lab-based swarms of anopheles Gambiae mosquitoes using 3d-video tracking. Sci. Rep. 13(1), 8745 (2023)37253765
111. Pellicciotta N Paoluzzi M Buonomo D Frangipane G Angelani L Di Leonardo R Colloidal transport by light induced gradients of active pressure. Nat. Commun. 2023 1 14
Pellicciotta, N., Paoluzzi, M., Buonomo, D., Frangipane, G., Angelani, L., Di Leonardo, R.: Colloidal transport by light induced gradients of active pressure. Nat. Commun. 1, 14 (2023). 10.1038/s41467-023-39974-5
112. Shankar S Raju V Mahadevan L Optimal transport and control of active drops Proc. Nat. Acad. Sci. 2022 35 119
Shankar, S., Raju, V., Mahadevan, L.: Optimal transport and control of active drops. Proc. Nat. Acad. Sci. 35, 119 (2022). 10.1073/pnas.2121985119
113. Davis LK Proesmans K Fodor E Active matter under control: insights from response theory Phys. Rev. X 2024 14 011012
Davis, L.K., Proesmans, K., Fodor, E.: Active matter under control: insights from response theory. Phys. Rev. X 14, 011012 (2024). 10.1103/PhysRevX.14.011012
114. Frim, A.G., DeWeese, M.R.: Shortcut engineering of active matter: run-and-tumble particles (2023)
115. Manacorda A Puglisi A Lattice model to derive the fluctuating hydrodynamics of active particles with inertia Phys. Rev. Lett. 2017 119 208003 29219378
Manacorda, A., Puglisi, A.: Lattice model to derive the fluctuating hydrodynamics of active particles with inertia. Phys. Rev. Lett. 119, 208003 (2017). 10.1103/PhysRevLett.119.20800329219378
116. Scholz C Jahanshahi S Ldov A Löwen H Inertial delay of self-propelled particles Nat. Commun. 2018 1 9
Scholz, C., Jahanshahi, S., Ldov, A., Löwen, H.: Inertial delay of self-propelled particles. Nat. Commun. 1, 9 (2018). 10.1038/s41467-018-07596-x
117. Löwen H Inertial effects of self-propelled particles: from active Brownian to active Langevin motion J. Chem. Phys. 2019 152 4 040901
Löwen, H.: Inertial effects of self-propelled particles: from active Brownian to active Langevin motion. J. Chem. Phys. 152(4), 040901 (2019)
118. Attanasi A Cavagna A Del Castello L Giardina I Grigera TS Jelić A Melillo S Parisi L Pohl O Shen E Viale M Information transfer and behavioural inertia in starling flocks Nat. Phys. 2014 10 9 691 696
Attanasi, A., Cavagna, A., Del Castello, L., Giardina, I., Grigera, T.S., Jelić, A., Melillo, S., Parisi, L., Pohl, O., Shen, E., Viale, M.: Information transfer and behavioural inertia in starling flocks. Nat. Phys. 10(9), 691–696 (2014). 10.1038/nphys3035
119. Cavagna A Del Castello L Giardina I Grigera T Jelic A Melillo S Mora T Parisi L Silvestri E Viale M Walczak AM Flocking and turning: a new model for self-organized collective motion J. Stat. Phys. 2014 158 3 601 627
Cavagna, A., Del Castello, L., Giardina, I., Grigera, T., Jelic, A., Melillo, S., Mora, T., Parisi, L., Silvestri, E., Viale, M., Walczak, A.M.: Flocking and turning: a new model for self-organized collective motion. J. Stat. Phys. 158(3), 601–627 (2014). 10.1007/s10955-014-1119-3
120. Cavagna A Conti D Creato C Del Castello L Giardina I Grigera TS Melillo S Parisi L Viale M Dynamic scaling in natural swarms Nat. Phys. 2017 13 9 914 918
Cavagna, A., Conti, D., Creato, C., Del Castello, L., Giardina, I., Grigera, T.S., Melillo, S., Parisi, L., Viale, M.: Dynamic scaling in natural swarms. Nat. Phys. 13(9), 914–918 (2017). 10.1038/nphys4153
