# Fluctuations of dynamical observables in linear diffusions with time delay:
a Riccati-based approach

M.L. Rosinberg<EMAIL_ADDRESS>Laboratoire de Physique Théorique de la
Matière Condensée, CNRS-UMR 7600, Sorbonne Université, 4 place Jussieu, 75252
Paris Cedex 05, France  G. Tarjus  Laboratoire de Physique Théorique de la
Matière Condensée, CNRS-UMR 7600, Sorbonne Université, 4 place Jussieu, 75252
Paris Cedex 05, France  T. Munakata  Department of Applied Mathematics and
Physics, Graduate School of Informatics, Kyoto University, Kyoto 606-8501,
Japan

###### Abstract

Our current understanding of fluctuations of dynamical (time-integrated)
observables in non-Markovian processes is still very limited. A major obstacle
is the lack of an appropriate theoretical framework to evaluate the associated
large deviation functions. In this paper we bypass this difficulty in the case
of linear diffusions with time delay by using a Markovian embedding procedure
that introduces an infinite set of coupled differential equations. We then
show that the generating functions of current-type observables can be computed
at arbitrary finite time by solving matrix Riccati differential equations
(RDEs) somewhat similar to those encountered in optimal control and filtering
problems. By exploring in detail the properties of these RDEs and of the
corresponding continuous-time algebraic Riccati equations (CAREs), we identify
the generic fixed point towards which the solutions converge in the long-time
limit. This allows us to derive the explicit expressions of the scaled
cumulant generating function (SCGF), of the pre-exponential factors, and of
the effective (or driven) process that describes how fluctuations are created
dynamically. Finally, we describe the special behavior occurring at the limits
of the domain of existence of the SCGF, in connection with fluctuation
relations for the heat and the entropy production.

###### Contents

  1. I Introduction
  2. II Model and observables
    1. II.1 Langevin equation and linear chain trick
    2. II.2 Dynamical observables
    3. II.3 Stability of the non-equilibrium steady state
  3. III Matrix Riccati equations and their properties
    1. III.1 Basic equations
    2. III.2 Properties of the Riccati differential equations (RDEs)
      1. III.2.1 Global existence of the solutions
      2. III.2.2 Solutions of the RDEs and asymptotic behavior
    3. III.3 Fixed points of the Riccati flows
  4. IV Time evolution of the moment generating functions
    1. IV.1 Finite-time divergences
    2. IV.2 Long-time behavior
      1. IV.2.1 Domain of existence of the scaled-cumulant generating function (SCGF)
      2. IV.2.2 Expressions of the SCGF and of the pre-exponential prefactors
      3. IV.2.3 Relation with the spectral problem for the tilted generators
      4. IV.2.4 Effective process
      5. IV.2.5 Role of (temporal) boundary terms
    3. IV.3 Behavior of the SCGF at the limits of its interval of definition
      1. IV.3.1 Attraction to a non-extremal solution of the CARE
      2. IV.3.2 Oscillatory behavior
  5. V The special case $\lambda=1$
    1. V.1 Preliminaries
    2. V.2 Integral fluctuation theorem (IFT) for the heat
    3. V.3 Asymptotic fluctuation relation for the entropy production
      1. V.3.1 Acausal dynamics and Jacobian
      2. V.3.2 Calculation of the SCGF
  6. VI Conclusion
  7. A Tilted generators
  8. B Riccati differential equations
    1. B.0.1 Solution of Eqs. (37) and (39)
    2. B.0.2 Explicit expressions of the matrices $K_{\rm o}(\lambda)$
    3. B.0.3 Derivation of Eq. (49)
  9. C Spectral properties of the Hamiltonian matrices
    1. C.0.1 Characteristic polynomial
    2. C.0.2 Eigenvectors
  10. D Derivation of Eq. (77)
  11. E Solutions of the CAREs (86) for $\lambda=0$
  12. F Derivation of the expression of the SCGF $\mu_{\sigma}(\lambda)$ for $\lambda=1$

##  I Introduction

A common problem in nonequilibrium statistical mechanics consists in
estimating the statistics of a time-integrated observable, such as the heat or
the entropy production for a system in contact with a heat reservoir and
driven by an external force. Under quite general conditions, the probability
density of such an observable obeys a large deviation principle in the limit
of large integration times and the fluctuations are then characterized by the
so-called rate or large deviation function (LDF). In the case of Markovian
stochastic dynamics, a powerful mathematical framework has been developed to
compute this quantity and the closely related scaled cumulant generating
function (SCGF) (for reviews, see, e.g., Refs. T2009 ; J2020 ).

However, many stochastic systems in biology, physics and technology are
affected by memory effects and long-range temporal correlations. These occur
in particular in the presence of feedback loops when the time lag between the
signal detection and the system response (or the control operation) makes the
dynamics inherently non-Markovian AM2008 ; A2010 ; JPSS2010 ; M2015 ; B2021 .
The theoretical description of the fluctuations of a dynamical observable then
becomes problematic. The major difficulty is that the time evolution of the
generating function can no longer be formulated as a linear partial
differential equation T2018 . Accordingly, the SCGF can no longer be obtained
as the largest eigenvalue of the corresponding “tilted” generator T2018 ;
CT2015 .

The main goal of this paper is to bypass this obstacle in the case of linear
diffusions with time delay. In this respect, the paper is a sequel of previous
work in which we studied dynamical fluctuations for an underdamped particle
trapped in a harmonic potential and submitted to a position-dependent, time-
delayed feedback force RTM2017 . Such a model describes the motion of
feedback-cooled mechanical resonators, e.g., a microcantilever in the vicinity
of its fundamental mode resonance P2007 ; M2012 ; KK2016 . The most
significant outcome of Ref. RTM2017 , recently tested experimentally DRLK2022
, was that the delay strongly affects the regime of large deviations and that
the fluctuations of heat, work, and entropy production in the nonequilibrium
steady state are quite different (whereas their expectation values are
identical). This feature cannot be rationalized from the sole knowledge of the
LDF: one needs to determine the complete asymptotic behavior of the
probability distributions and generating functions, including the pre-
exponential factors. This nontrivial task (which could not be achieved in Ref.
RTM2017 ) is fulfilled in the present paper.

To this end, and to make the problem mathematically tractable, we use a
procedure known in the literature as the “linear chain trick”, which consists
in replacing the discrete delay by the larger class of gamma-distributed
delays McD1978 . It is indeed a well-known fact in the theory of delay
integro-differential equations that an equivalent system of ordinary
differential equations is obtained whenever the delay is gamma-distributed
S2011 . This procedure, which requires one to introduce auxiliary variables,
is widely used in the context of biological modeling, population dynamics or
evolutionary systems C1979 ; PW1986 ; MD1986 ; McD1989 ; HK2019 . The discrete
delay is recovered when the number of auxiliary variables goes to infinity
LK2019 . Besides the fact that a distributed delay is often more likely to
capture reality than a discrete one (which justifies studying the properties
of such a system per se LHK2021 ), the bonus is that the dynamics of the
augmented system is Markovian. This allows us to work within the standard
framework of Markov processes.

More generally, we propose a method to calculate the generating functions
beyond the large-deviation regime, i.e., for stochastic trajectories of
arbitrary duration. Owing to the linearity of the dynamics and of the current-
type observables, the calculation boils down to solving continuous matrix
Riccati differential equations (RDEs) similar to those encountered in optimal
control and filtering problems. We can thus benefit from the extensive
mathematical literature devoted to the analysis and the numerical solution of
such equations (see, e.g., Refs. AFIJ2003 ; K2010  and references therein),
including for large-size systems BM2004 . However, there is a crucial
difference with the standard situation treated in linear quadratic (LQ)
optimal control which significantly complicates the theoretical analysis. In
particular, the solutions of the RDEs (and in turn the generating functions)
may diverge in a finite time or converge to different fixed points. These
features can be missed when only focusing on the spectral problem for the
dominant eigenvalue of the tilted generator. Our study thus requires a
detailed (and occasionally rather involved) exploration of the properties of
the RDEs and of the corresponding continuous algebraic Riccati equations
(CAREs) in order to anticipate the various possible scenarios. The reward is
that many of the results presented in this work are applicable beyond the
specific case of time-delayed Langevin equations and can be used to compute
the fluctuations of any linear current-type observables in multi-dimensional
linear diffusions111While completing the writing of this paper- which took
much longer than expected- we learned of Refs. DB2023 ; DBT2023  that also use
Riccati differential equations to study dynamical large deviations of linear
diffusions. Similarities and differences with our approach will be discussed
in the text (see in particular Sec. IV.2.5). A preliminary account of the
present work was presented orally at a conference in June 2022 ConfNice .. We
therefore hope that the present analysis will not only provide a better
understanding of the influence of memory effects on dynamical fluctuations,
but more generally will be useful for the application of large deviation
theory to nonequilibrium stochastic systems.

The content of the paper is the following:

In Sec. II.A we present the stochastic underdamped model with a discrete time
delay and we introduce the linear chain trick that makes it possible to
replace the original non-Markovian dynamics by an infinite set of coupled
linear equations without delay. We then define in Sec. II.B the three linear
currents (work, heat, and entropy production) whose fluctuations are commonly
studied in the framework of stochastic thermodynamics in connection with
fluctuation theorems (see, e.g., Refs. S2012 ; PP2021  and references
therein). We also briefly discuss in Sec. II.C how the delay affects the
stability of the nonequilibrium steady state (NESS).

Section III is mainly devoted to the study of the matrix Riccati equations
that play a major role in this work. In Sec. III.1 we first derive the exact
form of the moment generating functions of the fluctuating observables in
terms of real symmetric matrices that are solutions of Riccati differential
equations (RDEs). A specific feature of our treatment is that, for a given
observable, we study together the generating function conditioned on the
initial state (solution of a “backward” PDE) and the generating function
conditioned on the final state (solution of a “forward” PDE). Although this
modus operandi may appear redundant at first sight, it will turn out to be
very useful for analyzing the long-time behavior. In Sec. III.2 we then
investigate in detail the properties of the RDEs, making heavy use of the
concepts and methods available in the mathematical literature. We first
discuss the global existence of the solutions (Sec. III.2.1) and then provide
a closed-form representation of these solutions in terms of the eigenvalues
and eigenvectors of the associated Hamiltonian matrices (Sec. III.2.2). Next,
in Sec. III.3, we use this method to construct explicitly all fixed points of
the RDEs, which are solutions of the corresponding CAREs. The so-called
“maximal” solution is singled out as it is generically the fixed point towards
which the solutions of the RDEs converge asymptotically.

Section IV is the central and longest part of the paper in which we study the
time evolution of the generating functions. To make it more concrete, the
theoretical analysis is illustrated by numerical results obtained for the
gamma-distributed delay. In Sec. IV.1, we first investigate the domain of
existence of the generating functions at finite time. In Sec. IV.2, we then
focus on the long-time behavior. We first determine the domain of existence of
the SCGF (Sec. IV.2.1) and then derive the explicit expressions of the SCGF
and of the pre-exponential factors in terms of the maximal solutions of the
CAREs (Sec. IV.2.2). We also give a representation of the SCGF in terms of an
integral of the spectral density of the process, which shows the connection
between the Riccati-based approach and general results in the mathematical
literature for quadratic observables of stationary Gaussian processes BD1997 ;
BGR1997 ; GRZ1999 ; ZS2023 . In Sec. IV.2.3, we make contact with the standard
spectral problem for the dominant eigenvalue of the tilted generators, which
allows us in Sec. IV.2.4 to characterize the so-called effective or driven
process that describes how fluctuations are created dynamically in the long-
time limit. In Sec. IV.2.5, we then discuss the role of temporal boundary
terms in relation with the recent work of De Buisson and Touchette DBT2023 .
Finally, in Sec. IV.3, we discuss the nontrivial behavior occurring at the
limits of the domain of existence of the SCGF. We show that the solution of
the RDEs may be attracted to a non-maximal solution of the CARE (Sec. IV.3.1)
or may oscillate between two fixed points (Sec. IV.3.2). In both cases the
SCGF displays a positive jump discontinuity.

Finally, in Sec. V, we focus on the special case $\lambda=1$ in connection
with fluctuation relations for the heat and the entropy production. In
particular, we provide the analytical proof of the conjecture relating the
fluctuations of the entropy production at large times to the “Jacobian”
contribution induced by the breaking of causality in the backward process
RTM2017 ; MR2014 ; RMT2015 .

We end the paper in Section VI with a brief summary of the main results.
Several technical calculations and additional details are given in the
Appendices. A summary of the main notations used in this work is also provided
at the end of the paper.

##  II Model and observables

###  II.1 Langevin equation and linear chain trick

As in previous works RTM2017 ; MR2014 ; RMT2015 , we consider a Brownian
particle of mass $m$ trapped in a harmonic potential and immersed in a thermal
environment with viscous damping $\gamma$ and temperature $T$. The dynamical
evolution is governed by the one-dimensional underdamped Langevin equation

$\displaystyle m\dot{v}_{t}=-\gamma v_{t}-kx_{t}+F_{fb}(t)+\sqrt{2\gamma
T}\xi_{t}\ ,$ (1)

where $k$ is the spring constant and $\xi_{t}$ is a zero-mean Gaussian white
noise with unit variance (throughout the paper Boltzmann’s constant is set to
unity). $F_{fb}(t)$ is a feedback control force which is originally taken
proportional to the position of the particle at the time $t-\tau$,

$\displaystyle F_{fb}(t)=k^{\prime}x_{t-\tau}\ ,$ (2)

where $\tau>0$ is the time delay. We generally assume that the feedback is
positive ($k^{\prime}>0$). Eq. (1) accurately describes the motion of the
levitated nanoparticle studied in the experiments of Refs. DRLK2022 ;
DGALK2020 . We stress that the non-Markovian character of the dynamics results
from the feedback and not from the interaction with the environment. By
choosing the inverse angular resonance frequency $\omega_{0}^{-1}=\sqrt{m/k}$
as the unit of time and $x_{c}=\sqrt{2\gamma T}/k$ as the unit of length, the
Langevin equation takes the dimensionless form RMT2015

$\displaystyle\dot{v}_{t}=-\frac{1}{Q_{0}}v_{t}-x_{t}+\frac{g}{Q_{0}}x_{t-\tau}+\xi_{t}\
,$ (3)

where $Q_{0}=\omega_{0}\tau_{0}$ is the quality factor of the oscillator
($\tau_{0}=m/\gamma$ is the viscous relaxation time) and
$g=k^{\prime}/(\gamma\omega_{0})=(k^{\prime}/k)Q_{0}$ is the gain of the
feedback loop. The dynamics is thus fully characterized by the three
dimensionless parameters $Q_{0}$, $g$ and $\tau$.

In order to apply the linear chain trick, we need to smooth the discrete delay
kernel $\delta(t-\tau)$ and replace Eq. (3) by

$\displaystyle\dot{v}_{t}$
$\displaystyle=-\frac{1}{Q_{0}}v_{t}-x_{t}+\frac{g}{Q_{0}}\int_{-\infty}^{t}ds\>g_{n}(t-s,n/\tau)x_{s}+\xi_{t}\
,$ (4)

where

$\displaystyle g_{j}(t,a)=\frac{a^{j}}{(j-1)!}t^{j-1}e^{-at}\ ,\ \ t\geq 0\ ,$
(5)

is the probability density function of the gamma distribution (more exactly,
the Erlang distribution I2013 ). Note that the lower limit of the integral in
Eq. (4) is sent to $-\infty$ since we will only be interested in the steady-
state regime. At the lowest order, the kernel $g_{1}(t,1/\tau)$ describes an
exponentially fading memory with a decay rate $\tau^{-1}$ (or a low pass
filter with bandwidth $\tau^{-1}$ in another language). For $n>1$, the kernel
has a maximum around $t=\tau$ and the peak becomes sharper as $n$ increases
(see e.g. Fig. 7.1 in S2011  or Fig. 2 in LHK2021 ). The discrete delay is
recovered in the limit $n\to\infty$. In the frequency domain (i.e., in Fourier
space222We here define the Fourier transform of a function $f(t)$ by
$f(\omega)=\int_{-\infty}^{+\infty}dt\>f(t)e^{i\omega t}$.) this simply
amounts to approximating the delay function $e^{i\omega\tau}$ as note1

$\displaystyle e^{i\omega\tau}\approx\frac{1}{(1-i\omega\tau/n)^{n}}\ .$ (6)

The Erlang density functions $g_{j}(t,a)$ satisfy the recursion relation
$dg_{j}(t,a)/dt=a[g_{j-1}(t,a)-g_{j}(t,a)]$ for $j\geq 1$ which is the basis
of the linear chain trick. Eq. (4) is then equivalent to the set of $n+1$
differential equations

$\displaystyle\dot{v}(t)$
$\displaystyle=-\frac{1}{Q_{0}}v(t)-x(t)+\frac{g}{Q_{0}}x_{n}(t)+\xi(t),$
$\displaystyle\dot{x}_{j}(t)$
$\displaystyle=\frac{n}{\tau}[x_{j-1}(t)-x_{j}(t)]\ ,\ \ j=1\cdots n\ ,$ (7)

where the auxiliary dynamical variables $x_{j}(t)$ are defined by

$\displaystyle x_{j}(t)$
$\displaystyle=\int_{-\infty}^{t}ds\>g_{j}(t-s,n/\tau)x(s)\ ,$ (8)

with $x_{0}(t)\equiv x(t)$.

Thanks to this alternative representation of the dynamics, we are now dealing
with a Markov process in the enlarged space
$\\{v(t),x(t),x_{1}(t),...x_{n}(t)\\}$, whereas the marginal dynamics of
$x_{t}$ of course remains non-Markovian. In general, the auxiliary variables
do not represent actual physical degrees of freedom but this may be the case
for $n$ small, in particular $n=1$ (see, e.g., Ref. CLSV2021 ). Note also that
$x_{j}(t)\approx x(t-j\tau/n)$ as $n\to\infty$ LK2019 , so that the linear
chain trick in the large-$n$ limit may be interpreted as a discretization of
the trajectory of the particle in the time interval $[t-\tau,t]$. This
illustrates the well-known fact that a delay-differential equation such as Eq.
(3) defines an infinite-dimensional dynamical system (since an infinite number
of initial conditions -actually, a function - is needed to uniquely specify
the time evolution) L2010 .

To simplify the forthcoming analysis, it is convenient to recast the set of
Eqs. (II.1) into a matrix form by introducing the $n+2$-dimensional vector
${\bf u}$ with components $u_{1}=v,u_{2}=x,u_{3}=x_{1},...,u_{n+2}=x_{n}$.
Defining the drift matrix

$\displaystyle A=\begin{bmatrix}-1/Q_{0}&-1&0&0&...&g/Q_{0}\\\
1&0&0&0&...&0\\\ 0&n/\tau&-n/\tau&0&...&0\\\ 0&0&n/\tau&-n/\tau&...&0&\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&...&0&n/\tau&-n/\tau&0&\\\
0&...&0&0&n/\tau&-n/\tau&\end{bmatrix}\ ,$ (9)

the equations in (II.1) become

$\displaystyle\dot{\bf u}_{t}=A{\bf u}_{t}+{\boldsymbol{\xi}}_{t}\ ,$ (10)

where ${\boldsymbol{\xi}}_{t}=(\xi_{t},0,0...0)^{T}$.

###  II.2 Dynamical observables

Assuming that the system has reached a NESS, we are interested in studying the
fluctuations of three time-integrated stochastic currents. These are (in
reduced units):

a) the work done by the feedback force during the time window $[0,t]$,

$\displaystyle\beta{\cal
W}_{t}=\frac{2g}{Q_{0}^{2}}\int_{0}^{t}x_{n}(t^{\prime})\circ dx(t^{\prime})\
,$ (11)

where $\beta=(k_{B}T)^{-1}$ and $\circ$ denotes the Stratonovich product,

b) the corresponding heat dissipated into the environment S2010

$\displaystyle\beta{\cal Q}_{t}$
$\displaystyle=\frac{2}{Q_{0}}\int_{0}^{t}[\frac{1}{Q_{0}}v(t^{\prime})-\xi(t^{\prime})]\circ
dx(t^{\prime})$ $\displaystyle=\beta{\cal
W}_{t}-\frac{2}{Q_{0}}\int_{0}^{t}\left[x(t^{\prime})\circ
dx(t^{\prime})+v(t^{\prime})\circ dv(t^{\prime})\right]$
$\displaystyle=\beta{\cal W}_{t}-\Delta E\ ,$ (12)

where $\Delta E=(1/Q_{0})[x^{2}_{t}-x^{2}_{0}+v^{2}_{t}-v^{2}_{0}]$ is the
change in the internal energy of the system,

c) the entropy production (EP)

$\displaystyle\Sigma_{t}$ $\displaystyle=\beta{\cal
Q}_{t}+\ln\frac{p(x_{0},v_{0})}{p(x_{t},v_{t})}\ ,$ (13)

where $\beta{\cal Q}_{t}$ is the entropy change in the medium and $p(x,v)$ is
the stationary PDF (see Eq. (26) below). Specifically,

$\displaystyle\Sigma_{t}$ $\displaystyle=\beta{\cal
W}_{t}+\frac{1}{Q_{0}}\big{[}(\frac{T}{T_{x}}-1)(x^{2}_{t}-x^{2}_{0})+(\frac{T}{T_{v}}-1)(v^{2}_{t}-v^{2}_{0})\big{]}\
,$ (14)

where $T_{x}=(2T/Q_{0})\langle x_{t}^{2}\rangle$ and $T_{v}=(2T/Q_{0})\langle
v_{t}^{2}\rangle$ are the configurational and kinetic temperatures of the
system, respectively RMT2015  (angle brackets indicate a steady-state
average).

If the dynamics were Markovian, the second term in Eq. (13) would correspond
to the change in the Shannon entropy of the system C1999 ; S2005  and
$\Sigma_{t}$ would be the total stochastic EP in the time interval $[0,t]$.
Here, $\Sigma_{t}$ defines the “apparent” EP that an observer unaware of the
existence of the non-Markovian feedback would regard as the total EP. This
quantity can be extracted from the stochastic trajectories collected in
experiments, as done in Ref. DRLK2022 . Note in passing that the linearity of
the Langevin equation (3) implies that the probabilities of a trajectory and
its time reversal are Gaussian and equal in the NESS. In consequence, the log
ratio of these two probabilities, which is usually taken as the definition of
EP in stochastic thermodynamics (see Ref. PP2021  and references therein), is
zero and does not properly accounts for the irreversible character of the
feedback process RMT2015 .

Since the three observables only differ by temporal boundary terms, they share
the same average rate in the NESS,

$\displaystyle\langle\beta\dot{\cal W}_{t}\rangle=\langle\beta\dot{\cal
Q}_{t}\rangle=\langle\dot{\Sigma}_{t}\rangle=\frac{2}{Q_{0}}[\frac{1}{Q_{0}}\langle
v_{t}^{2}\rangle-\langle\xi_{t}v_{t}\rangle]=\frac{1}{Q_{0}}(\frac{T_{v}}{T}-1)\
.$ (15)

This expression illustrates the fact that the exchange of heat with the
thermal environment in an underdamped system proceeds through the kinetic
energy of the system independently of the form of the potential function S2010
. On the other hand, the fluctuations of the observables may differ, as shown
experimentally in Ref. DRLK2022  in the case of the discrete delay. In
particular, one has

$\displaystyle\langle e^{-\beta{\cal Q}_{t}}\rangle=e^{t/Q_{0}}$ (16)

at all times RTM2016 , whereas it is conjectured that

$\displaystyle\langle e^{-\Sigma_{t}}\rangle\sim e^{\dot{S}_{\cal J}t}$ (17)

as $t\to\infty$, where $\dot{S}_{\cal J}$ is a nontrivial quantity related to
the Jacobian originating from time reversal MR2014 ; RTM2017 .

In the following it will be convenient to collectively define the three
stochastic currents by

$\displaystyle{\cal A}_{t}=\int_{0}^{t}{\bf g}_{\rm o}({\bf
u}_{t^{\prime}})\circ d{\bf u}_{t^{\prime}}\ ,$ (18)

where333With the present definition (see also Ref. RTM2017 ), ${\cal A}_{t}/t$
is intensive in time and converges in probability to a constant as
$t\to\infty$. This differs from the definition adopted in Refs. CT2015 ; T2018
where ${\cal A}_{t}$ itself is the intensive quantity.

$\displaystyle{\bf g}_{\rm o}({\bf u})=B_{\rm o}{\bf u}\ ,$ (19)

and $B_{\rm o}$ is a matrix of dimension $(n+2)\times(n+2)$. (Henceforth, the
subscript $\rm o$ refers to the observable, with ${\rm o}=w,q$ and $\sigma$
for the work, the heat, and the entropy production, respectively.) From Eqs.
(11)-(14), we have

$\displaystyle B_{\rm o}$ $\displaystyle=B_{w}+S_{\rm o}\ ,$ (20)

where

$\displaystyle B_{w}=\frac{2g}{Q_{0}^{2}}\begin{bmatrix}0&0&0&...&0\\\
0&0&0&...&1\\\ 0&0&0&...&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&...&0&\end{bmatrix}\ ,$ (21)

and

$\displaystyle S_{q}=-\frac{2}{Q_{0}}\begin{bmatrix}1&0&0&...&0\\\
0&1&0&...&0\\\ 0&0&0&...&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&...&0&\end{bmatrix}\ ,$ (22)

$\displaystyle
S_{\sigma}=\frac{2}{Q_{0}}\begin{bmatrix}\frac{T}{T_{v}}-1&0&0&...&0\\\
0&\frac{T}{T_{x}}-1&0&...&0\\\ 0&0&0&...&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&0&...&0&\end{bmatrix}\ .$ (23)

For convenience, we also define the matrix $S_{w}$ as the null matrix. Note
that the matrices $B_{\rm o}$ characterizing the observables have the same
anti-symmetric part. This will have an important consequence in the following.

###  II.3 Stability of the non-equilibrium steady state

In the present study, we assume that the system has reached a stable NESS. For
given values of the parameters $(Q_{0},g,n)$, this requires to choose the
delay $\tau$ appropriately. Indeed, as is well known, a time delay induces a
complex dynamical behavior N2001 . In particular, there may be a series of
stability switches as $\tau$ increases, corresponding to
destabilizing/stabilizing Hopf bifurcations DVGW1995 . As usual with linear
stochastic systems, the boundaries of the domain of stability in the parameter
space can be determined by computing the roots of the characteristic
polynomial $p_{A}(s)$ of the drift matrix $A$. The loss of stability is then
associated with the occurrence of a root with a positive real part. From Eq.
(9), a straightforward calculation yields

$\displaystyle
p_{A}(s)\equiv\det(sI_{n+2}-A)=(\frac{n}{\tau})^{n}[(s^{2}+\frac{s}{Q_{0}}+1)(1+\frac{s\tau}{n})^{n}-\frac{g}{Q_{0}}]\
.$ (24)

We thus have to deal with a polynomial of degree $n+2$, for which there exist
powerful root-finding algorithms. As
$\lim_{n\to\infty}(1+\frac{s\tau}{n})^{n}=e^{s\tau}$ and

$\displaystyle\lim_{n\to\infty}(\frac{n}{\tau})^{-n}p_{A}(s)=(s^{2}+\frac{s}{Q_{0}}+1)e^{s\tau}-\frac{g}{Q_{0}}\
,$ (25)

this may be compared with the task of solving a transcendental equation in the
case of the discrete delay (see Appendix C of Ref. RMT2015 ).

However, the calculation of the stability diagram in the whole parameter space
$(Q_{0},g,\tau,n)$ is a formidable task which is beyond the scope of the
present work. We will content ourselves with the example shown in Fig. 1 that
illustrates the influence of $n$. This corresponds to a case for which only
two stability domains exist for $n$ finite. In this figure, the kinetic
temperature $T_{v}=(2T/Q_{0})\langle v^{2}\rangle$ (with $\langle
v^{2}\rangle$ given by the element $(1,1)$ of the covariance matrix $\Sigma$,
see Eq. (27) below), is plotted as a function of $\tau$, and the loss of
stability of the stationary state is signaled by the divergence of
$T_{v}~{}$RMT2015 . It can be seen that the width of the unstable region
increases with $n$ and that only the first stability domain survives for
$n=\infty$444Moreover, the critical delay $\tau_{c,n}^{(1)}$ corresponding to
the first destabilizing Hopf bifurcation is numerically found to decreases to
its limit $\tau^{(1)}_{c,\infty}$ as $1/n$.. This is in line with the general
lore that a discrete delay is more destabilizing than a distributed delay
S2011 ; B1989 ; CJ2009 ; C2010 ; ABG2020  (see also Ref. LHK2021 ).

Figure 1:  (Color on line) Kinetic temperature $T_{v}$ in the NESS as a
function of $\tau$ for $Q_{0}=4$, $g=3.9$, and different values of $n$: $n=25$
(red solid line), $n=50$ (blue), $n=75$ (green), $n=100$ (black), and
$n=\infty$ (black dashed line). For $n=25$, the stationary state is stable for
$0\leq\tau<2.79$ and $\tau>8.07$. The second stability region is displaced to
larger $\tau$ as $n$ increases and no longer exists for $n=\infty$ (discrete
delay). In the latter case, the system is stable for
$0\leq\tau<\tau_{c,\infty}^{(1)}\simeq 2.53$ only.

Due to the linearity of Eq. (10), the probability distribution function (pdf)
in the NESS is a multivariate Gaussian distribution,

$\displaystyle p({\bf
u})=\frac{1}{\sqrt{(2\pi)^{n+2}\det\Sigma}}e^{-\frac{1}{2}{\bf
u}^{T}.\Sigma^{-1}.{\bf u}}\ ,$ (26)

where $\Sigma$, the $(n+2)\times(n+2)$ covariance matrix (not to be confused
with the entropy production defined above), is solution of the Lyapunov
equation G2004

$\displaystyle A\Sigma+\Sigma A^{T}=-D\ ,$ (27)

with $D={\bf d}.{\bf d}^{T}$ has only one non-zero element with ${\bf
d}^{T}=(1,0,0...0)$. To solve this equation for large $n$, one can use the
property that the eigenvalues of $\Sigma$, ordered as
$\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{n+2}$, decay very fast as their
order increases. Therefore, although $\Sigma$ is full-rank in the stability
region, the numerical rank is very low. This feature is commonly encountered
in large-scale Lyapunov equations when the matrix in the right-hand side of
the equation has a low rank P2000 ; ASZ2002 . Such equations typically arise
in the study of the controllability and observability of linear time-invariant
(LTI) dynamical systems, with the matrix in the right-hand side having a rank
equal to the number of inputs or outputs in the system B2021 . One can then
use accurate low-rank approximations to the solution and consider large values
of $n$ BS2013 ; S2016 . When extremely high precision is required, one can
also use the closed-form expression of $\Sigma$ in terms of the Cauchy matrix
built from the eigenvalues of $A$ P2000 ; ASZ2002 ; S2016 .

##  III Matrix Riccati equations and their properties

###  III.1 Basic equations

We now come to the heart of the matter and derive the Riccati differential
equations (RDEs) that will allow us to investigate the time evolution of the
moment generating functions $G_{\rm o,\lambda}(t)\equiv\langle
e^{-\lambda{\cal A}_{t}}\rangle$ in the steady state and to extract the
complete asymptotic form at large time:

$\displaystyle G_{\rm o,\lambda}(t)\sim g_{\rm o}(\lambda)e^{\mu_{\rm
o}(\lambda)t}\ ,$ (28)

where

$\displaystyle\mu_{\rm o}(\lambda)\equiv\lim_{t\to\infty}\frac{1}{t}\ln G_{\rm
o,\lambda}(t)$ (29)

is the SCGF and $g_{\rm o}(\lambda)$ is the pre-exponential factor.

As will be discussed in detail in the following (see Sec. IV.2.1), the domain
of definition of the SCGF (i.e., the values of $\lambda$ for which $\mu_{\rm
o}(\lambda)<\infty$) may depend on the observable. This is due to rare but
large fluctuations of the initial or final points of the stochastic
trajectories that induce singularities in the pre-exponential factors.
Therefore, the knowledge of both $\mu_{\rm o}(\lambda)$ and $g_{\rm
o}(\lambda)$ is required to obtain the large-time behavior of the pdf $P_{\rm
o}(a,t)\equiv\langle\delta({\cal A}_{t}-at)\rangle$ characterized by the rate
function

$\displaystyle I_{\rm o}(a)\equiv-\lim_{t\to\infty}\frac{1}{t}\ln P_{\rm
o}(a,t)\ .$ (30)

Large deviation theory tells us that if $\mu_{\rm o}(\lambda)$ exists and is
differentiable, then $I_{\rm o}(a)$ is given by the Legendre-Fenchel transform
T2009 ; J2020

$\displaystyle I_{\rm o}(a)=\max_{\lambda}[-\lambda a-\mu_{\rm o}(\lambda]\ .$
(31)

In order to compute $G_{\rm o,\lambda}(t)$, it is convenient to start from the
restricted generating function $G_{\rm o,\lambda}({\bf u},t|{\bf
u}_{0})=\langle e^{-\lambda{\cal A}_{t}}\rangle_{{\bf u}_{0},{\bf u}}$, where
the initial and final configurations of the trajectories of duration $t$ are
fixed at ${\bf u}_{0}$ and ${\bf u}$, respectively. By definition, $G_{\rm
o,\lambda}({\bf u},0|{\bf u}_{0})=\delta({\bf u}-{\bf u}_{0})$ and

$\displaystyle G_{\rm o,\lambda}(t)$ $\displaystyle=\int d{\bf u}_{0}\>p({\bf
u}_{0})\int d{\bf u}\>G_{\rm o,\lambda}({\bf u},t|{\bf u}_{0})$ (32)

provided the integrals over ${\bf u}_{0}$ and ${\bf u}$ converge. This
suggests considering the time evolution of

$\displaystyle G^{r}_{\rm o,\lambda}({\bf u}_{0},t)$ $\displaystyle=\int d{\bf
u}\>G_{\rm o,\lambda}({\bf u},t|{\bf u}_{0})$ (33)

and

$\displaystyle G^{l}_{\rm o,\lambda}({\bf u},t)$ $\displaystyle=\int d{\bf
u}_{0}\>p({\bf u}_{0})G_{\rm o,\lambda}({\bf u},t|{\bf u}_{0})$ (34)

separately. They satisfy the initial conditions

$\displaystyle G^{r}_{\rm o,\lambda}({\bf u}_{0},0)=1$ (35)

and

$\displaystyle G^{l}_{\rm o,\lambda}({\bf u},0)=p({\bf u})\ .$ (36)

(The superscript r and l denote “right” and “left”, respectively. This
notation will be justified later when considering the long-time limit.)

Standard application of the Feynman-Kac formula (see, e.g., Ref. M2005  for a
pedagogical review) shows that $G^{r}_{\rm o,\lambda}({\bf u}_{0},t)$ evolves
in time according to the backward Fokker-Planck equation

$\displaystyle\partial_{t}G^{r}_{\rm o,\lambda}({\bf u}_{0},t)={\cal L}_{\rm
o,\lambda}G^{r}_{\rm o,\lambda}({\bf u}_{0},t)\ ,$ (37)

where ${\cal L}_{\rm o,\lambda}$ is the so-called tilted (or biased) generator
given by CT2015 ; T2018

$\displaystyle{\cal L}_{\rm o,\lambda}={\bf F}\cdot({\bf\nabla}-\lambda{\bf
g}_{\rm o})+\frac{1}{2}(\nabla-\lambda{\bf g}_{\rm o})\cdot
D(\nabla-\lambda{\bf g}_{\rm o})\ ,$ (38)

with ${\bf F}=A{\bf u}$ in the case at hand. Explicit expressions of ${\cal
L}_{\rm o,\lambda}$ for ${\rm o}=w,q,\sigma$ are given in Appendix A.
Likewise, $G^{l}_{\rm o,\lambda}({\bf u},t)$ satisfies the forward Fokker-
Planck equation

$\displaystyle\partial_{t}G^{l}_{\rm o,\lambda}({\bf u},t)={\cal L}_{\rm
o,\lambda}^{{\dagger}}G^{l}_{\rm o,\lambda}({\bf u},t)\ ,$ (39)

where ${\cal L}_{\rm o,\lambda}^{\dagger}$ is the adjoint of ${\cal L}_{\rm
o,\lambda}$.

Solving such linear partial differential equations beyond the long-time limit
requires one to determine the whole spectrum of the operators ${\cal L}_{\rm
o,\lambda}$ and ${\cal L}_{\rm o,\lambda}^{\dagger}$ and the associated
eigenfunctions, which is a daunting or even impossible task. In the present
case, however, both the drift ${\bf F}$ and the vector function ${\bf g}_{\rm
o}$ depend linearly on ${\bf u}$ so that we can anticipate that the solutions
of Eqs. (37) and (39) with initial conditions (35) and (36), respectively, are
just multivariate Gaussians. Specifically, we show in Appendix B.0.1 that

$\displaystyle G^{r}_{\rm o,\lambda}({\bf u}_{0},t)$
$\displaystyle=\exp\left(-\frac{1}{2}{\bf u}_{0}^{T}C^{r}_{\rm
o}(\lambda,t){\bf u}_{0}+\int_{0}^{t}f^{r}_{\rm
o}(\lambda,t^{\prime})dt^{\prime}\right)$ (40a) $\displaystyle G^{l}_{\rm
o,\lambda}({\bf u},t)$
$\displaystyle=\left[(2\pi)^{n+2}\det\Sigma\right]^{-1/2}\exp\left(-\frac{1}{2}{\bf
u}^{T}C^{l}_{\rm o}(\lambda,t){\bf u}+\int_{0}^{t}f^{l}_{\rm
o}(\lambda,t^{\prime})dt^{\prime}\right)\ ,$ (40b)

with

$\displaystyle f^{r}_{\rm o}(\lambda,t)$
$\displaystyle=-\frac{1}{2}\mbox{Tr}[D(C_{\rm o}^{r}(\lambda,t)+\lambda B_{\rm
o})]$ (41a) $\displaystyle f^{l}_{\rm o}(\lambda,t)$
$\displaystyle=-\frac{1}{2}\mbox{Tr}[D(C_{\rm o}^{l}(\lambda,t)-\lambda B_{\rm
o})]-\mbox{Tr}(A)\ .$ (41b)

$C^{r}_{\rm o}(\lambda,t)$ and $C^{l}_{\rm o}(\lambda,t)$ are symmetric
matrices of dimension $(n+2)\times(n+2)$ which are solutions of the matrix
differential equations

$\displaystyle\dot{C}_{\rm o}^{r}(\lambda,t)$ $\displaystyle={\cal R}_{\rm
o,\lambda}\big{[}C_{\rm o}^{r}(\lambda,t)\big{]}$ (42a)
$\displaystyle\dot{C}_{\rm o}^{l}(\lambda,t)$ $\displaystyle={\cal R}_{\rm
o,\lambda}\big{[}-C_{\rm o}^{l}(\lambda,t)\big{]}$ (42b)

with initial conditions

$\displaystyle C_{\rm o}^{r}(\lambda,0)$ $\displaystyle=0$ (43a)
$\displaystyle C_{\rm o}^{l}(\lambda,0)$ $\displaystyle=C\ ,$ (43b)

respectively, where $C$ is the inverse of the covariance matrix $\Sigma$. In
these equations, ${\cal R}_{\rm o,\lambda}$ is a quadratic Riccati operator
acting on a general $(n+2)\times(n+2)$ matrix $X$ as

$\displaystyle{\cal R}_{\rm o,\lambda}[X]=A_{\rm o}(\lambda)^{T}X+XA_{\rm
o}(\lambda)-XDX+K_{\rm o}(\lambda)\ ,$ (44)

with

$\displaystyle A_{\rm o}(\lambda)=A-\lambda DB_{\rm o}\ ,$ (45)

and

$\displaystyle K_{\rm o}(\lambda)=\lambda(A^{T}B_{\rm o}+B_{\rm
o}^{T}A)-\lambda^{2}B^{T}_{\rm o}DB_{\rm o}\ .$ (46)

The explicit expressions of the symmetric matrices $K_{\rm o}(\lambda)$ are
given in Appendix B.0.2. Furthermore, since $D_{ij}=\delta_{i1}\delta_{j1}$
and $\mbox{Tr}(A)=-(1/Q_{0}+n^{2}/\tau)$ in the present model, Eqs. (41)
become

$\displaystyle f^{r}_{\rm o}(\lambda,t)$ $\displaystyle=-\frac{1}{2}(C_{\rm
o,11}^{r}(\lambda,t)+\lambda B_{\rm o,11})$ (47a) $\displaystyle f^{l}_{\rm
o}(\lambda,t)$ $\displaystyle=-\frac{1}{2}(C_{\rm o,11}^{l}(\lambda,t)-\lambda
B_{\rm o,11})+\frac{1}{Q_{0}}+\frac{n^{2}}{\tau}\ .$ (47b)

By integrating $G^{r}_{\rm o,\lambda}({\bf u}_{0},t)p({\bf u}_{0})$ over ${\bf
u}_{0}$ and $G^{l}_{\rm o,\lambda}({\bf u},t)$ over ${\bf u}$, we then obtain
two expressions of the generating function $G_{\rm o,\lambda}(t)$,

$\displaystyle G_{\rm o,\lambda}(t)$
$\displaystyle=\Big{[}\frac{\det\big{(}C^{r}_{\rm o}(\lambda,t)+C\big{)}}{\det
C}\Big{]}^{-1/2}\exp\big{(}\int_{0}^{t}f^{r}_{\rm
o}(\lambda,t^{\prime})dt^{\prime}\big{)}$ (48a)
$\displaystyle=\Big{[}\frac{\det C^{l}_{\rm o}(\lambda,t)}{\det
C}\Big{]}^{-1/2}\>\exp\big{(}\int_{0}^{t}f^{l}_{\rm
o}(\lambda,t^{\prime})dt^{\prime}\big{)}\ .$ (48b)

By construction, these two expressions555Eq. (48b) is similar to the
expression of the generating function derived by C. Kwon et al. KNP2011  via a
path-integral method for the non-equilibrium work in linear diffusion systems.
The role of the matrix $C^{l}(\lambda,t)$ is played by a matrix
$\tilde{A}(\tau,\lambda)$ that obeys a non-linear matrix equation similar to
our Eq. (42b) (with the matrix $\Lambda$ playing the role of the matrix
$K_{w}(\lambda)$). Likewise, Eq. (48a) corresponds to Eq. (43) in the recent
paper of J. du Buisson and H. Touchette DBT2023 , with $-(1/2)C_{\rm
o}^{r}(\lambda,t)$ replaced by the matrix $B_{k}(t)$. From the dictionary
$\lambda\to-k$, $A\to-M$, $B_{\rm o}\to\Gamma$, one can easily check that the
RDE (42a) corresponds to Eq. (60) in Ref. DBT2023  (with $\Gamma^{T}=-\Gamma$
in Eq. (60) because $\Gamma$ is assumed to be antisymmetric). give the same
result as long as the solutions of Eqs. (42) exist (see the discussion in the
next section). But for the generating function $G_{\rm o,\lambda}(t)$ to be
finite it is mandatory that the matrices $C+C_{\rm o}^{r}(\lambda,t)$ and
$C_{\rm o}^{l}(\lambda,t)$ are positive definite. As we shall see later in
Sec. IV.1, this crucial condition is not always satisfied. Moreover, the
equivalence between Eqs. (48a) and (48b) may result from a nontrivial
mathematical mechanism.

Finally, let us note that an alternative expression of $G_{\rm o,\lambda}(t)$
is available when the matrix $C_{\rm o}^{l}(\lambda,t^{\prime})$ is invertible
for all $t^{\prime}\in[0,t]$. If so, a few manipulations detailed in Appendix
B.0.3 lead to666A similar, but slightly more complicated expression of the
generating function can be obtained in terms of the inverse of the matrix
$C+C^{r}_{\rm o}(\lambda,t)$, when the latter exists.

$\displaystyle G_{\rm o,\lambda}(t)$
$\displaystyle=\exp\left(-\frac{1}{2}\int_{0}^{t}\mbox{Tr}\big{(}K_{\rm
o}(\lambda)\Sigma^{l}_{\rm o}(\lambda,t^{\prime})+\lambda DB_{\rm
o}\big{)}dt^{\prime}\right)\ ,$ (49)

where $\Sigma^{l}_{\rm o}(\lambda,t)$, the inverse matrix of $C_{\rm
o}^{l}(\lambda,t)$, is solution of the complementary RDE

$\displaystyle\frac{\partial}{\partial t}\Sigma^{l}_{\rm o}(\lambda,t)=A_{\rm
o}(\lambda)\Sigma^{l}_{\rm o}(\lambda,t)+\Sigma^{l}_{\rm
o}(\lambda,t)A^{T}_{\rm o}(\lambda)-\Sigma^{l}_{\rm o}(\lambda,t)K_{\rm
o}(\lambda)\Sigma^{l}_{\rm o}(\lambda,t)+D\ ,$ (50)

with initial condition

$\displaystyle\Sigma^{l}_{\rm o}(\lambda,0)=\Sigma\ .$ (51)

In particular, using the fact that $B_{q,11}=-2/Q_{0}$ and $K_{q}(1)=0$ [Eq.
(B14)], Eq. (49) readily yields

$\displaystyle G_{q,\lambda=1}(t)$ $\displaystyle=e^{t/Q_{0}}\ ,$ (52)

which is the universal IFT for the fluctuating heat in underdamped Langevin
processes RTM2016 . This result is also directly obtained from Eq. (42a) since
the unique solution of the initial value problem is
$C^{r}_{q}(\lambda=1,t)=0$, which implies that
$f^{r}_{q}(\lambda=1,t)=1/Q_{0}$ from Eq. (47a).

###  III.2 Properties of the Riccati differential equations (RDEs)

RDEs similar to Eqs. (42), as well as the corresponding continuous algebraic
Riccati equations (CAREs) whose solutions are stationary solutions of the RDEs
(see Sec. III.3 below), appear in many branches of applied mathematics, the
most prominent application being linear optimal control and filtering problems
B2021 . Within this framework, several important issues have been extensively
discussed in the literature such as the global existence of the solutions as
one varies the coefficients of the differential equation or the initial data,
the convergence toward a particular solution of the corresponding CARE as
$t\to\infty$, and the mechanism of attraction (see Ref. AFIJ2003  and
references therein). In particular, it is a standard result that the solution
of the RDE (which is unique for a given initial condition) exists for
$t\in[0,\infty)$ and is symmetric, positive semidefinite if the source term
(in the present case, the matrix $K_{\rm o}(\lambda)$ in Eq. (44)), the
quadratic term (i.e., the term $XDX$), and the initial condition are positive
semidefinite. With additional conditions on the coefficients, it is also
proven that the solution converges monotonically to the maximal solution
$X^{+}$ of the CARE (to be defined later) which turns out to be the unique
symmetric positive semidefinite solution777From the viewpoint of control
theory, the important feature is that the maximal solution $X^{+}$ is
“stabilizing”, i.e., all the eigenvalues of the “closed-loop” matrix $A_{\rm
o}(\lambda)-DX^{+}$ have negative real parts (see e.g. Ref. K2010 )..

Unfortunately, it is readily seen from Eqs. (B13)-(B15) that neither
$K_{w}(\lambda)$ nor $K_{\sigma}(\lambda)$ are positive semidefinite ($n$
eigenvalues are equal to $0$, one is positive, and one is negative). Only
$K_{q}(\lambda)\geq 0$ for $0\leq\lambda\leq 1$. This is an essential
difference with the standard situation treated in linear quadratic (LQ)
optimal control, and it significantly complicates the present study note3 .
The consequences that will be explored in the rest of this paper and
illustrated numerically are the following:

1) the solutions of Eqs. (42) may exhibit a finite-time escape phenomenon,
which means that they may blow up in a finite time,

2) they may fail to converge (i.e., the solution may oscillate),

3) they may converge to a solution of the CARE which is not the maximal
solution,

4) the maximal solution is not automatically positive semidefinite.

On the positive side, the Riccati operator defined by Eq. (44) has a
remarkable property that holds for arbitrary matrices $A$, $B_{\rm o}$, and
$D$ symmetric. Indeed, if $S$ is a symmetric matrix, then

$\displaystyle{\cal R}_{\rm o,\lambda}[X]={\cal R}^{\prime}_{\rm
o,\lambda}[X-\lambda S]\ ,$ (53)

where ${\cal R}^{\prime}_{\rm o,\lambda}$ is the modified operator obtained by
changing $B_{\rm o}$ into $B^{\prime}_{\rm o}=B_{\rm o}+S$. Therefore, when
the matrices $B_{\rm o}$ characterizing the various observables only differ by
their symmetric part $S_{\rm o}$, which is the case for the three observables
${\cal W}_{t},{\cal Q}_{t},\Sigma_{t}$ under consideration, one has

$\displaystyle{\cal R}_{\rm o,\lambda}[X]={\cal R}_{\rm
o^{\prime},\lambda}[X+\lambda(S_{\rm o}-S_{\rm o^{\prime}})]\ .$ (54)

This “invariance” property allows one to compute all matrices $C^{r}_{\rm
o^{\prime}}(\lambda,t)$ and $C^{l}_{\rm o^{\prime}}(\lambda,t)$ by solving the
RDEs corresponding to a single operator ${\cal R}_{\rm o,\lambda}$ but with
different initial conditions. Specifically, by letting $X_{\rm
o^{\prime}}(\lambda,t)$ be the solution of the RDE $\dot{X}_{\rm
o^{\prime}}={\cal R}_{\rm o,\lambda}[X_{\rm o^{\prime}}]$ with initial
condition

$X_{\rm o^{\prime}}(\lambda,0)=\lambda(S_{\rm o^{\prime}}-S_{\rm o})$ (55)

we obtain that

$C^{r}_{\rm o^{\prime}}(\lambda,t)=X_{\rm
o^{\prime}}(\lambda,t)-\lambda(S_{\rm o^{\prime}}-S_{\rm o})$ (56)

and, similarly, if $X_{\rm o^{\prime}}(\lambda,t)$ is the solution of the RDE
$\dot{X}_{\rm o^{\prime}}={\cal R}_{\rm o,\lambda}[-X_{\rm o^{\prime}}]$ with
initial condition

$X_{\rm o^{\prime}}(\lambda,0)=C-\lambda(S_{\rm o^{\prime}}-S_{\rm o})$ (57)

we have that

$C^{l}_{\rm o^{\prime}}(\lambda,t)=X_{\rm
o^{\prime}}(\lambda,t)+\lambda(S_{\rm o^{\prime}}-S_{\rm o}).$ (58)

####  III.2.1 Global existence of the solutions

We begin our study of the solutions of the RDEs (42) by briefly discussing the
issue of their global existence. First of all, we note that the non-negativity
of the diffusion matrix $D$ implies that the solutions are bounded from above
by the solutions of the corresponding Lyapunov differential equations (i.e.,
Eqs. (42) with $D=0$) note4 . Since these upper bounds do not blow up in
finite time, the finite-time escape phenomenon, when it occurs, is due to the
absence of a lower bound and manifests itself by the divergence of the
smallest eigenvalue of $C^{r}_{\rm o}(\lambda,t)$ or $C^{l}_{\rm
o}(\lambda,t)$ toward $-\infty$. This is a crucial observation because it
means that this eigenvalue, which is initially positive (as $C+C^{r}_{\rm
o}(\lambda,0)=C^{l}_{\rm o}(\lambda,0)=C$), vanishes before diverging to
$-\infty$. Therefore, this finite-time escape phenomenon is always preceded by
the divergence of the generating function $G_{\rm o,\lambda}(t)$888Note that
the consistency between the two expressions of the generating function [Eqs.
(48a) and (48b)] imposes that the determinants of $C^{r}_{\rm o}(\lambda,t)+C$
and $C^{l}_{\rm o}(\lambda,t)$ vanish simultaneously, as will be observed
later in the numerical examples (see, e.g., Fig. 10 in Sec. IV.2.1)..

Let us focus on the case $0\leq\lambda\leq 1$ for which we can take advantage
of the semi-positiveness of the source term $K_{q}(\lambda)$ to draw definite
conclusions about the global existence of the solutions of the RDEs and of the
generating function. It suffices to consider the “right” matrices $C^{r}_{\rm
o}(\lambda,t)$ since the condition $C+C^{r}_{\rm o}(\lambda,t)>0$ for all
$t\in[0,\infty)$ implies that $C^{l}_{\rm o}(\lambda,t)>0$ (otherwise, the two
expressions (48a) and (48b) of $G_{\rm o,\lambda}(t)$ would not be
consistent).

First of all, since $C^{r}_{q}(\lambda,0)=0$ and the RDE $\dot{X}={\cal
R}_{q,\lambda}[X]$ has all the properties of the RDEs encountered in LQ
optimal control AFIJ2003 , we can immediately assert that the matrix
$C^{r}_{q}(\lambda,t)$ exists and is positive semidefinite on $[0,\infty)$.
Moreover, $\dot{C}^{r}_{q}(\lambda,0)=K_{q}(\lambda)\geq 0$ and thus
$C^{r}_{q}(\lambda,t)$ is monotonically non-decreasing note80 . Hence
$C^{r}_{q}(\lambda,t)+C>0$ and we conclude that the generating function
$G_{q,\lambda}(t)$ is always finite.

Reaching a conclusion about the existence and positiveness of the matrices
$C+C^{r}_{w}(\lambda,t)$ and $C+C^{r}_{\sigma}(\lambda,t)$ is less
straightforward. We first note from Eqs. (22) and (23) that $S_{q}\leq 0$ and
$S_{\sigma}-S_{q}\geq 0$. We then exploit the order-preserving property of
Riccati differential equations which states that the solutions depend
monotonically on the initial values note8 . Accordingly, the solutions
$X_{w}(\lambda,t)$ and $X_{\sigma}(\lambda,t)$ of the RDE $\dot{X}={\cal
R}_{q,\lambda}[X]$, which correspond respectively to the initial conditions
$X_{w}(\lambda,0)=-\lambda S_{q}\geq 0$ and
$X_{\sigma}(\lambda,0)=\lambda(S_{\sigma}-S_{q})\geq 0$, satisfy
$X_{w}(\lambda,t)\geq C^{r}_{q}(\lambda,t)$ and $X_{\sigma}(\lambda,t)\geq
C^{r}_{q}(\lambda,t)$. From the “invariance” relations (56) and (58) (with
$S_{w}=0$) we then obtain

$\displaystyle C_{w}^{r}(\lambda,t)=X_{w}(\lambda,t)+\lambda S_{q}\geq
C^{r}_{q}(\lambda,t)+\lambda S_{q}$ (59)

and

$\displaystyle
C_{\sigma}^{r}(\lambda,t)=X_{\sigma}(\lambda,t)-\lambda(S_{\sigma}-S_{q})\geq
C^{r}_{q}(\lambda,t)-\lambda(S_{\sigma}-S_{q})\ .$ (60)

Furthermore, it can be shown that

$\displaystyle\det[C+\lambda S_{q}]$ $\displaystyle=\det[I_{n+2}+\lambda
S_{q}\Sigma]\>\det
C=(1-\lambda\frac{T_{v}}{T})(1-\lambda\frac{T_{x}}{T})\>\det C$ (61)

and

$\displaystyle\det[C+\lambda(S_{q}-S_{\sigma})]=(1-\lambda)^{2}\>\det C\ .$
(62)

As a result, $C+\lambda S_{q}>0$ for $\lambda<\min(T/T_{x},T/T_{v})$ and
$C+\lambda(S_{q}-S_{\sigma})>0$ for $\lambda<1$ (since all eigenvalues of
these matrices are positive for $\lambda=0$ and remain positive as long as the
respective determinants do not vanish). From this set of inequalities we
conclude that

$\displaystyle C^{r}_{w}(\lambda,t)+C\geq C^{r}_{q}(\lambda,t)+C+\lambda
S_{q}>0\ ,\>\>0\leq\lambda<\min(1,T/T_{x},T/T_{v})$ (63)

and

$\displaystyle C^{r}_{\sigma}(\lambda,t)+C\geq
C^{r}_{q}(\lambda,t)+C-\lambda(S_{\sigma}-S_{q})>0\ ,\>\>0\leq\lambda<1\ .$
(64)

In consequence, the corresponding generating functions $G_{w,\lambda}(t)$ and
$G_{\sigma,\lambda}(t)$ are always finite for values of $\lambda$ within the
above ranges. We stress that these conditions are sufficient but not
necessary. Note also that the inequality (63) still holds for $\lambda=1$ if
both $T_{x}$ and $T_{v}$ are smaller than $T$. On the other hand, the strict
inequality (64) is replaced by $C^{r}_{\sigma}(1,t)+C\geq 0$. This special but
important case $\lambda=1$ will be treated in detail in Sec. V.

####  III.2.2 Solutions of the RDEs and asymptotic behavior

Various methods for solving RDEs are available in the literature, including
for large-scale problems BM2004 . Here, we will use the classical approach
that consists in transforming each quadratic differential equation into a
linear system of first-order Hamiltonian differential equations of double size
AFIJ2003 . We can then obtain a closed-form representation of the solution
which is suitable for analyzing the asymptotic behavior and the dependence on
the initial condition. Although this procedure is standard, it is worthwhile
to replicate the derivation in the case at hand.

Consider first Eq. (42a) with initial condition (43a). It can be checked by
direct substitution that the solution can be expressed as

$\displaystyle C^{r}_{{\rm o}}(\lambda,t)$ $\displaystyle=V^{r}_{\rm
o}(\lambda,t)U^{r}_{\rm o}(\lambda,t)^{-1}\ ,$ (65)

where the matrices $U^{r}_{\rm o}(\lambda,t)$ and $V^{r}_{\rm o}(\lambda,t)$
are solutions of the linear system

$\displaystyle\begin{bmatrix}\dot{U}^{r}_{\rm o}(\lambda,t)\\\
\dot{V}^{r}_{\rm o}(\lambda,t)\end{bmatrix}=H_{\rm
o}^{r}(\lambda)\begin{bmatrix}U^{r}_{\rm o}(\lambda,t)\\\ V^{r}_{\rm
o}(\lambda,t)\end{bmatrix}\ ,$ (66)

with

$\displaystyle H_{\rm o}^{r}(\lambda)=\begin{bmatrix}-A_{\rm o}(\lambda)&D\\\
K_{\rm o}(\lambda)&A_{\rm o}^{T}(\lambda)&\end{bmatrix}\ ,$ (67)

and initial condition

$\displaystyle\begin{bmatrix}U^{r}_{\rm o}(\lambda,0)\\\ V^{r}_{\rm
o}(\lambda,0)\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I_{n+2}\\\
0\end{bmatrix}\ .$ (68)

As a consequence,

$\displaystyle\begin{bmatrix}U^{r}_{\rm o}(\lambda,t)\\\ V^{r}_{\rm
o}(\lambda,t)\end{bmatrix}=e^{H^{r}_{\rm
o}(\lambda)t}\begin{bmatrix}I_{n+2}\\\ 0\end{bmatrix}\ .$ (69)

The existence of the solution of Eq. (42a) for all $t^{\prime}\in[0,t]$
ensures that the corresponding matrix $U^{r}_{\rm o}(\lambda,t^{\prime})$ is
invertible in this interval. Conversely, if $U^{r}_{\rm
o}(\lambda,t^{\prime})$ is nonsingular for all $t^{\prime}\in[0,t]$, then
$C^{r}_{\rm o}(\lambda,t^{\prime})$ exists in the same interval and is given
by Eq. (65).

A similar transformation holds for the solution $C^{l}_{\rm o}(\lambda,t)$ of
Eq. (42b), with $U^{r}_{\rm o}(\lambda,t),V^{r}_{\rm o}(\lambda,t)$ replaced
by $U^{l}_{\rm o}(\lambda,t),V^{l}_{\rm o}(\lambda,t)$ and Eqs. (69) replaced
by

$\displaystyle\begin{bmatrix}U^{l}_{\rm o}(\lambda,t)\\\ V^{l}_{\rm
o}(\lambda,t)\end{bmatrix}=e^{H^{l}_{\rm
o}(\lambda)t}\begin{bmatrix}I_{n+2}\\\ C\end{bmatrix}\ .$ (70)

with

$\displaystyle H^{l}_{\rm o}(\lambda)=\begin{bmatrix}A_{\rm o}(\lambda)&D\\\
K_{\rm o}(\lambda)&-A_{\rm o}^{T}(\lambda)&\end{bmatrix}\ .$ (71)

The $2(n+2)\times 2(n+2)$ matrices $H^{r}_{\rm o}(\lambda)$ and $H^{l}_{\rm
o}(\lambda)$ are Hamiltonian matrices which play a central role in the
forthcoming analysis999We recall that a real Hamiltonian matrix $H$ of
dimension $2(n+2)\times 2(n+2)$ satisfies the equation
$\left[\begin{smallmatrix}0&-I_{n+2}\\\
I_{n+2}&0\end{smallmatrix}\right]H=-H^{T}\left[\begin{smallmatrix}0&-I_{n+2}\\\
I_{n+2}&0\end{smallmatrix}\right]$, where $I_{n+2}$ is the $(n+2)\times(n+2)$
identity matrix. The eigenvalues of $H$ then come in quadruples: if
$s\in\mathbb{C}$ is an eigenvalue, so are $\bar{s},-s$ and $-\bar{s}$, where
an overline indicates complex conjugation.. Their spectral properties are
investigated in Appendix C. Observe in particular that $H^{r}_{\rm
o}(\lambda)$ and $H^{l}_{\rm o}(\lambda)$ are invertible and that the
eigenvalue spectrum only depends on the anti-symmetric part of the matrices
$B_{\rm o}$ (a property that is shared by all linear current-type
observables). In the present case, the anti-symmetric part is the same for the
three observables [Eq. (20)], and owing to the fact that
$D_{ij}=\delta_{1i}\delta_{1j}$, the characteristic polynomial of the
Hamiltonian matrices is shown in Appendix C.0.1 to be given by

$\displaystyle(-1)^{n}p_{H}(\lambda,s)=p_{A}(s)p_{A}(-s)-\frac{2\lambda
g}{Q_{0}^{2}}(\frac{n}{\tau})^{2n}\>s[(1-\frac{s\tau}{n})^{n}-(1+\frac{s\tau}{n})^{n}]\
,$ (72)

where $p_{A}(s)$, the characteristic polynomial of $A$, is given by Eq. (24).

For the sake of simplicity, we will assume in the following that the matrices
$H^{r}_{\rm o}(\lambda)$ and $H^{l}_{\rm o}(\lambda)$ are diagonalizable. This
is not an essential assumption but it simplifies the presentation as we can
then avoid to deal with the more complicated Jordan forms of the Hamiltonian
matrices AFIJ2003 . More important is the fact that the values of $\lambda$
will be restricted to a certain open interval ${\cal
D}_{H}=(\lambda_{\min},\lambda_{\max})$ for which $H^{r}_{\rm o}(\lambda)$ and
$H^{l}_{\rm o}(\lambda)$ have no purely imaginary eigenvalues101010One may
have $\lambda_{\min}=-\infty$ but $\lambda_{\max}$ is positive and finite.. As
will be shown in Sec. IV.2.3, this restriction is justified by the fact that
the values of $\lambda$ outside ${\cal D}_{H}$ are irrelevant for the
determination of the rate function $I(a)$.

The condition $\lambda\in{\cal D}_{H}$ has two significant consequences:

$\bullet$ The Hamiltonian matrices are dichotomically separable AFIJ2003 ,
with $n+2$ eigenvalues
$s^{+}_{1}(\lambda),s^{+}_{2}(\lambda),...s^{+}_{n+2}(\lambda)$ having a
positive real part and $n+2$ eigenvalues
$s_{1}^{-}(\lambda)=-s_{1}^{+}(\lambda),s_{2}^{-}(\lambda)=-s_{2}^{+}(\lambda),...s_{n+2}^{-}(\lambda)=-s^{+}_{n+2}(\lambda)$
having a negative real part (the eigenvalues are counted with their
multiplicities and the ordering is arbitrary unless otherwise specified).

$\bullet$ Exact mathematical statements ensure that the corresponding
algebraic Riccati equations have real symmetric solutions (see the discussion
in the next section).

Let us consider the matrix $H_{\rm o}^{r}(\lambda)$. As a result of the above
assumptions, we can introduce a basis-change matrix $W^{r}_{\rm o}(\lambda)$
that diagonalizes $H_{\rm o}^{r}(\lambda)$ such that

$\displaystyle H_{\rm o}^{r}(\lambda)=W^{r}_{\rm o}(\lambda)\
\begin{bmatrix}J(\lambda)&0\\\ 0&-J(\lambda)&\end{bmatrix}[W^{r}_{\rm
o}(\lambda)]^{-1}\ ,$ (73)

where
$J(\lambda)=\mbox{diag}\big{(}s^{+}_{1}(\lambda),s^{+}_{2}(\lambda),...s^{+}_{n+2}(\lambda)\big{)}$.
Eq. (69) then becomes

$\displaystyle\begin{bmatrix}U^{r}_{\rm o}(\lambda,t)\\\ V^{r}_{\rm
o}(\lambda,t)\end{bmatrix}=W^{r}_{\rm
o}(\lambda)\begin{bmatrix}e^{J(\lambda)t}&0\\\
0&e^{-J(\lambda)t}&\end{bmatrix}[W^{r}_{\rm
o}(\lambda)]^{-1}\begin{bmatrix}I_{n+2}\\\ 0\end{bmatrix}\ .$ (74)

We next partition the matrix $W^{r}_{\rm o}(\lambda)$ into $4$ blocks of size
$(n+2)\times(n+2)$,

$\displaystyle W^{r}_{\rm o}(\lambda)=\begin{bmatrix}W^{r,11}_{\rm
o}(\lambda)&W^{r,12}_{\rm o}(\lambda)\\\ W^{r,21}_{\rm
o}(\lambda)&W^{r,22}_{\rm o}(\lambda)&\end{bmatrix}\ ,$ (75)

where the first $n+2$ columns are the eigenvectors relative to the eigenvalues
$s_{i}^{+}(\lambda)$ and the remaining columns are the eigenvectors relative
to the eigenvalues $s_{i}^{-}(\lambda)$.

Now, suppose that the submatrix $W^{r,22}_{\rm o}(\lambda)$ is nonsingular, so
that the matrix

$\displaystyle T^{r}_{\rm o}(\lambda)=$ $\displaystyle-[W^{r,22}_{\rm
o}(\lambda)]^{-1}W^{r,21}_{\rm o}(\lambda)$ (76)

exists. Simple manipulations, similar to those performed in Refs. K1973 ;
FJ1996  and detailed in Appendix D, then lead to a representation of the
solution of the RDE (42a) that only involves negative exponentials111111We
leave to the reader the proof that Eq. (77) does not depend on the choice of
the basis-change matrix $W^{r}_{\rm o}(\lambda)$.:

$\displaystyle C^{r}_{\rm o}(\lambda,t)$ $\displaystyle=[W^{r,21}_{\rm
o}(\lambda)+W^{r,22}_{\rm o}(\lambda)P^{r}_{\rm o}(\lambda,t)][W^{r,11}_{\rm
o}(\lambda)+W^{r,12}_{\rm o}(\lambda)P^{r}_{\rm o}(\lambda,t)]^{-1}\ ,$ (77)

where

$\displaystyle P^{r}_{\rm o}(\lambda,t)=e^{-J(\lambda)t}T^{r}_{\rm
o}(\lambda)e^{-J(\lambda)t}\ .$ (78)

Two important features of the solution are revealed by Eq. (77):

$\bullet$ $C^{r}_{\rm o}(\lambda,t)$ diverges at time $t$ if the matrix
$W^{r,11}_{\rm o}(\lambda)+W^{r,12}_{\rm o}(\lambda)P^{r}_{\rm o}(\lambda,t)$
is singular. Finite-time singularities in the solution are thus poles
corresponding to $\det[W^{r,11}_{\rm o}(\lambda)+W^{r,12}_{\rm
o}(\lambda)P^{r}_{\rm o}(\lambda,t)]=0$. (We recall that the generating
function $G_{\rm o,\lambda}(t)$ diverges before the solution or the RDE blows
up.)

$\bullet$ $C^{r}_{\rm o}(\lambda,t)$ converges towards the matrix
$\hat{C}^{r,+}_{\rm o}(\lambda)\equiv W^{r,21}_{\rm o}(\lambda)[W^{r,11}_{\rm
o}(\lambda)]^{-1}$ as $t\to\infty$. (This is true even if finite-time
singularities are present121212In fact, solving the linear system of ODEs (66)
instead of the original RDE is a practical method to bypass singularities,
which is needed in certain applications, in particular for boundary-value
problems AMR1988 ..)

As will be shown in the next subsection, $\hat{C}^{r,+}_{\rm o}(\lambda)$ is
the so-called “maximal” (real symmetric) solution of the corresponding
algebraic equation (86a) and its existence is certified LR1995  (which means
that the matrix $W^{r,11}_{\rm o}(\lambda)$ is nonsingular). On the other
hand, the assumption that $W^{r,22}_{\rm o}(\lambda)$ is nonsingular, which is
required for Eq. (77) to be meaningful, may not be satisfied for some values
of $\lambda$. In other words, the condition $\det W^{r,22}_{\rm
o}(\lambda)\neq 0$ ensures that the initial matrix $C^{r}_{\rm
o}(\lambda,0)=0$ belongs to the basin of attraction of $\hat{C}^{r,+}_{\rm
o}(\lambda)$. Otherwise, Eq. (77) is no longer valid and the solution of the
RDE (42a) goes to another limit or may fail to converge (i.e., oscillates), as
will be discussed in Sec. IV.3.2.

Likewise, the solution of Eq. (42b) can be represented as

$\displaystyle C^{l}_{\rm o}(\lambda,t)$ $\displaystyle=[W^{l,21}_{\rm
o}(\lambda)+W^{l,22}_{\rm o}(\lambda)P^{l}_{\rm o}(\lambda,t)][W^{l,11}_{\rm
o}(\lambda)+W^{l,12}_{\rm o}(\lambda)P^{l}_{\rm o}(\lambda,t)]^{-1}\ ,$ (79)

where

$\displaystyle P^{l}_{\rm o}(\lambda,t)=e^{-J(\lambda)t}T^{l}_{\rm
o}(\lambda)e^{-J(\lambda)t}\ ,$ (80)

and

$\displaystyle T^{l}_{\rm o}(\lambda)=$ $\displaystyle-[W^{l,22}_{\rm
o}(\lambda)-CW^{l,12}_{\rm o}(\lambda)]^{-1}[W^{l,21}_{\rm
o}(\lambda)-CW^{l,11}_{\rm o}(\lambda)]\ .$ (81)

This requires that the matrix $W^{l,22}_{\rm o}(\lambda)-CW^{l,12}_{\rm
o}(\lambda)$ is nonsingular. If true, $C^{l}_{\rm o}(\lambda,t)$ converges
asymptotically toward the matrix $\hat{C}^{l,+}_{\rm o}(\lambda)\equiv
W^{l,21}_{\rm o}(\lambda)[W^{l,11}_{\rm o}(\lambda)]^{-1}$ which is the
maximal (real symmetric) solution of the corresponding CARE and whose
existence is also guaranteed (implying that $W^{l,11}_{\rm o}(\lambda)$ is
nonsingular). Moreover, it is easily seen from the structure of the
Hamiltonian matrices $H_{\rm o}^{r}(\lambda)$ and $H_{\rm o}^{l}(\lambda)$
that

$\displaystyle\begin{bmatrix}W_{\rm o}^{l,11}&W_{\rm o}^{l,12}\\\ W_{\rm
o}^{l,21}&W_{\rm o}^{l,22}&\end{bmatrix}=\begin{bmatrix}-W_{\rm
o}^{r,12}&-W_{\rm o}^{r,11}\\\ W_{\rm o}^{r,22}&W_{\rm
o}^{r,21}&\end{bmatrix}\ .$ (82)

As a result, $\hat{C}^{l,+}_{\rm o}(\lambda)=-W^{r,22}_{\rm
o}(\lambda)[W^{r,12}_{\rm o}(\lambda)]^{-1}$.

For future reference, it is instructive to rewrite the conditions of
convergence towards $\hat{C}^{r,+}_{\rm o}(\lambda)$ and $\hat{C}^{l,+}_{\rm
o}(\lambda)$ as

$\displaystyle\det W^{r,22}_{\rm o}(\lambda)$ $\displaystyle=\det
W^{l,21}_{\rm o}(\lambda)=\det\hat{C}^{l,+}_{\rm o}(\lambda)\>\det
W^{l,11}_{\rm o}(\lambda)\neq 0\ ,$ (83)

and

$\displaystyle\det[W^{l,22}_{\rm o}(\lambda)-CW^{l,12}_{\rm
o}(\lambda)]=\det[W^{r,21}_{\rm o}(\lambda)+CW^{r,11}_{\rm
o}(\lambda)]=\det[\hat{C}^{r,+}_{\rm o}(\lambda)+C]\det W^{r,11}_{\rm
o}(\lambda)\neq 0\ .$ (84)

As $\det W^{l,11}_{\rm o}(\lambda)$ and $\det W^{r,11}_{\rm o}(\lambda)$ are
$\neq 0$, this implies the following equivalences:

$\displaystyle\det\hat{C}^{l,+}_{\rm o}(\lambda)\neq
0\>\>\Leftrightarrow\>\>C^{r}_{\rm o}(\lambda,t)\to\hat{C}_{\rm
o}^{r,+}(\lambda)$ (85a) $\displaystyle\det[\hat{C}^{r,+}_{\rm
o}(\lambda)+C]\neq 0\>\>\Leftrightarrow\>\>C^{l}_{\rm
o}(\lambda,t)\to\hat{C}_{\rm o}^{l,+}(\lambda)$ (85b)

These dual relations will be used again and again in the following. They are
one of the main reasons for which it is fruitful to study the generating
functions $G^{r}_{\rm o,\lambda}({\bf u}_{0},t)$ and $G^{l}_{\rm
o,\lambda}({\bf u},t)$ together.

Finally, we stress that the representations (77) and (79) of the solutions of
the RDEs (42) not only reveal the most significant features of the solutions
but are also useful for numerical calculations. Indeed, as shown in Appendix
C.0.2, we have explicit expressions of the basis-change matrices
$W^{r}_{w}(\lambda)$ and $W^{l}_{w}(\lambda)$ as a function of the eigenvalues
of the Hamiltonian matrices. We can then only consider the RDE corresponding
the operator ${\cal R}_{w,\lambda}$ and compute all matrices $C^{r}_{\rm
o}(\lambda,t)$ and $C^{l}_{\rm o}(\lambda,t)$ by changing the initial
conditions and using relations (56) and (58).

###  III.3 Fixed points of the Riccati flows

How does one know that the matrices $\hat{C}^{r,+}_{\rm o}(\lambda)$ and
$\hat{C}^{l,+}_{\rm o}(\lambda)$ exist and what are their properties? To
answer these questions, we now consider the stationary versions of the RDEs
(42),

$\displaystyle{\cal R}_{\rm o,\lambda}[\hat{C}^{r}_{\rm o}(\lambda)]=0$ (86a)
$\displaystyle{\cal R}_{\rm o,\lambda}[-\hat{C}^{l}_{\rm o}(\lambda)]=0\ ,$
(86b)

which are referred to as continuous-time algebraic Riccati equations (CAREs)
in the context of optimal control. These equations may have no solutions at
all or multiple solutions, including complex and non-symmetric ones, and we
first recall how these solutions, in particular real symmetric ones, can be
built.

It is known that there is a one-to-one correspondence between the solutions of
a CARE and certain invariant subspaces of the associated Hamiltonian matrix
(see Refs. BLW1991 ; LR1995 ; BIM2012  for reviews). Let us consider for
instance Eq. (86a) and denote a solution by $\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)$ (then $-\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ is a solution of
Eq. (86b)). A direct calculation yields

$\displaystyle H^{r}_{\rm o}(\lambda)\begin{bmatrix}I_{n+2}\\\
\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)&\end{bmatrix}=\begin{bmatrix}I_{n+2}\\\
\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)&\end{bmatrix}[-A_{\rm
o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)]\ ,$ (87)

which shows that the columns of the matrix
$\left[\begin{smallmatrix}I_{n+2}\\\ \hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)&\end{smallmatrix}\right]$ span a graph invariant subspace131313The
graph of a matrix $X$ is defined as the $(n+2)$-dimensional subspace
$\displaystyle G(X)=\mbox{Im}\begin{bmatrix}I_{n+2}\\\
X&\end{bmatrix}\in\mathbb{C}^{2(n+2)}$ where Im denotes the image or column
space of a matrix LR1995 . A subspace of $\mathbb{C}^{2(n+2)}$ is called a
graph subspace if it has the form $G(X)$ for some $X$. of $H^{r}_{\rm
o}(\lambda)$ and the eigenvalues of $-A_{\rm
o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ are a subset of the
eigenvalues of $H^{r}_{\rm o}(\lambda)$. This simple fact leads to the
following characterization of the solutions note20 : Each solution
$\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ corresponds to a set ${\cal
S}_{\lambda}^{(\alpha)}$ of $n+2$ eigenvalues
$s_{\alpha_{1}},s_{\alpha_{2}},...s_{\alpha_{n+2}}$ of $H^{r}_{\rm
o}(\lambda)$ [specifically, the eigenvalues of $-A_{\rm
o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$] and $n+2$ associated
eigenvectors ${\bf e}^{r}_{\alpha_{j}}=\left[\begin{smallmatrix}{\bf
y}^{r}_{\alpha_{j}}\\\ {\bf z}^{r}_{\alpha_{j}}\end{smallmatrix}\right]$ (with
${\bf y}^{r}_{\alpha_{j}},{\bf z}^{r}_{\alpha_{j}}\in\mathbb{C}^{n+2}$), such
that

$\displaystyle\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)=Z^{r,(\alpha)}_{\rm
o}(\lambda)[Y^{r,(\alpha)}_{\rm o}(\lambda)]^{-1}\ ,$ (88)

where $Y^{r,(\alpha)}_{\rm o}(\lambda)=[{\bf y}^{r}_{\alpha_{1}},{\bf
y}^{r}_{\alpha_{2}},...{\bf y}^{r}_{\alpha_{n+2}}]$ and $Z^{r,(\alpha)}_{\rm
o}(\lambda)=[{\bf z}^{r}_{\alpha_{1}},{\bf z}^{r}_{\alpha_{2}},...{\bf
z}^{r}_{\alpha_{n+2}}]$. The invertibility of the $(n+2)\times(n+2)$ matrix
$Y^{r,(\alpha)}_{\rm o}(\lambda)$ is the condition ensuring that the solution
$\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ exists, and conversely. The solutions
of Eq. (86b) can be characterized in the same manner, with the index $r$
replaced by $l$ and $A_{\rm o}(\lambda)$ replaced by $-A_{\rm o}(\lambda)$ in
Eq. (87). (Recall that $H_{\rm o}^{r}(\lambda)$ and $H_{\rm o}^{r}(\lambda)$
have the same eigenvalue spectrum so that $\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)$ and $\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$ correspond to the
same subset ${\cal S}_{\lambda}^{(\alpha)}$ of eigenvalues of $H_{\rm
o}^{r}(\lambda)$ and $H_{\rm o}^{l}(\lambda)$.)

Interestingly, the solutions of Eqs. (86) for an observable $\rm o$’ can be
readily obtained from the solutions for the observable $\rm o$. This results
from the invariance property (54) of the Riccati operator ${\cal R}_{\rm
o,\lambda}$, which yields

$\displaystyle\hat{C}^{r,(\alpha)}_{\rm
o^{\prime}}(\lambda)=\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)+\lambda(S_{\rm
o}-S_{\rm o^{\prime}})=\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)+\lambda(B_{\rm
o}-B_{\rm o^{\prime}})$ (89a) $\displaystyle\hat{C}^{l,(\alpha)}_{\rm
o^{\prime}}(\lambda)=\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)-\lambda(S_{\rm
o}-S_{\rm o^{\prime}})=\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)-\lambda(B_{\rm
o}-B_{\rm o^{\prime}})\ .$ (89b)

As a result,

$\displaystyle\hat{C}^{r,(\alpha)}_{\rm
o^{\prime}}(\lambda)+\hat{C}^{l,(\alpha)}_{\rm
o^{\prime}}(\lambda)=\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)+\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)\ .$ (90)

Note that $\hat{C}^{r,(\alpha)}_{\rm o^{\prime}}(\lambda)$ and
$\hat{C}^{l,(\alpha)}_{\rm o^{\prime}}(\lambda)$ correspond to the same set
${\cal S}_{\lambda}^{(\alpha)}$ of eigenvalues as $\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)$ and $\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$, which is the set of
eigenvalues of the matrices $-A_{\rm o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)$ and $A_{\rm o}(\lambda)+D\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$.
Indeed, from the definition of $A_{\rm o}(\lambda)$ [Eq. (45)], one has

$\displaystyle-A_{\rm o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$
$\displaystyle=-A+\lambda DB_{\rm o}+D[\hat{C}^{r,(\alpha)}_{\rm
o^{\prime}}(\lambda)-\lambda(B_{\rm o}-B_{\rm o^{\prime}}]=-A_{\rm
o^{\prime}}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o^{\prime}}(\lambda)$ (91a)
$\displaystyle A_{\rm o}(\lambda)+D\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$
$\displaystyle=A-\lambda DB_{\rm o}+D[\hat{C}^{l,(\alpha)}_{\rm
o^{\prime}}(\lambda)+\lambda(B_{\rm o}-B_{\rm o^{\prime}})]=A_{\rm
o^{\prime}}(\lambda)+D\hat{C}^{l,(\alpha)}_{\rm o^{\prime}}(\lambda)\ .$ (91b)

Another interesting consequence of Eqs. (89) is that

$\displaystyle\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$
$\displaystyle=\hat{C}^{r,(\alpha)}_{\rm o,antisym}(\lambda)-\lambda B_{\rm
o,sym}$ (92a) $\displaystyle\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$
$\displaystyle=\hat{C}^{l,(\alpha)}_{\rm o,antisym}(\lambda)+\lambda B_{\rm
o,sym}\ ,$ (92b)

where $B_{\rm o,sym}\equiv(B_{\rm o}+B^{T}_{\rm o})/2$ is the symmetric part
of the matrix $B_{\rm o}$ and $\hat{C}^{r,(\alpha)}_{\rm o,antisym}(\lambda)$
and $\hat{C}^{l,(\alpha)}_{\rm o,antisym}(\lambda)$ are the solutions of the
CAREs (86) associated with $B_{\rm o,antisym}\equiv(B_{\rm o}-B^{T}_{\rm
o})/2$. Consequently,

$\displaystyle\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)+\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)=\hat{C}^{r,(\alpha)}_{\rm
o,antisym}(\lambda)+\hat{C}^{l,(\alpha)}_{\rm o,antisym}(\lambda)\ .$ (93)

In accordance with Eqs. (41), we associate to each solutions
$\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ and $\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)$ of the CAREs (86) the scalar functions

$\displaystyle f^{r,(\alpha)}_{\rm o}(\lambda)$
$\displaystyle=-\frac{1}{2}\mbox{Tr}[D(\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)+\lambda B_{\rm o})]$ (94a) $\displaystyle f^{l,(\alpha)}_{\rm
o}(\lambda)$ $\displaystyle=-\frac{1}{2}\mbox{Tr}[D(\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)-\lambda B_{\rm o})]-\mbox{Tr}(A)\ .$ (94b)

These two functions are actually equal and independent of the observables.
This follows from the fact that ${\cal
S}_{\lambda}^{(\alpha)}=\\{s_{\alpha_{i}}\\}_{i=1}^{n+2}$ is the spectrum of
both $-A_{\rm o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ and $A_{\rm
o}(\lambda)+D\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$, whatever $\rm o$, as we
just noticed. As a consequence,

$\displaystyle\sum_{i=1}^{n+2}s_{\alpha_{i}}(\lambda)$
$\displaystyle=\mbox{Tr}[-A_{\rm o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)]=\mbox{Tr}[A_{\rm o}(\lambda)+D\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)]\ .$ (95)

Using the definition of $A_{\rm o}(\lambda)$ [Eq. (45)], this yields
$f^{r,(\alpha)}_{\rm o}(\lambda)=f^{l,(\alpha)}_{\rm
o}(\lambda)=f^{(\alpha)}(\lambda)$ with

$\displaystyle
f^{(\alpha)}(\lambda)=-\frac{1}{2}[\mbox{Tr}(A)+\sum_{i=1}^{n+2}s_{\alpha_{i}}(\lambda)]\
.$ (96)

So far, we have considered all solutions of Eq. (86) (assuming that they
exist). However, in the present context we are only interested in real
symmetric solutions. It turns out that the existence of such solutions is
ensured due to the following two properties of the matrices involved in the
Riccati operator ${\cal R}_{\rm o,\lambda}$ defined by Eq. (44) note6 :

i) the matrix $D$ is positive semidefinite,

ii) the pair of matrices $(A_{\rm o}(\lambda),D)$ is controllable141414A
standard result in control theory is that a pair $(A,B)$, where $A$ is a
$n\times n$ matrix and $B$ is a $n\times m$ matrix, is controllable if the
rank of the $n\times nm$ matrix $[B,AB,A^{2}B,...A^{n-1}B]$ is equal to $n$
B2021 . It is easy to show that the pair $(A_{\rm o}(\lambda),D)$ satisfies
this crucial property by using the equivalent PBH (Popov-Belevitch-Hautus)
test which states that there must be no nonzero left eigenvector ${\bf v}^{T}$
of $A_{\rm o}(\lambda)$ such that ${\bf v}^{T}\cdot D=0$. Since
$D_{ij}=\delta_{i1}\delta_{j1}$, the second condition imposes that ${\bf
v}^{T}=(0,v_{2},v_{3},...v_{n+2})^{T}$. The first condition then yields
$v_{2}=v_{3}...=v_{n+2}=0$, as can be readily checked by inspection..

These properties also ensure that note7  :

$\bullet$ Each real symmetric solution $\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$
or $\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)$ corresponds to a set ${\cal
S}^{(\alpha)}_{\lambda}$ of eigenvalues of $H^{r}_{\rm o}(\lambda)$ and
$H^{l}_{\rm o}(\lambda)$ such that $s\in{\cal S}^{(\alpha)}_{\lambda}$ implies
$\bar{s}\in{\cal S}^{(\alpha)}_{\lambda}$ and $-\bar{s}\notin{\cal
S}^{(\alpha)}_{\lambda}$.

$\bullet$ The matrices $\hat{C}^{r,+}_{\rm o}(\lambda)=W^{r,21}_{\rm
o}(\lambda)[W^{r,11}_{\rm o}(\lambda)]^{-1}$ and $\hat{C}^{l,+}_{\rm
o}(\lambda)=W^{l,21}_{\rm o}(\lambda)[W^{l,11}_{\rm o}(\lambda)]^{-1}$
introduced previously exist and are the maximal real symmetric solutions of
Eqs. (86a) and (86b) with respect to the positive definiteness ordering. They
are obtained by taking ${\cal S}^{(\alpha)}_{\lambda}$ to be the set ${\cal
S}_{\lambda}^{+}=\\{s^{+}_{1},s^{+}_{2},...s^{+}_{n+2}\\}$ of eigenvalues with
a positive real part (see the partitioning of the basis change matrices
$W^{r}_{\rm o}(\lambda)$ and $W^{l}_{\rm o}(\lambda)$). Likewise, the matrices
$C^{r,-}_{\rm o}(\lambda)=W_{\rm o}^{r,22}(\lambda)[W_{\rm
o}^{r,12}(\lambda)]^{-1}$ and $C^{l,-}_{\rm o}(\lambda)=W_{\rm
o}^{l,22}(\lambda)[W_{\rm o}^{l,12}(\lambda)]^{-1}$ exist and are the minimal
solutions obtained from the set ${\cal
S}_{\lambda}^{-}=\\{s^{-}_{1},s^{-}_{2},...s^{-}_{n+2}\\}$ of eigenvalues with
a negative real part. This means that all other real symmetric solutions
$\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ (resp. $\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)$) are such that $\hat{C}^{r,-}_{\rm
o}(\lambda)\leq\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)\leq\hat{C}^{r,+}_{\rm
o}(\lambda)$ (resp. $\hat{C}^{l,-}_{\rm
o}(\lambda)\leq\hat{C}^{l,(\alpha)}_{\rm o}(\lambda)\leq\hat{C}^{l,+}_{\rm
o}(\lambda)$.

The maximal solutions $\hat{C}^{r,+}_{\rm o}(\lambda)$ and $\hat{C}^{l,+}_{\rm
o}(\lambda)$ can be shown to be analytic functions of $\lambda$ RR1988 , but
at variance with the common situation in LQ optimal control LR1995  these
matrices are not necessarily positive semidefinite. On the other hand, since

$\displaystyle\hat{C}^{l,+}_{\rm o}(\lambda)$
$\displaystyle=-\hat{C}^{r,-}_{\rm o}(\lambda)$
$\displaystyle\hat{C}^{l,-}_{\rm o}(\lambda)$
$\displaystyle=-\hat{C}^{r,+}_{\rm o}(\lambda)$ (97)

from the relation between the basis change matrices $W_{\rm o}^{r}(\lambda)$
and $W_{\rm o}^{l}(\lambda)$ [Eq. (82)], we have

$\displaystyle\hat{C}^{r,+}_{\rm o}(\lambda)+\hat{C}^{l,+}_{\rm
o}(\lambda)=\hat{C}^{r,+}_{\rm o}(\lambda)-\hat{C}^{r,-}_{\rm o}(\lambda)>0\
.$ (98)

As shown in Appendix E, $\hat{C}^{r,+}_{\rm o}(\lambda)$ and
$\hat{C}^{l,+}_{\rm o}(\lambda)$ have another important property: They are the
only solutions of the CAREs (86) that satisfy $\hat{C}^{r,+}_{\rm o}(0)=0$ and
$\hat{C}^{l,+}_{\rm o}(0)=C$. Since $G^{r}_{\rm o,\lambda=0}({\bf u}_{0},t)=1$
and $G^{l}_{\rm o,\lambda=0}({\bf u},t)=p({\bf u})$ for all ${\bf u}_{0},{\bf
u}$ and $t$ from the definition of the generating functions [Eqs. (33) and
(34)], these conditions must be obeyed by the solutions of the RDEs (42a) and
(42b), respectively151515It is suggested in Ref. DBT2023  that the behavior at
$\lambda=0$ can be used to select the correct asymptotic fixed point of the
differential equation among all the solutions of the corresponding CARE.
However, this is not a viable procedure in general as it assumes that these
solutions are known analytically, which is only true for simple two-
dimensional systems. By showing that the maximal solution is the one that
satisfies the exact condition at $\lambda=0$ we overcome this obstacle.
Furthermore, we must keep in mind that the solution of the RDE may not be a
continuous function of $\lambda$, as will be seen in Sec. IV.3..

Finally, we introduce the function $f^{+}(\lambda)$ associated with the
maximal solutions of the CAREs which will play a prominent role in the
following. Since the maximal solutions correspond to the set ${\cal
S}_{\lambda}^{+}$, we have

$\displaystyle
f^{+}(\lambda)=-\frac{1}{2}[\mbox{Tr}(A)+\sum_{i=1}^{n+2}s^{+}_{i}(\lambda)]$
(99)

from Eq. (96). $f^{+}(\lambda)$ stands out among the functions
$f^{(\alpha)}(\lambda)$ associated with the other solutions of Eqs. (86)
because of the following three properties:

i) it is an analytic function of $\lambda$ in the interval ${\cal D}_{H}$
(since $\hat{C}^{r,+}_{\rm o}(\lambda)$ and $\hat{C}^{l,+}_{\rm o}(\lambda)$
are analytic);

ii) it satisfies

$\displaystyle f^{+}(0)=0\ ,$ (100)

as can be readily seen by inserting $\hat{C}^{r,+}_{\rm o}(0)=0$ into Eq.
(94a)161616Alternatively, one can use the fact that
$p_{H}(0,s)=(-1)^{n}p_{A}(s)p_{A}(-s)$ from Eq. (72). Therefore, the zeros of
$p_{H}(0,s)$ with a positive real part are the zeros of $p_{A}(-s)$ and
$\sum_{i=1}^{n+2}s^{+}_{i}(0)=-\mbox{Tr}(A)$.;

iii) it obeys the inequality

$\displaystyle f^{(\alpha)}(\lambda)>f^{+}(\lambda),\;{\rm for}\;{\cal
S}^{(\alpha)}_{\lambda}\neq{\cal S}^{+}_{\lambda}.$ (101)

Indeed, a set ${\cal S}_{\lambda}^{(\alpha)}\neq{\cal S}_{\lambda}^{+}$
contains $m^{(\alpha)}\neq 0$ eigenvalues $s^{-}_{\alpha_{i}}$ with a negative
real part and $n+2-m^{(\alpha)}$ eigenvalues $s^{+}_{\alpha_{i}}$ with a
positive real part. In consequence,

$\displaystyle\sum_{i=1}^{n+2}s_{\alpha_{i}}(\lambda)$
$\displaystyle=\sum_{i=1}^{m^{(\alpha)}}s^{-}_{\alpha_{i}}(\lambda)+\sum_{i=m^{(\alpha)}+1}^{n+2-m^{(\alpha)}}s^{+}_{\alpha_{i}}(\lambda)=\sum_{i=1}^{n+2}s^{+}_{i}(\lambda)-2\sum_{i=1}^{m^{(\alpha)}}s^{-}_{\alpha_{i}}(\lambda)\
,$ (102)

where we have used the symmetry of the eigenvalues with respect to the
imaginary axis. Inserting the above into Eq. (96) and using Eq. (99), we then
obtain

$\displaystyle
f^{(\alpha)}(\lambda)-f^{+}(\lambda)=-\sum_{i=1}^{m^{(\alpha)}}s^{-}_{\alpha_{i}}(\lambda)>0\
.$ (103)

An alternative expression of $f^{+}(\lambda)$ in terms of the spectral density
of the process will be derived in Sec. IV.2.2.

##  IV Time evolution of the moment generating functions

Equipped with the explicit representations of the time-dependent solutions
$C^{r}_{\rm o}(\lambda,t)$ and $C^{l}_{\rm o}(\lambda,t)$ of the RDEs (42)
[Eqs. (77) and (79)], with the definition of the maximal solutions
$\hat{C}^{r,+}_{\rm o}(\lambda)$ and $\hat{C}^{l,+}_{\rm o}(\lambda)$ of the
CAREs (86), and with the expression of the function $f^{+}(\lambda)$, we are
now in position to study the time evolution of the moment generating function
$G_{\rm o,\lambda}(t)$ given by Eqs. (48) and in particular their long-time
behavior. We will then explicitly compute the SCGF $\mu_{\rm o}(\lambda)$ and
the sub-exponential prefactors $g_{\rm o}(\lambda)$. We recall that the values
of $\lambda$ are restricted to the interval ${\cal
D}_{H}=(\lambda_{\min},\lambda_{\max})$ for which the Hamiltonian matrices
have no purely imaginary eigenvalues and the existence of real symmetric
solutions of the CAREs is ensured171717This does not mean that real symmetric
solutions of the RDEs do not exist for $\lambda\not\in{\cal D}_{H}$, but they
blow up at some finite time. Note also that the extremal solutions of the
CAREs may also exist when the Hamiltonian matrices have purely imaginary
eigenvalues. However, this requires that the partial multiplicities of these
eigenvalues are all even note30 , and one can show that it is not the case for
the model under study..

To make it more concrete, the forthcoming discussion is illustrated by
numerical results obtained for a gamma-distributed delay with $n=5$ (which
generates a system of $n+2=7$ dynamical variables). Although this is still far
from the case of a discrete delay, the gamma distribution (5) is already
markedly peaked around $t=\tau$ and further investigations show that the
qualitative behavior of the fluctuations does not change significantly for
larger values of $n$. We choose (rather arbitrarily) $Q_{0}=2$, so that the
harmonic oscillator is in a moderate underdamped regime, and a feedback gain
$g=1.5$. We then vary $\tau$ in order to illustrate different types of
behavior. For these values of the parameters, the system reaches a stable
stationary state for all values of $\tau$. The absence of
destabilizing/stabilizing Hopf bifurcations such as those shown in Fig. 1 is
irrelevant to our discussion.

###  IV.1 Finite-time divergences

As mentioned in Sec. III.1, the generating function $G_{\rm o,\lambda}(t)$
given by Eqs. (48) may diverge in a finite time, depending on the observable
$\rm o$ and the value of $\lambda$. In some cases one can prove from the
outset that such divergence does not take place [see for instance Eqs. (63)
and (64)]. In general, however, one needs to numerically compute $C^{r}_{\rm
o}(\lambda,t)$ from Eq. (77) or $C^{l}_{\rm o}(\lambda,t)$ from Eq. (79). By
increasing (resp. decreasing) $\lambda$ from $0$ to $\lambda_{\max}$ (resp.
from $0$ to $\lambda_{\min}$) at a given $t$, one can then determine the first
value of $\lambda$ for which $\det(C^{r}_{\rm o}(\lambda,t)+C)=\det C^{l}_{\rm
o}(\lambda,t)=0$ (we recall that the two determinants vanish simultaneously
since the expressions (48a) and (48b) of $G_{\rm o,\lambda}(t)$ are consistent
as long as the matrices $C^{r}_{\rm o}(\lambda,t)$ and $C^{l}_{\rm
o}(\lambda,t)$ exist). This defines a time-dependent interval ${\cal D}_{\rm
o}(t)=\big{(}\lambda_{\rm o}^{-}(t),\lambda_{\rm o}^{+}(t)\big{)}$ such that
$C+C^{r}_{\rm o}(\lambda,t)>0$ and $C^{l}_{\rm o}(\lambda,t)>0$ for all
$\lambda\in{\cal D}_{\rm o}(t)$. Accordingly, $G_{\rm o,\lambda}(t)$ is finite
for $\lambda\in{\cal D}_{\rm o}(t)$ and diverges at $\lambda=\lambda_{\rm
o}^{-}(t)$ and $\lambda=\lambda_{\rm o}^{+}(t)$. Since the very definition of
the SCGF requires that $\lim_{t\to\infty}(1/t)\ln G_{\rm
o,\lambda}(t)<\infty$, the domain of existence of the SCGF is ${\cal D}_{\rm
o}(\infty)$.

Figure 2:  (Color on line) (a) Plots of $\det C^{l}_{w}(\lambda,t)/\det C$
(solid lines) and $\det(C^{r}_{w}(\lambda,t)+C)/\det C$ (dashed lines) versus
$\lambda$ for $t=3,5,10$ and $\tau=1$. (b) Corresponding moment generating
function $G_{w,\lambda}(t)$. Figure 3:  (Color on line) (a) Plots of $\det
C^{l}_{w}(\lambda,t)/\det C$ (solid lines) and $\det(C^{r}_{\rm
o}(\lambda,t)+C)/\det C$ (dashed lines) versus $\lambda$ for $t=3,5,10$ and
$\tau=3$. (b) Corresponding moment generating function $G_{w,\lambda}(t)$.

Two typical examples of the evolution of the determinants with $\lambda$ at
different times are shown in Figs. 2 and 3 together with the corresponding
moment generating functions. The calculations are performed for the observable
${\cal W}_{t}$ but the same type of results are obtained for the other
observables. Inspection of the characteristic polynomial $p_{H}(s)$ reveals
that $\lambda_{\min}=-472.080$ and $\lambda_{\max}\simeq 1.019$ for $\tau=1$
and $\lambda_{\min}\simeq-5.118$ and $\lambda_{\max}\simeq 1.665$ for
$\tau=3$. Moreover, it is found that $T_{x}/T\approx 1.663$, $T_{v}/T\approx
0.570$ for $\tau=1$ and $T_{x}/T\approx 1.152$, $T_{v}/T\approx 0.857$ for
$\tau=3$. We can thus infer from Eq. (63) that the matrices
$C^{r}_{w}(\lambda,t)+C$ and $C^{l}_{w}(\lambda,t)$ are always positive
definite and the generating function $G_{w,\lambda}(t)$ is always finite for
values of $\lambda$ in the range $0\leq\lambda\lesssim 0.601$ for $\tau=1$ and
$0\leq\lambda\lesssim 0.868$ for $\tau=3$.

By comparing Figs. 2 and 3, we can see at once that the width of the interval
${\cal D}_{w}(t)$ does not decrease monotonically with time for $\tau=3$, in
contrast to the case $\tau=1$. (Note in passing that $G_{w,\lambda}(t)$ does
not satisfy the symmetry $\lambda\leftrightarrow 1-\lambda$; see the
discussion in Ref. RTM2017 .) The difference between these two cases is even
more manifest if we plot the evolution of the determinants with $t$ at fixed
$\lambda$, as done in Fig. 4 (for brevity, we only show the behavior of $\det
C^{l}_{w}(\lambda,t)$ and focus on negative values of $\lambda$).

Figure 4:  (Color on line) Time evolution of $\det C^{l}_{w}(\lambda,t)/\det
C$ for $\tau=1$ (a) and $\tau=3$ (b). The values of $\lambda$ decrease from
top to bottom and $\det C^{l}_{w}(\lambda,t)$ is always positive for
$\lambda\gtrsim-0.963$ (red curve in panel a) or $\lambda\gtrsim-1.51$ (red
curve in panel b). The parts of the curves corresponding to $\det
C^{l}_{w}(\lambda,t)<0$ are plotted as dashed lines.

In both cases, there is a critical value of $\lambda$ above which $\det
C^{l}_{w}(\lambda,t)>0$ at all times and the generating function is always
finite. For smaller values of $\lambda$, Figs. 4(a) and 4(b) depict two
different scenarios. For $\tau=1$, the determinant vanishes at a unique time
$t^{-}_{w}(\lambda)$ which increases monotonically to infinity as $\lambda$
approaches the critical value. In this case $t^{-}_{w}(\lambda)$ is just the
inverse function of $\lambda_{w}^{-}(t)$. For $\tau=3$, the curves cross the
$t$ axis a number of times (as the smallest eigenvalue of the matrix changes
its sign), so that the determinant is positive in a certain time range, then
negative, then positive again, etc. The function $t^{-}_{w}(\lambda)$ is thus
multivalued and the fact that the positive parts of the curves disappear at
different times explains the non-monotonic behavior observed in Fig. 3 as $t$
varies.

A similar behavior is observed for $\lambda>0$ and this eventually leads to
Figs. 5(a) and 5(b) which describe how the domain of existence of
$G_{w,\lambda}(t)$ evolves with time. Note that $\lambda_{w}^{+}(t)>0.601$ for
$\tau=1$ and $\lambda_{w}^{+}(t)>0.868$ for $\tau=3$, in agreement with the
predictions of Eq. (63).

Figure 5:  (Color on line) Evolution of the domain of existence of the
generating function $G_{w,\lambda}(t)$, the interval ${\cal
D}_{w}(t)=\big{(}\lambda^{+}_{w}(t),\lambda^{-}_{w}(t)\big{)}$, for $\tau=1$
(a) and $\tau=3$ (b). The matrices $C+C^{r}_{w}(\lambda,t)$ and
$C_{w}^{l}(\lambda,t)$ are positive definite and $G_{w,\lambda}(t)$ is finite
in the unshaded regions. At fixed $t$, $G_{w,\lambda}(t)$ diverges as
$\lambda\to\lambda^{+}_{w}(t)$ and $\lambda\to\lambda^{-}_{w}(t)$ (solid red
lines). The black dashed lines indicate the limits $\lambda_{w1}$ and
$\lambda_{w2}$ of the interval $\hat{\cal D}_{w}$ (see Fig. 8 below).

We stress that the generating function $G_{w,\lambda}(t)$ is finite at time
$t$ if $\det C_{w}^{l}(\lambda,t)>0$ whatever the sign of the determinant for
$t^{\prime}<t$. Otherwise, one would not obtain the same values of the
thresholds $\lambda_{w}^{-}(t)$ and $\lambda_{w}^{+}(t)$ by varying $\lambda$
at fixed $t$ or varying $t$ at fixed $\lambda$, and the results displayed in
Fig. 4(b) would not be consistent with those in Fig. 3. We thus disagree with
the statement made by C. Kwon et al in Ref. KNP2011  that the generating
function181818We recall that our expression (48b) of the generating function
is similar to the one derived in Ref. KNP2011 . is well defined only when the
determinant is positive for all $t^{\prime}<t$, which implies that the
threshold is “frozen” beyond a critical time191919We are indebted to Marco
Zamparo for providing us with a simple example illustrating that this
assertion is erroneous.. For $\tau=3$, this would amount to rejecting the non-
monotonic dependence of $\lambda^{+}_{w}(t)$ and $\lambda^{-}_{w}(t)$ on time
displayed in Fig. 5(b). We therefore dispute the claim that the generating
function jumps discontinuously to infinity beyond this critical time (see Fig.
2(b) in Ref. KNP2011 ) and we challenge the existence of the so-called
“dynamic phase transitions” discussed in Refs. KNP2011 ; NKP2013 ; KNP2013 .

In fact, the behavior of $\lambda_{w}^{-}(t)$ and $\lambda_{w}^{+}(t)$ with
$t$ can be directly checked by looking at the pdf $P_{t}(\beta{\cal
W}_{t}={\cal W})$. Indeed, the divergences of $G_{w,\lambda}(t)$ at
$\lambda^{\pm}_{w}(t)$ signal that $P_{t}({\cal W})$ has exponential tails.
More precisely, the numerical solution of the RDE (42a) reveals that $\det
C_{w}^{l}(\lambda,t)\sim\det(C_{w}^{r}(\lambda,t)+C)\sim(\lambda-\lambda^{\pm}_{w}(t))$
as $\lambda\to\lambda^{\pm}_{w}(t)$, so that, from Eqs. (48),

$\displaystyle G_{w,\lambda}(t)\sim(\lambda-\lambda^{\pm}_{w}(t))^{-1/2}$
(104)

and, thus KNP2011 ,

$\displaystyle P_{t}({\cal W})\sim\frac{e^{\lambda^{\mp}_{w}(t){\cal
W}}}{\sqrt{|{\cal W}|}}\ ,\>\>{\cal W}\to\pm\infty\ .$ (105)

Figure 6:  (Color on line) The pdf $P_{t}({\cal W})$ versus ${\cal W}$ for
$\tau=3$ at different times: $t=3$ (blue stars), $5$ (red circles), and $10$
(black squares). $P_{t}({\cal W})$ is obtained from numerical simulations of
$N=5.10^{6}$ samples. The straight lines on the right-hand side correspond to
the asymptotic form (105) with the numerical values of $\lambda_{w}^{-}(t)$
taken from Fig. 5(b).

In Fig. 6, the pdf obtained by numerically integrating the set of dynamical
equations (II.1) for $n=5$ and $\tau=3$ is plotted in the semi-log scale. It
is manifest that the slopes of the tails of the pdf do not vary monotonically
with time and they are very well described by the asymptotic form (105) (see
the solid lines on the right-hand side).

Figure 7: (Color on line) (a) Time evolution of the variance $\langle
W_{t}^{2}\rangle_{c}$ for $\tau=1$ and $\tau=3$. The circles represent the
results of the numerical simulation. (b) The width
$l_{w}(t)=\lambda_{w}^{+}(t)-\lambda_{w}^{-}(t)$ of the interval ${\cal
D}_{w}(t)$ for $\tau=3$ (red dashed line) is compared to $[\langle
W_{t}^{2}\rangle_{c}-\mu_{w}^{\prime\prime}(0)t]^{-1}$ (solid black line).

Remarkably, the nontrivial behavior of $\lambda_{w}^{\pm}(t)$ for $\tau=3$
also affects the variance $\langle W^{2}_{t}\rangle_{c}=\langle
W_{t}^{2}\rangle-\langle W_{t}\rangle^{2}=\partial^{2}\ln
G_{w,\lambda}(t)/\partial^{2}\lambda|_{\lambda=0}$. As shown in Fig. 7, the
non-monotonic variation of $\langle W^{2}_{t}\rangle_{c}$ for $\tau=3$ is
related to the evolution of the width
$l_{w}(t)=\lambda_{w}^{+}(t)-\lambda_{w}^{-}(t)$ of the interval ${\cal
D}_{w}(t)$. Since the variance is larger when $l_{w}(t)$ is smaller and
$\langle W_{t}^{2}\rangle_{c}\sim\mu_{w}^{\prime\prime}(0)t$ as $t\to\infty$
(where $\mu_{w}(\lambda)=f^{+}(\lambda)$ is the SCGF - see Sec. IV.2.3 below),
this relation is made more visible by comparing $l_{w}(t)$ to the inverse of
$\langle W_{t}^{2}\rangle_{c}-\mu_{w}^{\prime\prime}(0)t$, as done in Fig.
7(b). So the pdf $P_{t}({\cal W})$ does not necessarily become more
distributed as $t$ increases, as could be naively expected.

###  IV.2 Long-time behavior

####  IV.2.1 Domain of existence of the scaled-cumulant generating function
(SCGF)

We now focus on the behavior for $t\to\infty$. As can be seen in Fig. 5, the
thresholds $\lambda_{w}^{\pm}(t)$ converge to finite values
$\lambda_{w}^{\pm}(\infty)$. More generally, for an observable $\rm o$, ${\cal
D}_{\rm o}(t)\to{\cal D}_{\rm o}(\infty)=\big{(}\lambda_{\rm
o}^{-}(\infty),\lambda_{\rm o}^{+}(\infty)\big{)}$ with ${\cal D}_{\rm
o}(\infty)\subseteq{\cal D}_{H}$. On the other hand, the analysis of Sec.
III.2.2 tells us that $C^{r}_{\rm o}(\lambda,t)\to\hat{C}^{r,+}(\lambda)$ and
$C^{l}_{\rm o}(\lambda,t)\to\hat{C}^{l,+}(\lambda)$ for generic values of
$\lambda\in{\cal D}_{H}$ (i.e., for values of $\lambda$ such that
$\det\hat{C}^{l,+}_{\rm o}(\lambda)\neq 0$ and $\det(\hat{C}^{r,+}_{\rm
o}(\lambda)+C)\neq 0$). Therefore, ${\cal D}_{\rm o}(\infty)$, the domain of
existence of the SCGF, coincides with the interval $\hat{\cal D}_{\rm
o}=(\lambda_{\rm o1},\lambda_{\rm o2})$ for which the matrices
$\hat{C}^{r,+}_{\rm o}(\lambda)+C$ and $\hat{C}^{l,+}_{\rm o}(\lambda)$ are
both positive definite. Generically,

${\cal D}_{\rm o}(\infty)=\hat{\cal D}_{\rm o}\subseteq{\cal D}_{H},$ (106)

with $\hat{\cal D}_{\rm o}={\cal D}_{H}$ when these matrices
$\hat{C}^{r,+}_{\rm o}(\lambda)+C$ and $\hat{C}^{l,+}_{\rm o}(\lambda)$ are
positive definite for all $\lambda\in{\cal D}_{H}$.

Figure 8:  (Color on line) $\det\hat{C}^{l,+}_{w}(\lambda)/\det C$ (solid
lines) and $\det(\hat{C}^{r,+}_{w}(\lambda)+C)/\det C$ (dashed lines) as a
function of $\lambda$ for $\tau=1$ (a) and $\tau=3$ (b) (in both panels, the
largest value of $\lambda$ is $\lambda_{\max}$). The arrows indicate the range
$\hat{\cal D}_{w}$ of values of $\lambda$ for which the matrices
$\hat{C}^{r,+}_{w}(\lambda)+C$ and $\hat{C}^{l,+}_{w}(\lambda)$ are both
positive definite.

This is verified numerically in Fig. 8 which shows the variation of the
determinants of $\hat{C}^{r,+}_{w}(\lambda)+C$ and
$\hat{C}^{l,+}_{w}(\lambda)$ with $\lambda$ for the same values of the
parameters as in Figs. 2-7. For $\tau=1$ and $\tau=3$, it is found that the
two determinants are positive for $\lambda\in(-0.963,0.649)$ and
$\lambda\in(-1.722,1.200)$, respectively, and these values are in excellent
agreement with the thresholds $\lambda_{w}^{-}(\infty)$ and
$\lambda_{w}^{+}(\infty)$ that can be extrapolated from Fig. 5. Note that
$\lambda_{w1}$, the lower limit of the interval $\hat{\cal D}_{w}$, here
corresponds to $\det\hat{C}^{l,+}_{w}(\lambda_{w1})=0$ (with
$\det(\hat{C}^{r,+}_{w}(\lambda_{w1})+C)\neq 0$) whereas $\lambda_{w2}$, the
upper limit, corresponds to $\det(\hat{C}^{r,+}_{w}(\lambda_{w2})+C)=0$ (with
$\det\hat{C}^{l,+}_{w}(\lambda_{w2})\neq 0$)202020This is not a general
feature and the opposite behavior is observed for other values of the
parameters (for instance $g=1$ and $\tau=3$) or for the other observables..

At first sight, the latter observation may seem to contradict the fact that
the determinants of the matrices $C^{l}_{w}(\lambda,t)$ and
$C^{r}_{w}(\lambda,t)+C$ vanish for the same values of $\lambda$, i.e.,
$\lambda_{w}^{-}(t)$ or $\lambda_{w}^{+}(t)$ (see Figs. 2 and 3). To solve
this puzzle we need to look more closely into the long-time behavior of the
determinants at the boundaries of the interval $\hat{\cal
D}_{w}=(\lambda_{w1},\lambda_{w2})$. This is done in Fig. 9 for $\tau=1$ where
for conciseness we only consider the behavior of
$\det(C^{r}_{w}(\lambda,t)+C)$ in the vicinity of $\lambda_{w1}\simeq-0.963$.

Figure 9:  (Color on line) $\det(C^{r}_{w}(\lambda,t)+C)/\det C$ (symbols) and
$\det(C^{l}_{w}(\lambda,t)/\det C$ (solid lines) as a function of $\lambda$
for $\tau=1$ and $t=10,15,20$. The solid and dashed black lines represent
$\det\hat{C}_{w}^{l,+}(\lambda)/\det C$ and
$\det(\hat{C}_{w}^{r,+}(\lambda)+C)/\det C$, respectively, like in Fig. 8.
Figure 10:  (Color on line) Time evolution of $\det C_{w}^{l}(\lambda,t)/\det
C$ (solid lines) and $\det(C^{r}_{w}(\lambda,t)+C)/\det C$ (dashed lines) for
$\tau=1$ and different values of $\lambda$ close to
$\lambda_{w1}\simeq-0.963$. The red lines correspond to
$\lambda=\lambda_{w1}$.

What this figure reveals is that $\det(C^{r}_{w}(\lambda,t)+C)$ strongly
deviates from $\det(\hat{C}^{r,+}_{w}(\lambda)+C)$ as $\lambda$ decreases
before vanishing at $\lambda=\lambda_{w}^{-}(t)$. This contrasts with the
behavior of $\det C^{l}_{w}(\lambda,t)$ that smoothly approaches
$\det\hat{C}^{l,+}_{w}(\lambda)$ as $t$ increases. In addition,
$\det(C^{r}_{w}(\lambda,t)+C)$ diverges to $-\infty$ for $\lambda$ very close
to (but smaller than) $\lambda_{w}^{-}(t)$, signaling that the solution of the
RDE (42a) blows up. Another striking feature is that the value of the
determinant for $\lambda=\lambda_{w1}$ (indicated by the vertical line in the
figure) is approximately independent of time (as the three representative
curves cross each other at almost the same point). This suggests that
$\lim_{t\to\infty}\det(C^{r}_{w}(\lambda_{w1},t)+C)$ is finite but differs
from $\det(\hat{C}^{r,+}_{w}(\lambda_{w1})+C)$. This is indeed confirmed by
Fig. 10 which shows the time evolution of $\det C_{w}^{l}(\lambda,t)$ and
$\det(C^{r}_{w}(\lambda,t)+C)$ for values of $\lambda$ close to
$\lambda_{w1}$. One can see that $\det(C^{r}_{w}(\lambda_{w1},t)+C)/\det
C\approx 1.17$ at large times whereas
$\det(\hat{C}^{r,+}_{w}(\lambda_{w1})+C)/\det C\approx 2.98$.

This phenomenon is rather subtle but should come as no surprise. Indeed, as
predicted by Eq. (85a), the condition $\det C_{\rm o}^{l,+}(\lambda_{w1})=0$
implies that the solution of the RDE (42a) does not go asymptotically to
$\hat{C}_{\rm o}^{r,+}(\lambda_{w1})$ and, as a result,
$\det(C^{r}_{w}(\lambda_{w1},t\to\infty)+C)\neq\det(C^{r,+}_{w}(\lambda_{w1})+C)$.
As a matter of fact, a similar behavior takes place in the vicinity of
$\lambda_{w2}$, with the roles of $C^{r}_{w}(\lambda,t)+C$ and
$C^{l}_{w}(\lambda,t)$ inverted212121This explains the abrupt decrease of
$\det C^{l}(\lambda,t)$ in the vicinity of $\lambda_{w}^{+}(t)\approx 0.65$
for $t=10$ in Fig. 2(a)., and it is also observed for $\tau=3$. More
generally, this is one of the two scenarios that may occur at the limits of
the domain of definition of the SCGF when $\hat{\cal D}_{\rm o}\subset{\cal
D}_{H}$ (i.e., $\lambda_{\rm o1}>\lambda_{\min}$ and $\lambda_{\rm
o2}<\lambda_{\max}$). This is a nontrivial issue which will be further
discussed in Sec. IV.3. An important case is $\lambda_{\rm o2}=1$, as it
concerns the fluctuation theorem (16) for the heat and the conjecture (17) for
the apparent entropy production.

####  IV.2.2 Expressions of the SCGF and of the pre-exponential prefactors

Now that the domain of existence of the SCGF has been identified as $\hat{\cal
D}_{\rm o}=(\lambda_{\rm o1},\lambda_{\rm o,2})$, we derive explicit
expressions for the SCGF $\mu_{\rm o}(\lambda)$ and the pre-exponential
factors $g_{\rm o}(\lambda)$ when $\lambda_{\rm o1}<\lambda<\lambda_{\rm o2}$.
Since $C^{r}_{\rm o}(\lambda,t)+C\to\hat{C}^{r,+}_{\rm o}(\lambda)+C>0$ and
$C^{l}_{\rm o}(\lambda,t)\to\hat{C}^{l,+}_{\rm o}(\lambda)>0$, we have

$\displaystyle f^{r}_{\rm o}(\lambda,t)\to f^{+}(\lambda)$ $\displaystyle
f^{l}_{\rm o}(\lambda,t)\to f^{+}(\lambda)\ ,$ (107)

and from Eqs. (40) we obtain the two asymptotic expressions

$\displaystyle G^{r}_{\rm o,\lambda}({\bf u}_{0},t)\sim g^{r}_{\rm
o}(\lambda)\exp\big{(}-\frac{1}{2}{\bf u}_{0}^{T}\hat{C}^{r,+}_{\rm
o}(\lambda){\bf u}_{0}\big{)}e^{f^{+}(\lambda)t}$ (108a) $\displaystyle
G^{l}_{\rm o,\lambda}({\bf u},t)\sim g^{l}_{\rm
o}(\lambda)\exp\big{(}-\frac{1}{2}{\bf u}^{T}\hat{C}^{l,+}_{\rm
o}(\lambda){\bf u}\big{)}e^{f^{+}(\lambda)t}\ ,$ (108b)

with

$\displaystyle g^{r}_{\rm o}(\lambda)$
$\displaystyle=\exp\left(\int_{0}^{\infty}dt\>[f^{r}_{\rm
o}(\lambda,t)-f^{+}(\lambda)]\right)$ (109a) $\displaystyle g^{l}_{\rm
o}(\lambda)$ $\displaystyle=\left(\frac{\det
C}{(2\pi)^{n+2}}\right)^{1/2}\exp\left(\int_{0}^{\infty}dt\>[f^{l}_{\rm
o}(\lambda,t)-f^{+}(\lambda)]\right)\ .$ (109b)

Integrating $G^{r}_{\rm o,\lambda}({\bf u}_{0},t)p({\bf u}_{0})$ over the
initial state ${\bf u}_{0}$ and $G^{l}_{\rm o,\lambda}({\bf u},t)$ over the
final state ${\bf u}$, we finally get

$\displaystyle G_{\rm o,\lambda}(t)\sim g_{\rm o}(\lambda)e^{f^{+}(\lambda)t}\
,$ (110)

with

$\displaystyle g_{\rm o}(\lambda)$
$\displaystyle=\Big{[}\frac{\det[C+\hat{C}^{r,+}_{\rm o}(\lambda)]}{\det
C}\Big{]}^{-1/2}\exp\left(\int_{0}^{\cal\infty}dt\>[f^{r}_{\rm
o}(\lambda,t)-f^{+}(\lambda)]\right)$ (111a)
$\displaystyle=\Big{[}\frac{\det\hat{C}^{l,+}_{\rm o}(\lambda)}{\det
C}\Big{]}^{-1/2}\exp\left(\int_{0}^{\cal\infty}dt\>[f^{l}_{\rm
o}(\lambda,t)-f^{+}(\lambda)]\right)\ .$ (111b)

Eq. (110) is a central result of our work as it shows that for
$\lambda\in\hat{\cal D}_{\rm o}$ the SCGF $\mu_{\rm o}(\lambda)$ defined in
Eq. (29) is equal to the function $f^{+}(\lambda)$ which is associated with
the maximal solutions of the CAREs (86). The SCGF is therefore the same for
the three observables ${\cal W}_{t},{\cal Q}_{t}$ and $\Sigma_{t}$ since the
corresponding matrices $B_{\rm o}$ have the same anti-symmetric part, which
implies that the Hamiltonian matrices have the same spectrum. On the other
hand, the interval $\hat{\cal D}_{\rm o}$ and the pre-exponential factor
$g_{\rm o}(\lambda)$ depend on the observable. In particular, the expressions
(111) of $g_{\rm o}(\lambda)$ diverge at the boundaries of $\hat{\cal D}_{\rm
o}$222222However, this does not mean that the actual prefactor diverges at
$\lambda_{\rm o1}$ or $\lambda_{\rm o2}$, as discussed in Sec. IV.3..

The SCGF can be computed from Eq. (99), i.e., from the eigenvalues of the
Hamiltonian matrices. This is a standard numerical task, even for large $n$,
but we now derive an equivalent and even more convenient expression as an
integral over frequency. To this end, we use the property that each symmetric
solution, $\hat{C}_{\rm o}^{r,(\alpha)}(\lambda)$ or $\hat{C}_{\rm
o}^{l,(\alpha)}(\lambda)$, of the CAREs gives rise to a factorization of the
characteristic polynomial $p_{H}(\lambda,s)$ of $H^{r}_{\rm o}(\lambda)$ and
$H^{l}_{\rm o}(\lambda)$ in the form K2010

$\displaystyle
p_{H}(\lambda,s)=(-1)^{n}q^{(\alpha)}(\lambda,s)q^{(\alpha)}(\lambda,-s)\ ,$
(112)

where $q^{(\alpha)}(\lambda,s)$ is the characteristic polynomial of the
matrices $-A_{\rm o}(\lambda)+D\hat{C}_{\rm o}^{r,(\alpha)}(\lambda)$ and
$A_{\rm o}(\lambda)+D\hat{C}_{\rm o}^{l,(\alpha)}(\lambda)$. (Owing to the
invariance relations (91), $q^{(\alpha)}(\lambda,s)$ does not depend on the
observable.) Eq. (112) follows from

$\displaystyle\begin{bmatrix}I_{n+2}&0\\\ \hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)&I_{n+2}&\end{bmatrix}^{-1}H^{r}_{\rm
o}(\lambda)\begin{bmatrix}I_{n+2}&0\\\ \hat{C}^{r,(\alpha)}_{\rm
o}(\lambda)&I_{n+2}&\end{bmatrix}=\begin{bmatrix}-A_{\rm
o}(\lambda)+D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)&D\\\ 0&[A_{\rm
o}(\lambda)-D\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)]^{T}&\end{bmatrix}$ (113)

and a similar relation involving $C^{l,(\alpha)}_{\rm o}(\lambda)$ and
$H^{l}_{\rm o}(\lambda)$. Focusing on the maximal solutions
$\hat{C}^{r,+}_{\rm o}(\lambda)$ and $\hat{C}^{l,+}_{\rm o}(\lambda)$, we then
define the rational function

$\displaystyle G^{+}(\lambda,s)=\frac{q^{+}(\lambda,s)}{p_{A}(s)}\ ,$ (114)

where $p_{A}(s)$ is the characteristic polynomial of the drift matrix $A$
given by Eq. (24). As a result, $G^{+}(\lambda,s)$ has $n+2$ zeros but no
poles in the closed right half plane (since all eigenvalues of $A$ have a
negative real part). In addition,

$\displaystyle
q^{+}(\lambda,s)=\prod_{i=1}^{n+2}(s-s_{i}^{+})=s^{n+2}-(\sum_{i}^{n+2}s_{i}^{+})s^{n+1}+...$
(115)

and

$\displaystyle p_{A}(s)=s^{n+2}-\mbox{Tr}(A)s^{n+1}+...\ ,$ (116)

so that the function $L^{+}(\lambda,s)$ defined through

$\displaystyle G^{+}(\lambda,s)=[1+L^{+}(\lambda,s)]^{-1}$ (117)

has a relative degree equal to $1$. (The relative degree is the difference
between the degrees of the polynomials in the denominator and in the
numerator. In control theory, $L^{+}_{\lambda}(s)$ would be interpreted as a
proper, scalar rational loop transfer function of a feedback system B2021 ;
AM2008 .) A standard application of Cauchy’s residue theorem, known as Bode’s
sensitivity integral in control theory B2021 ; AM2008 , then yields

$\displaystyle\frac{1}{2\pi
i}\int_{-i\infty}^{i\infty}ds\>\ln|G^{+}(\lambda,s)|=\sum_{i=1}^{n+2}s_{i}^{+}(\lambda)-\frac{1}{2}\kappa^{+}(\lambda)\
,$ (118)

where the integration is performed along the imaginary axis in the complex
$s$-plane, and

$\displaystyle\kappa^{+}(\lambda)$
$\displaystyle\equiv\lim_{s\to\infty}sL^{+}(\lambda,s)=\lim_{s\to\infty}s\frac{p_{A}(s)-q^{+}(\lambda,s)}{p_{A}(s)}$
$\displaystyle=\sum_{i=1}^{n+2}s_{i}^{+}(\lambda)-\mbox{Tr}(A)\ .$ (119)

Inserting Eqs. (118) and (IV.2.2) into Eq. (99), we obtain that

$\displaystyle f^{+}(\lambda)=-\frac{1}{2\pi
i}\int_{-i\infty}^{i\infty}ds\>\ln|G^{+}(\lambda,s)|,$

which thanks to Eqs. (112), (114), and (72) can then be rewritten as

$\displaystyle f^{+}(\lambda)$ $\displaystyle=-\frac{1}{4\pi
i}\int_{-i\infty}^{i\infty}ds\>\ln\frac{(-1)^{n}p_{H}(\lambda,s)}{|p_{A}(s)|^{2}}$
$\displaystyle=-\frac{1}{4\pi
i}\int_{-i\infty}^{i\infty}ds\>\ln\Big{[}1-\frac{2\lambda
g}{Q_{0}^{2}}(\frac{n}{\tau})^{2n}s\frac{(1-\frac{s\tau}{n})^{n}-(1+\frac{s\tau}{n})^{n}}{|p_{A}(s)|^{2}}\Big{]}\
.$ (121)

Note that the integral is finite as $p_{H}(\lambda,s)$ has no purely imaginary
roots for $\lambda\in{\cal D}_{H}$.

We can further transform Eq. (IV.2.2) by replacing the Laplace variable $s$ by
the frequency $\omega=is$ and introducing the response function
$\chi_{n}(\omega)$ defined in Fourier space by
$x(\omega)=\chi_{n}(\omega)\xi(\omega)$. It is readily found from Eq. (4) that

$\displaystyle\chi_{n}(\omega)=\Big{(}-\omega^{2}-\frac{i\omega}{Q_{0}}+1-\frac{g}{Q_{0}}(1-\frac{i\omega\tau}{n})^{-n}\Big{)}^{-1}\
,$ (122)

which can be also expressed as

$\displaystyle\chi_{n}(\omega)$
$\displaystyle=(\frac{n}{\tau})^{n}\frac{(1-\frac{i\omega\tau}{n})^{n}}{p_{A}(s=-i\omega)}$
(123)

by using Eq. (24). Eq. (IV.2.2) is then recast as

$\displaystyle f^{+}(\lambda)$
$\displaystyle=-\frac{1}{2}\int\frac{d\omega}{2\pi}\>\ln\Big{[}1-\frac{4\lambda
g}{Q_{0}^{2}}\omega\sin_{n}(\omega\tau)|\chi_{n}(\omega)|^{2}\Big{]}\ ,$ (124)

where we have defined the $n$-dependent sine-like function

$\displaystyle\sin_{n}(x)=\frac{1}{2i}[(1-\frac{ix}{n})^{-n}-(1+\frac{ix}{n})^{-n})]$
(125)

which satisfies $\lim_{n\to\infty}\sin_{n}(x)=\sin(x)$. As a consequence,
$f^{+}(\lambda)$ is simply expressed in terms of the spectral density
$|\chi_{n}(\omega)|^{2}$. This shows the connection between the Riccati
formalism and the results in the mathematical literature BD1997 ; BGR1997 ;
GRZ1999 ; ZS2023  for the SCGF of stationary Gaussian processes232323These
results concern quadratic observables instead of linear currents but one can
easily adapt the present derivation to such observables and obtain an integral
representation of the SCGF similar to one given in the first line of Eq.
(IV.2.2)..

Figure 11:  (Color on line) (a) Function $f^{+}(\lambda)$ computed from Eq.
(124) for $n=5$, $\tau=1$ and $\tau=3$. The parts of the curves in red
correspond to the interval $\hat{\cal D}_{w}$ for which the SCGF
$\mu_{w}(\lambda)$ is finite and equal to $f^{+}(\lambda)$ (see Fig. 8) and
the vertical dashed lines indicate the limits $\lambda_{\min}$ and
$\lambda_{\max}$ of ${\cal D}_{H}$, the interval of definition of
$f^{+}(\lambda)$ ($\lambda_{\min}\approx-472.8$ for $\tau=1$, which is outside
the scale of the figure). (b) Corresponding rate function $I_{w}(a)$. The
parts of the curves in red are obtained from the Legendre transform (126) and
the linear parts in black are obtained from Eqs. (127).

The behavior of $f^{+}(\lambda)$ with $\lambda$ for $n=5$, $\tau=1$ and
$\tau=3$ is shown in Fig. 11(a). This illustrates two important properties of
the function: i) $f^{+}(\lambda)$ has infinite slopes at the boundaries of
${\cal D}_{H}$ and ii) $f^{+}(\lambda)$ is convex in ${\cal D}_{H}$ (we
already know that this function is differentiable). Accordingly, the Legendre-
Fenchel transform (31) that relates the SCGF to the rate function reduces to
the usual Legendre transform. Specifically, defining $\lambda^{*}(a)$ as the
(unique) solution of $df^{+}(\lambda)/d\lambda|_{\lambda=\lambda^{*}}=-a$, one
has T2009

$\displaystyle
I(a)=-a\lambda^{*}-f^{+}(\lambda^{*})\>\>\>\mbox{for}\>\>\lambda^{*}(a)\in\hat{\cal
D}_{\rm o}\ .$ (126)

If $\hat{\cal D}_{\rm o}={\cal D}_{H}$ the function $I(a)$, which does not
depend on the observable, is asymptotically linear as $a\to-\infty$ (resp.
$a\to+\infty$) with slope $-\lambda_{\max}$ (resp. $-\lambda_{\min}$)242424Due
to the minus sign in our definition of $G_{\rm o,\lambda}(t)$ (which is the
same as in Ref. J2020 ), the behavior of the SCGF for $\lambda>0$ (resp.
$\lambda<0$) is relevant for the rate function for $a<\langle a\rangle$ (resp.
$a>\langle a\rangle$).. This means that the values of $\lambda$ outside ${\cal
D}_{H}$ are irrelevant to the determination of the rate function, which
justifies our initial choice to restrict the study of the solutions of the
Riccati equations to this interval.

On the other hand, if $\hat{\cal D}_{\rm o}\subset{\cal D}_{H}$ (where here
and below $\subset$ means a strict inclusion), the condition
$\lambda^{*}(a)\in\hat{\cal D}_{\rm o}$ defines two special values $a_{\rm
o}^{-}=-df^{+}(\lambda)/d\lambda|_{\lambda_{\rm o2}}$ and $a_{\rm
o}^{+}=-df^{+}(\lambda)/d\lambda|_{\lambda_{\rm o1}}$ beyond which the rate
function is no longer given by Eq. (126). Since $f^{+}(\lambda)$ has finite
slopes at the two cutoff $\lambda_{\rm o1}$ and $\lambda_{\rm o2}$, the
properties of Legendre transforms T2009  imply that the rate function (which
now depends on the observable) exhibits linear branches beyond $a_{\rm
o}^{\pm}$,

$\displaystyle I_{\rm o}(a)$ $\displaystyle=-a\lambda_{\rm
o2}-f^{+}(\lambda_{\rm o2})\>\>\>\>\>a\leq a_{\rm o}^{-}$ (127a)
$\displaystyle I_{\rm o}(a)$ $\displaystyle=-a\lambda_{\rm
o1}-f^{+}(\lambda_{\rm o1})\>\>\>\>\>a\geq a_{\rm o}^{+}\ .$ (127b)

This yields the curves $I_{w}(a)$ plotted in Fig. 11(b).

Finally, we consider the limit $n\to\infty$. It can be directly taken in Eq.
(IV.2.2), replacing $(1\pm s\tau/n)^{n}$ by $e^{\pm s\tau}$ and using Eq.
(25). This yields

$\displaystyle\lim_{n\to\infty}\frac{(-1)^{n}p_{H}(\lambda,s)}{p_{A}(s)p_{A}(-s)}=1-\frac{2\lambda
g}{Q_{0}^{2}}\frac{s(e^{-s\tau}-e^{s\tau})}{(s^{2}+\frac{s}{Q_{0}}+1-\frac{g}{Q_{0}}e^{-s\tau})(s^{2}-\frac{s}{Q_{0}}+1-\frac{g}{Q_{0}}e^{s\tau})}\
.$ (128)

Alternatively, from Eq. (124),

$\displaystyle f^{+}_{\infty}(\lambda)\equiv\lim_{n\to\infty}f^{+}(\lambda)$
$\displaystyle=-\frac{1}{2}\int\frac{d\omega}{2\pi}\>\ln\Big{[}1-\frac{4\lambda
g}{Q_{0}^{2}}\omega\sin(\omega\tau)|\chi_{\infty}(\omega)|^{2}\Big{]}\ ,$
(129)

with

$\displaystyle\chi_{\infty}(\omega)=(-\omega^{2}-\frac{i\omega}{Q_{0}}+1-\frac{g}{Q_{0}}e^{i\omega\tau})^{-1}\
.$ (130)

As expected, this result coincides with the expression of the SCGF for the
discrete delay obtained in our previous work RTM2017  by imposing periodic
boundary conditions on the solution of the Langevin equation (i.e.,
$x(0)=x(t))$ and expanding $x(t)$ in a Fourier series (this amounts to
assuming that boundary conditions can be neglected to leading order in $t$ in
the calculation of $\langle e^{-\lambda{\cal A}_{t}}\rangle$; see, e.g., Ref.
ZBCK2005  for a similar calculation). It is clear that Eq. (124) can be also
derived by using the same method for $n$ finite.

####  IV.2.3 Relation with the spectral problem for the tilted generators

As is well known, the SCGF of additive functionals such as the ones considered
in this work is given (in its domain of definition) by the dominant eigenvalue
$\mu_{\rm o}(\lambda)$ of the tilted generators ${\cal L}_{\rm o,\lambda}$ and
${\cal L}^{\dagger}_{\rm o,\lambda}$ T2009 ; J2020 ; CT2015 ; T2018 . This
results from the expansion of the restricted generating function $G_{\rm
o,\lambda}({\bf u},t|{\bf u}_{0})$ in a complete basis of bi-orthogonal
eigenfunctions, which yields asymptotically CT2015

$\displaystyle G_{\rm o,\lambda}({\bf u},t|{\bf u}_{0})\sim r_{\rm
o,\lambda}({\bf u}_{0})l_{\rm o,\lambda}({\bf u})e^{\mu_{\rm o}(\lambda)t}\ ,$
(131)

where $r_{\rm o,\lambda}({\bf u}_{0})$ and $l_{\rm o,\lambda}({\bf u})$ are
the right and left eigenfunctions associated with $\mu_{\rm o}(\lambda)$.
Finding the SCGF thus requires solving the spectral problem

$\displaystyle{\cal L}_{\rm o,\lambda}\,r_{\rm o,\lambda}({\bf u}_{0})$
$\displaystyle=\mu_{\rm o}(\lambda)\,r_{\rm o,\lambda}({\bf u}_{0})$ (132a)
$\displaystyle{\cal L}_{\rm o,\lambda}^{{\dagger}}\,l_{\rm o,\lambda}({\bf
u})$ $\displaystyle=\mu_{\rm o}(\lambda)\,l_{\rm o,\lambda}({\bf u})\ ,$
(132b)

with the eigenfunctions commonly normalized according to CT2015

$\displaystyle\int d{\bf u}\>l_{\rm o,\lambda}({\bf u})r_{\rm o,\lambda}({\bf
u})=1$ (133a) $\displaystyle\int d{\bf u}\>l_{\rm o,\lambda}({\bf u})=1\ .$
(133b)

(In particular, Eq. (133a) is inherited from the duality between the operators
${\cal L}_{\rm o,\lambda}$ and ${\cal L}^{\dagger}_{\rm o,\lambda}$ which
imposes the boundary condition $l_{\rm o,\lambda}({\bf u})r_{\rm
o,\lambda}({\bf u})=0$ at infinity T2018 .) Since we have shown above that the
SCGF $\mu_{\rm o}(\lambda)$ is equal to $f^{+}(\lambda)$ inside its domain of
definition $\hat{\cal D}_{\rm o}$, we already know the solution of the
spectral problem for the dominant eigenvalue,252525It may be confusing that
$f^{+}(\lambda)$ is the largest eigenvalue of the operators ${\cal L}_{\rm
o,\lambda}$ and ${\cal L}^{\dagger}_{\rm o,\lambda}$ but is smaller than of
all other functions $f^{(\alpha)}(\lambda)$ according to inequality (101).
There is no contradiction though because the corresponding matrices
$\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)$ and $\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)$ are not valid solutions of the spectral problem (132). Indeed,
the condition $\hat{C}^{r,(\alpha)}_{\rm o}(\lambda)+\hat{C}^{l,(\alpha)}_{\rm
o}(\lambda)>0$ is only satisfied by the maximal solutions of the CAREs [Eq.
98)]. Take for instance the matrices $\hat{C}^{r,-}_{\rm o}(\lambda)$ and
$\hat{C}^{l,-}_{\rm o}(\lambda)$ that are the minimal solutions of the CAREs.
From Eq. (96), $f^{-}(\lambda)$ is the largest of all functions
$f^{(\alpha)}(\lambda)$ since all eigenvalues of the Hamiltonian matrices in
the set ${\cal S}_{\lambda}^{-}$ have a negative real part. On the other hand,
relations (III.3) imply that $\hat{C}^{r,-}_{\rm
o}(\lambda)+\hat{C}^{l,-}_{\rm o}(\lambda)=-[\hat{C}^{r,+}_{\rm
o}(\lambda)+\hat{C}^{l,+}_{\rm o}(\lambda)]$ and therefore all the eigenvalues
of this matrix are negative. As a result, the condition (133a) is only
satisfied by the eigenfunctions given by Eqs. (134). and the eigenfunctions
corresponding to $f^{+}(\lambda)$ are readily obtained from the asymptotic
expressions (108) as

$\displaystyle r_{\rm o,\lambda}({\bf u}_{0})$ $\displaystyle=g_{\rm
o}^{r}(\lambda)\exp\left(-\frac{1}{2}{\bf u}_{0}^{T}\hat{C}^{r,+}_{\rm
o}(\lambda){\bf u}_{0}\right)$ (134a) $\displaystyle l_{\rm o,\lambda}({\bf
u})$ $\displaystyle=\frac{g_{\rm o}^{l}(\lambda)}{\int d{\bf u}_{0}\>p({\bf
u}_{0})r_{\rm o,\lambda}({\bf u}_{0})}\exp\left(-\frac{1}{2}{\bf
u}^{T}\hat{C}^{l,+}_{\rm o}(\lambda){\bf u}\right)\ ,$ (134b)

with $g_{\rm o}^{r}(\lambda)$ and $g_{\rm o}^{l}(\lambda)$ given by Eqs.
(109). These expressions can be used to derive a simple expression of the
prefactor $g_{\rm o}(\lambda)$. Inserting Eqs. (134) into Eq. (133a) and
performing a few manipulations, we obtain the two relations

$\displaystyle\exp\left(\int_{0}^{\cal\infty}dt\>[f^{r}_{\rm
o}(\lambda,t)-f^{+}(\lambda)]\right)=\Big{[}\frac{\det[\hat{C}^{r,+}_{\rm
o}(\lambda)+\hat{C}^{l,+}_{\rm o}(\lambda)]}{\det\hat{C}^{l,+}_{\rm
o}(\lambda)}\Big{]}^{1/2}$ (135a)
$\displaystyle\exp\left(\int_{0}^{\cal\infty}dt\>[f^{l}_{\rm
o}(\lambda,t)-f^{+}(\lambda)]\right)=\Big{[}\frac{\det[\hat{C}^{r,+}_{\rm
o}(\lambda)+\hat{C}^{l,+}_{\rm o}(\lambda)]}{\det[C+\hat{C}^{r,+}_{\rm
o}(\lambda)]}\Big{]}^{1/2}\ ,$ (135b)

and using Eqs. (111) we find

$\displaystyle g_{\rm o}(\lambda)=\left(\frac{\det C\>\det[\hat{C}^{l,+}_{\rm
o}(\lambda)+\hat{C}^{r,+}_{\rm o}(\lambda)]}{\det\hat{C}^{l,+}_{\rm
o}(\lambda)\>\det[\hat{C}^{r,+}_{\rm o}(\lambda)+C]}\right)^{1/2}\ .$ (136)

This latter expression is another significant result of this work. It has the
great advantage of no longer involving an integration over time; it suffices
to compute the maximal solutions of the CAREs (86) from Eq. (88), which is a
simple task262626As those in Eqs. (111), this expression diverges when
$C^{l,+}_{\rm o}(\lambda)$ and/or $\hat{C}^{r,+}_{\rm o}(\lambda)+C$ are no
longer positive definite, i.e., at the limits of the interval $\hat{\cal
D}_{\rm o}$. However, we stress again that this does not imply that the actual
prefactor is infinite at $\lambda_{\rm o1}$ or $\lambda_{\rm o2}$, as will be
shown in Sec. IV.3..

Finally, after using Eqs. (135) and the normalization integral (133b), the
expressions of $r_{\rm o,\lambda}({\bf u}_{0})$ and $l_{\rm o,\lambda}({\bf
u})$ can be rewritten as

$\displaystyle r_{\rm o,\lambda}({\bf u}_{0})$
$\displaystyle=\Big{[}\frac{\det[\hat{C}^{r,+}_{\rm
o}(\lambda)+\hat{C}^{l,+}_{\rm o}(\lambda)]}{\det\hat{C}^{l,+}_{\rm
o}(\lambda)}\Big{]}^{1/2}\exp\left(-\frac{1}{2}{\bf
u}_{0}^{T}\hat{C}^{r,+}_{\rm o}(\lambda){\bf u}_{0}\right)$ (137a)
$\displaystyle l_{\rm o,\lambda}({\bf u})$
$\displaystyle=\Big{[}\frac{\det\hat{C}^{l,+}_{\rm
o}(\lambda)}{(2\pi)^{n+2}}\Big{]}^{1/2}\exp\left(-\frac{1}{2}{\bf
u}^{T}\hat{C}^{l,+}_{\rm o}(\lambda){\bf u}\right)\ .$ (137b)

One can notice the similarity of these equations with Eqs. (55) and (56) in
Ref. KSD2011 . In this paper, the large-time expression of the generating
function for the heat flow in harmonic chains was computed by using finite-
time Fourier transforms. As already mentioned, the expression (124) of the
SCGF is easily obtained by this method. With more efforts, one can also
compute the eigenfunctions, i.e., the matrices $\hat{C}^{r,+}_{\rm
o}(\lambda)$ and $\hat{C}^{l,+}_{\rm o}(\lambda)$, and in turn the pre-
exponential factors $g_{\rm o}(\lambda)$ (Eq. (136) is similar to Eq. (59) in
Ref. KSD2011 ). However, this method has a major drawback which makes it
impracticable for multidimensional systems: Each component of the matrices is
given by an integral over frequency, so that the numerical computation becomes
more and more burdensome as the dimensionality increases272727As a matter of
fact, the method has only been applied to one-dimensional systems Sab2011 ..
Furthermore, one does not have access to the finite-time behavior of the
generating functions and one cannot study the special cases where the
solutions of the RDEs do not converge to $\hat{C}^{r,+}_{\rm o}(\lambda)$ and
$\hat{C}^{l,+}_{\rm o}(\lambda)$ (see Sec. V). For all these reasons, the
Riccati approach is much more illuminating and numerically effective.

####  IV.2.4 Effective process

From the knowledge of the eigenfunction $r_{\rm o,\lambda}({\bf u})$, one can
build the so-called effective or driven process that describes how
fluctuations of the observables are created dynamically in the long-time limit
CT2015 . By construction, a given fluctuation $a$ of the time-intensive
observable ${\cal A}_{t}/t$ in the original process is realized as a typical
value in the effective process.

Since $r_{\rm o,\lambda}({\bf u})$ is a multivariate Gaussian, this process is
again a linear diffusion,

$\displaystyle\dot{{\bf u}}_{t}=\hat{A}(\lambda){\bf
u}_{t}+{\boldsymbol{\xi}}_{t}\,,$ (138)

with the same diffusion matrix $D$ as the original process but with a modified
drift $\hat{\bf F}_{\lambda}=\hat{A}(\lambda){\bf u}$ given by CT2015

$\displaystyle\hat{\bf F}_{\lambda}$ $\displaystyle={\bf F}+D(\nabla\ln r_{\rm
o,\lambda}-\lambda{\bf g}_{\rm o})$ $\displaystyle=[A-D(\hat{C}^{r,+}_{\rm
o}(\lambda)+\lambda B_{\rm o})]{\bf u}\ .$ (139)

Hence, from the definition of the matrix $A_{\rm o}(\lambda)$ [Eq. (45)],

$\displaystyle\hat{A}(\lambda)=A_{\rm o}(\lambda)-D\hat{C}^{r,+}_{\rm
o}(\lambda).$ (140)

By definition of the maximal solution $\hat{C}^{r,+}_{\rm o}(\lambda)$, the
eigenvalues of $-\hat{A}(\lambda)$ belong to the subset ${\cal
S}_{\lambda}^{+}$ of eigenvalues of the Hamiltonian matrices [see Eq. (87)]
and therefore $-\hat{A}(\lambda)>0$. In addition, owing to Eqs. (91a) and
(92a), $\hat{A}(\lambda)$ only depends on the antisymmetric part of the matrix
$B_{\rm o}$. As a result, the driven process is the same for all observables
having the same antisymmetric part (but it is only defined for
$\lambda\in\hat{\cal D}_{\rm o}$, the domain of definition of the SCGF, which
depends on the observable).

The corresponding invariant density is then given by CT2015

$\displaystyle\hat{p}_{\lambda}({\bf u})$ $\displaystyle=r_{\rm
o,\lambda}({\bf u})l_{\rm o,\lambda}({\bf
u})=\big{[}(2\pi)^{n+2}\det\hat{\Sigma}(\lambda)\big{]}^{-1/2}\exp\left(-\frac{1}{2}{\bf
u}^{T}\hat{\Sigma}^{-1}(\lambda){\bf u}\right)\ ,$ (141)

with

$\displaystyle\hat{\Sigma}(\lambda)=[\hat{C}^{r,+}_{\rm
o}(\lambda)+\hat{C}^{l,+}_{\rm o}(\lambda)]^{-1}\ .$ (142)

Using the fact that $\hat{C}^{r,+}_{\rm o}(\lambda)$ and $\hat{C}^{l,+}_{\rm
o}(\lambda)$ are solutions of the CAREs (86a) and (86b), respectively, one can
easily verify that $\hat{A}(\lambda)$ and $\hat{\Sigma}(\lambda)$ are related
via the Lyapunov equation

$\displaystyle\hat{A}(\lambda)\hat{\Sigma}(\lambda)+\hat{\Sigma}(\lambda)\hat{A}^{T}(\lambda)=-D\
,$ (143)

as it must be. We recall that the sum $\hat{C}^{r,+}_{\rm
o}(\lambda)+\hat{C}^{l,+}_{\rm o}(\lambda)$ does not depend on the observable
due to Eq. (90) and is positive definite, which guarantees that the invariant
density exists and the driven process is ergodic.

The above equations are quite general and apply to any multidimensional linear
diffusions282828Eqs. (140) and (141) correspond to Eqs. (65) and (54) in Ref.
DBT2023  via the changes $\hat{A}(\lambda)\to-M_{k}$ and
$\hat{\Sigma}(\lambda)\to C_{k}$. but they take a simpler form for the model
under study owing to the fact that $D_{ij}=\delta_{i1}\delta_{j1}$. The
elements of the matrix $\hat{A}(\lambda)$ are then given by

$\displaystyle\hat{A}_{ij}(\lambda)=[A_{\rm
o}(\lambda)]_{ij}-[\hat{C}^{r,+}_{\rm o}(\lambda)]_{1j}\delta_{i1}.$ (144)

Moreover, since $\hat{A}(\lambda)$ does not depend on the observable, we can
choose $\rm o=w$, which yields

$\displaystyle\hat{A}_{ij}(\lambda)=A_{ij}-[\hat{C}^{r,+}_{w}(\lambda)]_{1j}\delta_{i1}\
.$ (145)

Therefore, only the first equation of the set of equations (II.1) is modified
and replaced by

$\displaystyle\dot{v}(t)$
$\displaystyle=-\big{(}\frac{1}{Q_{0}}+[\hat{C}_{w}^{r,+}(\lambda)]_{11}\big{)}v(t)-x(t)+\frac{g}{Q_{0}}x_{n}(t)-\sum_{j=0}^{n}\>[\hat{C}_{w}^{r,+}(\lambda)]_{1(j+2)}\>x_{j}(t)+\xi(t)$
$\displaystyle=-\big{[}\frac{1}{Q_{0}}-2f^{+}(\lambda)\big{]}v(t)-\big{[}1+2f^{+}(\lambda)(f^{+}(\lambda)-\frac{1}{Q_{0}})\big{]}x(t)+\frac{g}{Q_{0}}x_{n}(t)-\sum_{j=1}^{n}\>[\hat{C}_{w}^{r,+}(\lambda)]_{1(j+2)}\>x_{j}(t)+\xi(t)\
,$ (146)

where we have used that $[\hat{C}_{w}^{r,+}(\lambda)]_{11}=-2f^{+}(\lambda)$
[Eq. (94a)] and
$[\hat{C}_{w}^{r,+}(\lambda)]_{12}=2f^{+}(\lambda)(f^{+}(\lambda)-1/Q_{0})$ to
derive the second line292929More generally, the element $(12)$ of a general
solution of the CARE ${\cal
R}_{w,\lambda}[\hat{C}_{w}^{r,(\alpha)}(\lambda)]=0$ is equal to
$2f^{r,(\alpha)}(\lambda)[f^{r,(\alpha)}(\lambda)-1/Q_{0}]$. One also has
$[\hat{C}_{w}^{r,(\alpha)}(\lambda)]_{23}=(2\tau/n)f^{r,(\alpha)}(\lambda)\big{[}(f^{r,(\alpha)}(\lambda))^{2}-f^{r,(\alpha)}(\lambda)/Q_{0}+1\big{]}\big{[}f^{r,(\alpha)}(\lambda)-1/Q_{0}\big{]}$.
The other elements of the matrix $\hat{C}_{w}^{r,(\alpha)}(\lambda)$ are not
simply expressed in terms of the function $f^{r,(\alpha)}(\lambda)$.. After
inserting the definition of the auxiliary variables $x_{j}(t)$ [Eq. (8)], the
equation is finally recast as

$\displaystyle\dot{v}(t)=-\big{[}\frac{1}{Q_{0}}-2f^{+}(\lambda)\big{]}v(t)-\big{[}1+2f^{+}(\lambda)(f^{+}(\lambda)-\frac{1}{Q_{0}})\big{]}x(t)+\int_{-\infty}^{t}ds\>[\frac{g}{Q_{0}}\>g_{n}(t-s,n/\tau)-h_{n}(\lambda,t-s)]x(s)+\xi(t)\
,$ (147)

where

$\displaystyle
h_{n}(\lambda,t)=\sum_{j=1}^{n}\>[\hat{C}_{w}^{r,+}]_{1(j+2)}\>g_{j}(t,n/\tau)\
.$ (148)

Eq. (147) is another major result of this work. To our knowledge, this is the
first explicit example of an effective process for a non-Markovian Langevin
dynamics. Comparing with the original process governed by Eq. (4), we see that
atypical fluctuations of the observable are created in the long-time limit by
modifying not only the friction and the spring constants but also the memory
kernel.

The stationary density associated with Eq. (147) is obtained by tracing out
the auxiliary variables $u_{j}$ ($j\geq 2$) in Eq. (141), which leads to the
bivariate Gaussian distribution

$\displaystyle\hat{p}_{\lambda}(x,v)=\frac{1}{\pi
Q_{0}\sqrt{\frac{\hat{T}_{x}(\lambda)}{T}\frac{\hat{T}_{v}(\lambda)}{T}}}e^{-\frac{1}{Q_{0}}(\frac{T}{\hat{T}_{x}(\lambda)}x^{2}+\frac{T}{\hat{T}_{v}(\lambda)}v^{2})}$
(149)

characterized by the two $\lambda$-dependent temperatures
$\hat{T}_{x}(\lambda)=(2T/Q_{0})\langle
x^{2}\rangle_{\lambda}=(2T/Q_{0})\hat{\Sigma}_{22}(\lambda)$ and
$\hat{T}_{v}(\lambda)=(2T/Q_{0})\langle
v^{2}\rangle_{\lambda}=(2T/Q_{0})\hat{\Sigma}_{11}(\lambda)$ (here
$\langle...\rangle_{\lambda}$ denotes an average over stochastic trajectories
generated by the effective process in the stationary limit). Interestingly,
the variances $\langle x^{2}\rangle_{\lambda}$ and $\langle
v^{2}\rangle_{\lambda}$, and more generally all elements of the covariance
matrix $\hat{\Sigma}(\lambda)$, can be expressed in terms of the spectral
density $|\chi_{n}(\omega)|^{2}$, as the SCGF $f^{+}(\lambda)$.

To show this, we write the solution of the Lyapunov equation (143) as an
integral over frequency of the spectrum matrix. Since
$D_{ij}=\delta_{i1}\delta_{j1}$, we have G2004

$\displaystyle\hat{\Sigma}_{ij}(\lambda)=\int\frac{d\omega}{2\pi}\langle{\bf
u}(\omega){\bf
u}^{T}(-\omega)\rangle_{\lambda}=\int\frac{d\omega}{2\pi}\hat{R}_{i1}(\lambda,\omega)\hat{R}_{j1}(\lambda,-\omega)\
,$ (150)

where $\hat{R}(\lambda,t)=e^{\hat{A}(\lambda)t}$ ($t>0$) is the response or
Green’s function which reads in the Fourier space

$\displaystyle\hat{R}(\lambda,\omega)=\int_{0}^{\infty}e^{\hat{A}(\lambda)t}e^{i\omega
t}dt=-[\hat{A}(\lambda)+i\omega I_{n+2}]^{-1}\ .$ (151)

Thanks to Eq. (145), the first column of $\hat{R}(\lambda,\omega)$ does not
depend on $[\hat{C}^{r,+}_{\rm o}(\lambda)]_{1j}$ ($j=1,...,n+2)$, and from
the expression of the drift matrix $A$ [Eq. (9)] we find

$\displaystyle\hat{R}_{11}(\lambda,\omega)=\frac{i\omega(i\omega-\frac{n}{\tau})^{n}}{q^{+}(\lambda,s=i\omega)}$
(152)

and

$\displaystyle\hat{R}_{i1}(\lambda,\omega)=\frac{(-1)^{i-1}(\frac{n}{\tau})^{i-2}(i\omega-\frac{n}{\tau})^{n+2-i}}{q^{+}(\lambda,s=i\omega)}\
,\>\>\>i=2,3,...,n+2\ ,$ (153)

where $q^{+}(\lambda,s)=\prod_{i=1}^{n+2}[s-s_{i}^{+}(\lambda)]$ [Eq. (115)].
Inserting the above expressions into Eq. (150) and using Eqs. (72), (122),
(123) and (125), we then obtain

$\displaystyle\frac{\hat{T}_{x}(\lambda)}{T}=\frac{2}{Q_{0}}\hat{\Sigma}_{22}(\lambda)$
$\displaystyle=\frac{2}{Q_{0}}\int\frac{d\omega}{2\pi}\frac{(\omega^{2}+\frac{n^{2}}{\tau^{2}})^{n}}{(-1)^{n}p_{H,\lambda}(s=i\omega)}$
$\displaystyle=\frac{2}{Q_{0}}\int\frac{d\omega}{2\pi}\frac{1}{|\chi_{n}(\omega)|^{-2}-\frac{4\lambda
g}{Q_{0}^{2}}\omega\sin_{n}(\omega\tau)}$ (154a)
$\displaystyle\frac{\hat{T}_{v}(\lambda)}{T}=\frac{2}{Q_{0}}\hat{\Sigma}_{11}(\lambda)$
$\displaystyle=\frac{2}{Q_{0}}\int\frac{d\omega}{2\pi}\frac{\omega^{2}(\omega^{2}+\frac{n^{2}}{\tau^{2}})^{n}}{(-1)^{n}p_{H,\lambda}(s=i\omega)}$
$\displaystyle=\frac{2}{Q_{0}}\int\frac{d\omega}{2\pi}\frac{\omega^{2}}{|\chi_{n}(\omega)|^{-2}-\frac{4\lambda
g}{Q_{0}^{2}}\omega\sin_{n}(\omega\tau)}\ .$ (154b)

Likewise,

$\displaystyle\hat{\Sigma}_{n+2,1}(\lambda)=\langle
x_{n}(t)v(t)\rangle_{\lambda}$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{(\frac{n}{\tau})^{n}i\omega(i\omega+\frac{n}{\tau})^{n}}{(-1)^{n}p_{H}(\lambda,s=i\omega)}$
$\displaystyle=-\int\frac{d\omega}{2\pi}\frac{\omega\sin_{n}(\omega\tau)}{|\chi_{n}(\omega)|^{-2}-\frac{4\lambda
g}{Q_{0}^{2}}\omega\sin_{n}(\omega\tau)}\ ,$ (155)

and by comparing with the expression of $f^{+}(\lambda)$ [Eq. (124)], we find
that $(2g/Q_{0}^{2})\langle
x_{n}(t)v(t)\rangle_{\lambda}=-df^{+}(\lambda)/d\lambda$. This equality was
expected. It expresses that if one is interested by a particular fluctuation
in which the time-intensive observable
$(2g/Q_{0}^{2})(1/t)\int_{0}^{t}dt^{\prime}x_{n}(t^{\prime})\circ
dx(t^{\prime})$ takes the value $a$, this fluctuation is realized in the long-
time limit as a typical value in the effective process, with $\lambda$ given
by $a=-df^{+}(\lambda)/d\lambda$. (However, while the three observables ${\cal