June 02, 2026
The supplementary materials are devoted to provide detailed derivations in the main context.
In this section, we present detailed expressions of the linear response function under slow and weak driving.
Under slowly varying control protocols, linear response theory allows us to express the expectation value of a dynamical observable \(\langle A(t)\rangle\) in terms of response functions as \[\begin{align} \langle A(t)\rangle & \approx\langle A\rangle_{s}+\int_{-\infty}^{t}\boldsymbol{\chi}(A;t,t')\cdot[\boldsymbol{\lambda}(t)-\boldsymbol{\lambda}(t')]dt'\label{seq:linearresponse} \end{align}\tag{1}\] Expanding the control parameters around time \(t'\), \[\boldsymbol{\lambda}(t)\approx\boldsymbol{\lambda}(t')+\dot{\boldsymbol{\lambda}}(t')(t-t')+O((t-t')^{2}),\label{seq:protocolslow}\tag{2}\] and retaining terms up to first order, we obtain \[\langle A(t)\rangle\approx\langle A\rangle_{s}+\int_{-\infty}^{t}\boldsymbol{\chi}(A;t,t')\cdot\dot{\boldsymbol{\lambda}}(t')(t-t')dt'.\label{seq:linearapp}\tag{3}\] Here, the response function is defined as \(\boldsymbol{\chi}(A;t,t')\equiv\delta\langle A(t)\rangle/\delta\boldsymbol{\lambda}(t')\). Introducing the variable change \(t''=t-t'\), and expanding the control velocity as \[\dot{\boldsymbol{\lambda}}(t')=\dot{\boldsymbol{\lambda}}(t-t'')\approx\dot{\boldsymbol{\lambda}}(t)+O(\ddot{\boldsymbol{\lambda}}),\label{seq:protocolveco}\tag{4}\] we obtain a time-local approximation: \[\begin{align} \langle A(t)\rangle & \approx\langle A\rangle_{s}+\dot{\boldsymbol{\lambda}}(t)\cdot\int_{-\infty}^{t}\boldsymbol{\chi}(A;t,t')(t-t')dt'\nonumber \\ & =\langle A\rangle_{s}+\dot{\boldsymbol{\lambda}}(t)\cdot\int_{0}^{\infty}\boldsymbol{\chi}(A;t'',0)t''dt''\nonumber \\ & =\langle A\rangle_{s}+\boldsymbol{I}(A)\cdot\dot{\boldsymbol{\lambda}}(t),\label{seq:timelocallinear} \end{align}\tag{5}\] where \(\boldsymbol{I}(A)\equiv\int_{0}^{\infty}\boldsymbol{\chi}(A;t'',0)t''dt''\) represents the response kernel. This expression is equivalent to Eq.\(\;\)(5) in the main text.
The expectation value of a dynamical observable \(\langle A(t)\rangle\) can be expressed in the path-integral representation as \[\begin{align} \langle A(t)\rangle & =\Pi_{i}\iint D[\mathbf{v}_{i}]D[\boldsymbol{\xi}_{i}]A[t;\boldsymbol{\xi}_{i},\mathbf{v}_{i}]P_{v}[\mathbf{v}_{i}]P_{\xi}[\boldsymbol{\xi}_{i}],\label{seq:pathintegres} \end{align}\tag{6}\] where \(P_{\xi}[\boldsymbol{\xi}_{i}]\) and \(P_{v}[\mathbf{v}_{i}]\) represent the probability functionals of the thermal noise \(\boldsymbol{\xi}_{i}\) and the active noise \(\mathbf{v}_{i}\), respectively. Starting from the Langevin dynamics in Eq.\(\;\)(1), we insert the identity \[\begin{align} 1 & =\Pi_{i}\int D[\boldsymbol{\xi}_{i}]\delta(\dot{\mathbf{r}}_{i}+\gamma^{-1}\nabla_{i}U-\gamma^{-1}\boldsymbol{\xi}_{i}-\mathbf{v}_{i}),\label{seq:pathiden} \end{align}\tag{7}\] and transform the integration variables from \(\boldsymbol{\xi}_{i}\) to \(\mathbf{r}_{i}\), yielding
\[\begin{align} 1 & =\Pi_{i}\int D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]\delta(\dot{\mathbf{r}}_{i}+\gamma^{-1}\nabla_{i}U-\gamma^{-1}\boldsymbol{\xi}_{i}-\mathbf{v}_{i}),\label{seq:yakebi} \end{align}\tag{8}\] where \(J[\mathbf{r}_{i}]\equiv\mathrm{det}(\delta\boldsymbol{\xi}_{i}(t)/\delta\mathbf{r}_{i}(t'))\) is the corresponding Jacobian. Substituting this identity into the path integral and integrating out the thermal noise leads to \[\begin{align} \langle A(t)\rangle & =\Pi_{i}\iiint D[\mathbf{v}_{i}]D[\boldsymbol{\xi}_{i}]D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]A[t;\boldsymbol{\xi}_{i},\mathbf{v}_{i}]P_{v}[\mathbf{v}_{i}]P_{\xi}[\boldsymbol{\xi}_{i}]\delta(\dot{\mathbf{r}}_{i}+\gamma^{-1}\nabla_{i}U-\gamma^{-1}\boldsymbol{\xi}_{i}-\mathbf{v}_{i})\nonumber \\ & =\mathcal{N}^{-1}\Pi_{i}\iint D[\mathbf{v}_{i}]D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]P_{v}[\mathbf{v}_{i}]A[t;\mathbf{r}_{i},\mathbf{v}_{i}]e^{-\mathcal{S}[\mathbf{r}_{i},\mathbf{v}_{i}]},\label{seq:integthermalno} \end{align}\tag{9}\] where \(\mathcal{S}\equiv[\gamma/(4T)]\int dt(\dot{\mathbf{r}}_{i}+\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})^{2}\) denotes the Onsager-Machlup action [1], and \[\mathcal{N}\equiv\Pi_{i}\iint D[\mathbf{v}_{i}]D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]P_{v}[\mathbf{v}_{i}]e^{-\mathcal{S}[\mathbf{r}_{i},\mathbf{v}_{i}]}\label{seq:normconst}\tag{10}\] ensures normalization.
Differentiating the path-integral expression in Eq.\(\;\)(9 ) yields the response function as \[\begin{align} \boldsymbol{\chi}(A;t,t') & \equiv\frac{\delta\langle A(t)\rangle}{\delta\boldsymbol{\lambda}(t')}\nonumber \\ & =\Pi_{i}\iint D[\mathbf{v}_{i}]D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]P_{v}[\mathbf{v}_{i}]A[t;\mathbf{r}_{i},\mathbf{v}_{i}]\frac{\delta(\mathcal{N}^{-1}e^{-\mathcal{S}})}{\delta\boldsymbol{\lambda}(t')}\nonumber \\ & =-\Pi_{i}\iint D[\mathbf{v}_{i}]D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]P_{v}[\mathbf{v}_{i}]A[t;\mathbf{r}_{i},\mathbf{v}_{i}]\mathcal{N}^{-1}e^{-\mathcal{S}}(\frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')}-\langle\frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')}\rangle)\nonumber \\ & =-\langle A(t)\frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')}\rangle+\langle A(t)\rangle\langle\frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')}\rangle,\label{seq:resfunS} \end{align}\tag{11}\] where we have used \[\begin{align} \frac{\delta\mathcal{N}}{\delta\boldsymbol{\lambda}(t')} & =-\Pi_{i}\iint D[\mathbf{v}_{i}]D[\mathbf{r}_{i}]J[\mathbf{r}_{i}]P_{v}[\mathbf{v}_{i}]e^{-\mathcal{S}}\frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')}\nonumber \\ & =-\mathcal{N}\langle\frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')}\rangle.\label{seq:deffalphanoe} \end{align}\tag{12}\] For the unperturbed dynamics with constant control parameters, the functional derivative of the action becomes \[\begin{align} \frac{\delta\mathcal{S}}{\delta\boldsymbol{\lambda}(t')} & =\frac{1}{2T}\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})(\dot{\mathbf{r}}_{i}+\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})\nonumber \\ & =\frac{1}{2T}[\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})\cdot\dot{\mathbf{r}}_{i}+\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})\cdot(\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})]\nonumber \\ & =\frac{1}{2T}[\frac{d}{dt'}(\frac{\partial U}{\partial\boldsymbol{\lambda}})+\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})\cdot(\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})].\label{seq:deffsrespon} \end{align}\tag{13}\] Substituting this result into the response function finally gives \[\begin{align} \boldsymbol{\chi}(A;t,t') & =\frac{1}{2T}[-\frac{d}{dt'}\langle A(t)(\frac{\partial U}{\partial\boldsymbol{\lambda}})_{t'}\rangle+\langle A\rangle\frac{d}{dt'}\langle(\frac{\partial U}{\partial\boldsymbol{\lambda}})_{t'}\rangle-\langle A(t)\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})_{t'}\cdot(\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})_{t'}\rangle+\langle A\rangle\langle\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})\cdot(\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})\rangle].\label{seq:finalresponse} \end{align}\tag{14}\] For a single control parameter, \(\boldsymbol{\lambda}\to\lambda\), this expression reduces to the response formula derived by Davis et al. [2]. Equation\(\;\)(14 ) provides a generalized fluctuation–response relation for active matter under multidimensional control protocols and serves as the starting point for the thermodynamic geometric formulation developed below. The response kernel then follows as \[\begin{align} \boldsymbol{I}(A)\equiv & \int_{0}^{\infty}\boldsymbol{\chi}(A;t,0)tdt\nonumber \\ = & \frac{1}{2T}\int_{0}^{\infty}\{\langle A(t)(\frac{\partial U}{\partial\boldsymbol{\lambda}})_{t'=0}\rangle-\langle A(t)\rangle\langle(\frac{\partial U}{\partial\boldsymbol{\lambda}})_{t'=0}\rangle-t[\langle A(t)\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})_{t'=0}\cdot(\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})_{t'=0}\rangle\nonumber \\ & -\langle A\rangle\langle\nabla_{i}(\frac{\partial U}{\partial\boldsymbol{\lambda}})\cdot(\gamma^{-1}\nabla_{i}U-\mathbf{v}_{i})\rangle]\}.\label{seq:resknerl} \end{align}\tag{15}\]
According to stochastic thermodynamics [3], [4], the mean output work is defined as \[W\equiv-\int_{0}^{\tau}\langle\partial U/\partial\boldsymbol{\lambda}\rangle\cdot\dot{\boldsymbol{\lambda}}dt\label{seq:meanwork}\tag{16}\] and the mean heat absorbed from the environment is \[Q\equiv\sum_{i}\int_{0}^{\tau}\langle\dot{\mathbf{r}}_{i}\circ(-\gamma\dot{\mathbf{r}}_{i}+\boldsymbol{\xi}_{i})\rangle dt.\label{seq:meanheat}\tag{17}\] Using the first-law relation, the heat can be rewritten as \[Q=\Delta\langle U\rangle-\int_{0}^{\tau}(\langle\partial U/\partial\boldsymbol{\lambda}\rangle\cdot\dot{\boldsymbol{\lambda}}+\langle f\rangle)dt,\label{seq:firstlawre}\tag{18}\] where \(\Delta\langle U\rangle\equiv\langle U(\tau)\rangle-\langle U(0)\rangle_{s}\) denotes the mean potential difference. Under the linear response approximation, the relevant observables in Eqs.\(\;\)(16 ) and\(\;\)(17 ) can be expanded as \[\begin{align} \langle\partial U/\partial\boldsymbol{\lambda}\rangle & \approx\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}+\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot\dot{\boldsymbol{\lambda}}(t),\nonumber \\ \langle U(\tau)\rangle & \approx\langle U(\tau)\rangle_{s}+\boldsymbol{I}(U)\cdot\dot{\boldsymbol{\lambda}}(\tau),\nonumber \\ \langle f\rangle & \approx\langle f\rangle_{s}+\boldsymbol{I}(f)\cdot\dot{\boldsymbol{\lambda}}(t),\label{seq:linearresobser} \end{align}\tag{19}\] where \(\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\), \(\boldsymbol{I}(U)\), and \(\boldsymbol{I}(f)\) denote the corresponding response kernels. Substituting these expressions into the definitions of work and heat, we obtain \[W\approx-\int_{0}^{\tau}\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}\cdot\dot{\boldsymbol{\lambda}}dt-\int_{0}^{\tau}\dot{\boldsymbol{\lambda}}^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot\dot{\boldsymbol{\lambda}}dt\label{seq:linearmeanwork}\tag{20}\] and \[\begin{align} Q\approx & \boldsymbol{I}(U)\cdot\dot{\boldsymbol{\lambda}}(\tau)-\int_{0}^{\tau}\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}\cdot\dot{\boldsymbol{\lambda}}dt-\int_{0}^{\tau}\boldsymbol{I}(f)\cdot\dot{\boldsymbol{\lambda}}dt\nonumber \\ & -\int_{0}^{\tau}\dot{\boldsymbol{\lambda}}^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot\dot{\boldsymbol{\lambda}}dt-\int_{0}^{\tau}\langle f\rangle_{s}dt,\label{seq:linearmeanheat} \end{align}\tag{21}\] where \(\Delta\langle U\rangle_{s}=0\) for cyclic processes in the linear response approxiamtion. The results in Eqs.\(\;\)(20 ) and\(\;\)(21 ) directly correspond to Eqs\(\;\)(6) and\(\;\)(7) in the main text.
In this section, we prove the scaling relations of the mean work and heat within the linear response framework.
In linear response theory, the response kernel associated with an observable \(A\) is defined as \[\begin{align} \boldsymbol{I}(A)\equiv\int_{0}^{\infty}\boldsymbol{\chi}(A;t'',0)t'dt', & \;\text{\;\;}\boldsymbol{\chi}(A;t,0)\equiv\delta\langle A(t)\rangle/\delta\boldsymbol{\lambda}(0).\label{seq:responsekernel} \end{align}\tag{22}\] By construction, the response function \(\boldsymbol{\chi}\) characterizes the system’s response to an infinitesimal perturbation around a reference state with fixed control parameters. As a result, the kernel \(\boldsymbol{I}(A)\) depends only on the intrinsic dynamical properties of the system—such as the interaction potential, noise statistics, and activity parameters—and is independent of the driving rate \(\dot{\boldsymbol{\lambda}}\).
For the commonly adopted protocol design \(\boldsymbol{\lambda}(t)=\boldsymbol{\lambda}(t/\tau)\) [5], introducing the rescaled time \(u=t/\tau\), the control velocity scales as \[\dot{\boldsymbol{\lambda}}(t)=\frac{1}{\tau}\boldsymbol{\lambda}'(u),\text{\;\;\text{\;}}\boldsymbol{\lambda}'\equiv d\boldsymbol{\lambda}/du.\label{seq:velocityscaling}\tag{23}\] Accordingly, the linear-response correction behaves as \[\boldsymbol{I}(A)\cdot\dot{\boldsymbol{\lambda}}(t)=\frac{1}{\tau}\boldsymbol{I}(A)\cdot[\boldsymbol{\lambda}'](u),\label{seq:responsescaling}\tag{24}\] and its time integral satisfies \[\int_{0}^{\tau}\boldsymbol{I}(A)\cdot\dot{\boldsymbol{\lambda}}dt=\int_{0}^{1}\boldsymbol{I}(A)\cdot[\boldsymbol{\lambda}']du,\label{seq:integralscaling}\tag{25}\] which is independent of the protocol duration \(\tau\).
In contrast, the steady-state contribution \(\langle A\rangle_{s}\), evaluated at fixed control parameters, is independent of the driving rate. Its time integral therefore scales as \[\int_{0}^{\tau}\langle A\rangle_{s}dt=\tau\int_{0}^{1}\langle A\rangle_{s}du,\label{seq:tauplscaling}\tag{26}\] provided that \(\langle A\rangle_{s}\) remains finite along the protocol.
Physically, this separation reflects two distinct contributions: the steady-state term \(\langle A\rangle_{s}\) describes persistent nonequilibrium activity and gives rise to an extensive \(O(\tau)\) contribution, while the kernel \(\boldsymbol{I}(A)\) plays a role analogous to a transport coefficient, characterizing the intrinsic response of the active system to external perturbations and yielding a finite \(O(1)\) correction, while the driving rate \(\dot{\boldsymbol{\lambda}}\) acts as the conjugate thermodynamic force.
Based on the above analysis, the response relations in Eq.\(\;\)(19 ) take the form \[\begin{align} \langle\partial U/\partial\boldsymbol{\lambda}\rangle & \approx\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}+\frac{1}{\tau}\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot[\boldsymbol{\lambda}'](u),\nonumber \\ \langle U(\tau)\rangle & \approx\langle U(1)\rangle_{s}+\frac{1}{\tau}\boldsymbol{I}(U)\cdot[\boldsymbol{\lambda}'](1),\nonumber \\ \langle f\rangle & \approx\langle f\rangle_{s}+\frac{1}{\tau}\boldsymbol{I}(f)\cdot[\boldsymbol{\lambda}'](u).\label{seq:resuufscal} \end{align}\tag{27}\] Substituting these expressions into the approximate forms of work and heat in Eqs.\(\;\)(20 ) and\(\;\)(21 ), we obtain \[W\approx-\int_{0}^{1}\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}\cdot[\boldsymbol{\lambda}']du-\frac{1}{\tau}\int_{0}^{1}[\boldsymbol{\lambda}']^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot[\boldsymbol{\lambda}']du\label{seq:scalingstructurework}\tag{28}\] and \[\begin{align} Q\approx & -\int_{0}^{1}\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}\cdot[\boldsymbol{\lambda}']dt-\int_{0}^{1}\boldsymbol{I}(f)\cdot[\boldsymbol{\lambda}']du\nonumber \\ & -\tau\int_{0}^{1}\langle f\rangle_{s}du+\frac{1}{\tau}\{\boldsymbol{I}(U)\cdot[\boldsymbol{\lambda}'](1)-\int_{0}^{1}[\boldsymbol{\lambda}']^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot[\boldsymbol{\lambda}']du\}.\label{seq:scalingstructureheat} \end{align}\tag{29}\] These expressions reveal a universal scaling structure characterized by three distinct contributions: \[W=\Gamma-\frac{\Sigma}{\tau},\;\;\;Q\approx\Lambda-\frac{\Sigma}{\tau}-\Omega\tau,\label{seq:universcali}\tag{30}\] where the \(\tau\)-independent coefficients are \[\begin{gather} \Gamma\equiv-\int_{0}^{1}\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}\cdot[\boldsymbol{\lambda}']du,\;\;\;\boldsymbol{\lambda}'\equiv d\boldsymbol{\lambda}/du,\nonumber \\ \Sigma\equiv\int_{0}^{1}[\boldsymbol{\lambda}']^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot[\boldsymbol{\lambda}']du,\;\;\;\Omega\equiv\int_{0}^{1}\langle f\rangle_{s}du,\nonumber \\ \Lambda\equiv-\int_{0}^{1}\langle\partial U/\partial\boldsymbol{\lambda}\rangle_{s}\cdot[\boldsymbol{\lambda}']du-\int_{0}^{1}\boldsymbol{I}(f)\cdot[\boldsymbol{\lambda}']du.\label{eq:functiongasg} \end{gather}\tag{31}\] Here we have considered the fact that the boundary term \(\boldsymbol{I}(U)\cdot[\boldsymbol{\lambda}'](1)\) is subleading compared with the integrated dissipation term and can therefore be neglected for sufficiently large protocol duration \(\tau\). The leading \(O(1/\tau)\) contribution then arises from the quadratic form \[\Sigma\equiv\int_{0}^{1}[\boldsymbol{\lambda}']^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot[\boldsymbol{\lambda}']du,\label{seq:geometricdiss}\tag{32}\] which defines a geometric dissipation functional. The universal scaling in Eq.\(\;\)(30 ) reduces to Eq.\(\;\)(8) in the main text.
In this section, we consider a two-dimensional active Brownian particle, whose position vector is denoted by \(\mathbf{r}=(x,y)\), confined by the harmonic potential \[U(x,y,\boldsymbol{\lambda})=\frac{\lambda_{1}}{2}(x^{2}+y^{2})+\lambda_{2}xy,\label{seq:harmoncon}\tag{33}\] where \(\boldsymbol{\lambda}(t)\equiv(\lambda_{1},\lambda_{2})\) represents the set of externally controlled parameters. Stability of the confinement requires \(\lambda_{1}>|\lambda_{2}|\). The dynamics of the particle is governed by the Langevin equations \[\begin{align} \dot{x} & =-\lambda_{1}x-\lambda_{2}y+v\cos\theta+\xi_{x}(t),\nonumber \\ \dot{y} & =-\lambda_{2}x-\lambda_{1}y+v\sin\theta+\xi_{y}(t).\label{seq:langevinequhar} \end{align}\tag{34}\] Throughout this example, we set the friction coefficient \(\gamma=1\) for simplicity.
To evaluate the correlation functions entering the response kernels, we first solve the Langevin equations\(\;\)(34 ). Introducing the vector notation \(\mathbf{e}=(\cos\theta\text{\;}\sin\theta)^{T}\), the dynamics can be written compactly as \[\dot{\mathbf{r}}=-A\mathbf{r}+v\mathbf{e}+\boldsymbol{\xi},\label{seq:vectorlangevin}\tag{35}\] where \[A=\left(\begin{array}{cc} \lambda_{1} & \lambda_{2}\\ \lambda_{2} & \lambda_{1} \end{array}\right).\label{seq:harstrength}\tag{36}\] For fixed control parameters \(\boldsymbol{\lambda}\), the formal solution is \[\mathbf{r}(t)=\mathbf{r}(0)e^{-A(t-t_{0})}+\int_{t_{0}}^{t}dse^{-A(t-s)}(v\mathbf{e}(s)+\boldsymbol{\xi}(s)).\label{seq:formalsoultion}\tag{37}\] Taking the stationary limit \(t_{0}\to-\infty\), the contribution from the initial condition decays exponentially and vanishes, yielding \[\mathbf{r}(t)=\int_{t_{0}}^{t}dse^{-A(t-s)}(v\mathbf{e}(s)+\boldsymbol{\xi}(s)).\label{seq:stationlimit}\tag{38}\] The matrix \(A\) is diagonalizable with eigenvalues \(\alpha_{\pm}=\lambda_{1}\pm|\lambda_{2}|\), and can be decomposed as \[A=P\left(\begin{array}{cc} \alpha_{+} & 0\\ 0 & \alpha_{-} \end{array}\right)P^{T},\label{seq:Adecomn}\tag{39}\] where \[P=\frac{1}{\sqrt{2}}\left(\begin{array}{cc} 1 & 1\\ 1 & -1 \end{array}\right).\label{seq:eigmatrix}\tag{40}\] Using the spectral mapping theorem, the matrix exponential becomes \[e^{-At}=P\left(\begin{array}{cc} e^{-\alpha_{+}t} & 0\\ 0 & e^{-\alpha_{-}t} \end{array}\right)P^{T}=\frac{1}{2}\left(\begin{array}{cc} e^{-\alpha_{+}t}+e^{-\alpha_{+}t} & e^{-\alpha_{+}t}-e^{-\alpha_{+}t}\\ e^{-\alpha_{+}t}-e^{-\alpha_{+}t} & e^{-\alpha_{+}t}+e^{-\alpha_{+}t} \end{array}\right),\label{seq:ematrixdecom}\tag{41}\] which can be written explicitly as \[e^{-At}=\frac{1}{2}\left(\begin{array}{cc} e^{-\alpha_{+}t}+e^{-\alpha_{+}t} & e^{-\alpha_{+}t}-e^{-\alpha_{+}t}\\ e^{-\alpha_{+}t}-e^{-\alpha_{+}t} & e^{-\alpha_{+}t}+e^{-\alpha_{+}t} \end{array}\right).\label{seq:eacorrfunex}\tag{42}\]
Using the solution obtained above, the correlation functions required for evaluating the response kernels can be derived analytically.
The coordinate correlation functions are \[\begin{align} C_{xx}(t) & =\langle x(t)x(0)\rangle=\frac{K_{+}}{4(\lambda_{1}+|\lambda_{2}|)}e^{-\alpha_{+}t}+\frac{K_{-}}{4(\lambda_{1}-|\lambda_{2}|)}e^{-\alpha_{-}t},\nonumber \\ C_{xy}(t) & =\langle x(t)y(0)\rangle=\frac{K_{+}}{4(\lambda_{1}+|\lambda_{2}|)}e^{-\alpha_{+}t},\nonumber \\ C_{yx}(t) & =\langle y(t)x(0)\rangle=\frac{K_{+}}{4(\lambda_{1}+|\lambda_{2}|)}e^{-\alpha_{+}t},\nonumber \\ C_{yy}(t) & =\langle y(t)y(0)\rangle=\frac{K_{+}}{4(\lambda_{1}+|\lambda_{2}|)}e^{-\alpha_{+}t}+\frac{K_{-}}{4(\lambda_{1}-|\lambda_{2}|)}e^{-\alpha_{-}t},\label{seq:coordinatecorr} \end{align}\tag{43}\] where \(K_{\pm}\equiv2T+v^{2}\tau_{p}/(1-\tau_{p}^{2}\alpha_{\pm}^{2})\). The correlations between the particle coordinates and the self-propulsion velocity are \[\begin{align} D_{xx}(t) & =v\langle x(t)e_{x}(0)\rangle=\frac{R_{+}}{2}e^{-\alpha_{+}t}+\frac{R_{-}}{2}e^{-\alpha_{-}t},\nonumber \\ D_{xy}(t) & =v\langle x(t)e_{y}(0)\rangle=\frac{R_{+}}{2}e^{-\alpha_{+}t}-\frac{R_{-}}{2}e^{-\alpha_{-}t},\nonumber \\ D_{yx}(t) & =v\langle y(t)e_{x}(0)\rangle=\frac{R_{+}}{2}e^{-\alpha_{+}t}-\frac{R_{-}}{2}e^{-\alpha_{-}t},\nonumber \\ D_{yy}(t) & =v\langle y(t)e_{y}(0)\rangle=\frac{R_{+}}{2}e^{-\alpha_{+}t}+\frac{R_{-}}{2}e^{-\alpha_{-}t},\label{seq:xvcorr} \end{align}\tag{44}\] with \(R_{\pm}\equiv v^{2}\tau_{p}/(1-\tau_{p}^{2}\alpha_{\pm}^{2})\). The correlations between the self-propulsion velocity and the particle coordinates are \[\begin{align} E_{xx}(t) & =v\langle e_{x}(t)x(0)\rangle=\frac{v^{2}}{4}(\frac{1}{\alpha_{+}}+\frac{1}{\alpha_{-}})e^{-t/\tau_{p}},\nonumber \\ E_{xy}(t) & =v^{2}\langle e_{x}(t)y(0)\rangle=\frac{v^{2}}{4}(\frac{1}{\alpha_{+}}-\frac{1}{\alpha_{-}})e^{-t/\tau_{p}},\nonumber \\ E_{yx}(t) & =v^{2}\langle e_{y}(t)x(0)\rangle=\frac{v^{2}}{4}(\frac{1}{\alpha_{+}}-\frac{1}{\alpha_{-}})e^{-t/\tau_{p}},\nonumber \\ E_{yy}(t) & =v^{2}\langle e_{y}(t)y(0)\rangle=\frac{v^{2}}{4}(\frac{1}{\alpha_{+}}+\frac{1}{\alpha_{-}})e^{-t/\tau_{p}}.\label{seq:vxcorre} \end{align}\tag{45}\] Finally, the correlations between the thermal noise and the particle coordinates are \[\begin{align} F_{xx}(t) & =\langle\xi_{x}(t)x(0)\rangle=T(\frac{1}{\alpha_{+}}+\frac{1}{\alpha_{-}})\delta(t),\nonumber \\ F_{xy}(t) & =\langle\xi_{x}(t)y(0)\rangle=T(\frac{1}{\alpha_{+}}-\frac{1}{\alpha_{-}})\delta(t),\nonumber \\ F_{yx}(t) & =\langle\xi_{y}(t)x(0)\rangle=T(\frac{1}{\alpha_{+}}-\frac{1}{\alpha_{-}})\delta(t),\nonumber \\ F_{yy}(t) & =\langle\xi_{y}(t)y(0)\rangle=T(\frac{1}{\alpha_{+}}+\frac{1}{\alpha_{-}})\delta(t).\label{seq:thenoxcorr} \end{align}\tag{46}\]
Substituting the correlation functions derived in the previous subsection into Eq.\(\;\)(15 ), the response kernels associated with the generalized force \(\partial U/\partial\boldsymbol{\lambda}\) can be expressed as \[\begin{align} I_{1}(\partial U/\partial\lambda_{1})= & \frac{1}{2T}\int_{0}^{\infty}dt[(1+2\lambda_{1}t)(C_{xx}^{2}(t)+C_{xy}^{2}(t)+C_{yx}^{2}(t)+C_{yy}^{2}(t))+2\lambda_{2}t(C_{xx}(t)C_{xy}(t)+C_{yx}(t)C_{yy}(t))\nonumber \\ & -t(C_{xx}(t)D_{xx}(t)+C_{xy}(t)D_{xy}(t)+C_{yx}(t)D_{yx}(t)+C_{yy}(t)D_{yy}(t))],\nonumber \\ I_{1}(\partial U/\partial\lambda_{2})= & \frac{1}{2T}\int_{0}^{\infty}dt[(1+2\lambda_{1}t)(C_{xx}(t)C_{xy}(t)+C_{yx}(t)C_{yy}(t))+\lambda_{2}t(C_{xx}^{2}(t)+C_{xy}^{2}(t)+C_{yx}^{2}(t)+C_{yy}^{2}(t))\nonumber \\ & -t(C_{xy}(t)D_{xx}(t)+C_{xx}(t)D_{xy}(t)+C_{yy}(t)D_{yx}(t)+C_{yx}(t)D_{yy}(t))],\nonumber \\ I_{2}(\partial U/\partial\lambda_{1})= & \frac{1}{2T}\int_{0}^{\infty}dt[(1+2\lambda_{1}t)(C_{xx}(t)C_{yx}(t)+C_{xy}(t)C_{yy}(t))+2\lambda_{2}t(C_{xx}(t)C_{yy}(t)+C_{xy}(t)C_{yx}(t))\nonumber \\ & -t(C_{xx}(t)D_{yx}(t)+C_{yx}(t)D_{xx}(t)+C_{xy}(t)D_{yy}(t)+C_{yy}(t)D_{xy}(t))],\nonumber \\ I_{2}(\partial U/\partial\lambda_{2})= & \frac{1}{2T}\int_{0}^{\infty}dt[(1+2\lambda_{1}t)(C_{xx}(t)C_{yy}(t)+C_{xy}(t)C_{yx}(t))+2\lambda_{2}t(C_{xx}(t)C_{yx}(t)+C_{xy}(t)C_{yy}(t))\nonumber \\ & -t(C_{xy}(t)D_{yx}(t)+C_{yy}(t)D_{xx}(t)+C_{xx}(t)D_{yy}(t)+C_{yx}(t)D_{xy}(t))].\label{seq:IUlambdacorre} \end{align}\tag{47}\] Because all correlation functions consist of exponential modes, the above time integrals can be evaluated analytically. After straightforward calculations, we obtain \[\begin{align} I_{1}(\partial U/\partial\lambda_{1})= & \frac{1}{2T}\{\frac{K_{+}^{2}[1+(\lambda_{1}+\lambda_{2}/2)/(\lambda_{1}+|\lambda_{2}|)]}{8(\lambda_{1}+|\lambda_{2}|)^{3}}+\frac{K_{-}^{2}[1+\lambda_{1}/(\lambda_{1}-|\lambda_{2}|)]}{16(\lambda_{1}-|\lambda_{2}|)^{3}}+\frac{K_{+}K_{-}[1+\lambda_{2}/(4\lambda_{1})]}{4\lambda_{1}(\lambda_{1}^{2}-\lambda_{2}^{2})}\nonumber \\ & -\frac{K_{+}R_{+}}{8(\lambda_{1}+|\lambda_{2}|)^{3}}-\frac{K_{-}R_{-}}{16(\lambda_{1}-|\lambda_{2}|)^{3}}-\frac{K_{-}R_{+}}{16\lambda_{1}^{2}(\lambda_{1}-|\lambda_{2}|)}\},\nonumber \\ I_{1}(\partial U/\partial\lambda_{2})= & \frac{1}{2T}\{\frac{K_{+}^{2}[1+(\lambda_{1}+\lambda_{2})/(\lambda_{1}+|\lambda_{2}|)]}{16(\lambda_{1}+|\lambda_{2}|)^{3}}+\frac{K_{+}K_{-}[1+\lambda_{2}/(2\lambda_{1})]}{8\lambda_{1}(\lambda_{1}^{2}-\lambda_{2}^{2})}+\frac{K_{-}^{2}\lambda_{2}}{32(\lambda_{1}-|\lambda_{2}|)^{4}}\nonumber \\ & -\frac{K_{+}R_{+}}{8(\lambda_{1}+|\lambda_{2}|)^{3}}+\frac{K_{-}R_{-}}{16(\lambda_{1}-|\lambda_{2}|)^{3}}-\frac{K_{-}R_{+}}{16\lambda_{1}^{2}(\lambda_{1}-|\lambda_{2}|)}\},\nonumber \\ I_{2}(\partial U/\partial\lambda_{1})= & \frac{1}{2T}\{\frac{K_{+}^{2}[1+(\lambda_{1}+\lambda_{2})/(\lambda_{1}+|\lambda_{2}|)]}{16(\lambda_{1}+|\lambda_{2}|)^{3}}+\frac{K_{+}K_{-}[1+\lambda_{2}/(2\lambda_{1})]}{8\lambda_{1}(\lambda_{1}^{2}-\lambda_{2}^{2})}+\frac{K_{-}^{2}\lambda_{2}}{32(\lambda_{1}-|\lambda_{2}|)^{4}}\nonumber \\ & -\frac{K_{+}R_{+}}{8(\lambda_{1}+|\lambda_{2}|)^{3}}+\frac{K_{-}R_{-}}{16(\lambda_{1}-|\lambda_{2}|)^{3}}-\frac{K_{-}R_{+}}{16\lambda_{1}^{2}(\lambda_{1}-|\lambda_{2}|)}\},\nonumber \\ I_{2}(\partial U/\partial\lambda_{2})= & \frac{1}{2T}\{\frac{K_{+}^{2}[1+(\lambda_{1}+\lambda_{2})/(\lambda_{1}+|\lambda_{2}|)]}{16(\lambda_{1}+|\lambda_{2}|)^{3}}+\frac{K_{-}^{2}[1+\lambda_{1}/(\lambda_{1}-|\lambda_{2}|)]}{32(\lambda_{1}-|\lambda_{2}|)^{3}}+\frac{K_{+}K_{-}[1+\lambda_{2}/(2\lambda_{1})]}{8\lambda_{1}(\lambda_{1}^{2}-\lambda_{2}^{2})}\nonumber \\ & -\frac{K_{+}R_{+}}{8(\lambda_{1}+|\lambda_{2}|)^{3}}-\frac{K_{-}R_{-}}{16(\lambda_{1}-|\lambda_{2}|)^{3}}-\frac{K_{-}R_{+}}{16\lambda_{1}^{2}(\lambda_{1}-|\lambda_{2}|)}\}.\label{seq:IUlambdaan} \end{align}\tag{48}\] These expressions completely determine the dissipation metric and yield the rescaled dissipation \(\Sigma\equiv\int_{0}^{1}[\boldsymbol{\lambda}']^{T}\cdot\boldsymbol{I}(\partial U/\partial\boldsymbol{\lambda})\cdot[\boldsymbol{\lambda}']du\). Applying the same procedure to the active-force contribution \(f\), we obtain the response kernels \(\boldsymbol{I}(f)\) as \[\begin{align} I_{1}(f)= & \frac{1}{2T}\{-\frac{v^{2}}{2}[\frac{K_{+}}{2(\lambda_{1}+|\lambda_{2}|)(\alpha_{+}+1/\tau_{p})}+\frac{\lambda_{1}^{2}+\lambda_{2}|\lambda_{2}|}{\lambda_{1}+|\lambda_{2}|}\frac{\lambda_{1}K_{-}}{2(\lambda_{1}-|\lambda_{2}|)^{2}(\alpha_{-}+1/\tau_{p})}]\nonumber \\ & -\frac{v^{2}}{2}[\frac{(\lambda_{1}+\lambda_{2})^{2}K_{+}}{(\lambda_{1}+|\lambda_{2}|)^{2}(\alpha_{+}+1/\tau_{p})^{2}}+\frac{\lambda_{1}}{\lambda_{1}+|\lambda_{2}|}\frac{(\lambda_{1}^{2}+\lambda_{2}^{2}-2\lambda_{2}|\lambda_{2}|)K_{-}}{(\lambda_{1}-|\lambda_{2}|)^{2}(\alpha_{-}+1/\tau_{p})^{2}}]\nonumber \\ & +\frac{\lambda_{1}^{2}+\lambda_{2}^{2}}{(\lambda_{1}^{2}-\lambda_{2}^{2})^{2}}\frac{v^{4}\tau_{p}}{4}(1+\lambda_{1}\tau_{p})-\frac{\lambda_{1}|\lambda_{2}|}{(\lambda_{1}^{2}-\lambda_{2}^{2})^{2}}\frac{v^{4}\tau_{p}^{2}\lambda_{2}}{2}+\frac{\lambda_{1}^{2}+\lambda_{2}^{2}}{(\lambda_{1}^{2}-\lambda_{2}^{2})^{2}}2v^{2}D\nonumber \\ & +\frac{v^{2}}{2}[\frac{(\lambda_{1}+\lambda_{2})(K_{+}+2R_{+})}{2(\lambda_{1}+|\lambda_{2}|)(\alpha_{+}+1/\tau_{p})^{2}}+\frac{\lambda_{1}K_{-}+2(\lambda_{1}-\lambda_{2})R_{-}}{2(\lambda_{1}-|\lambda_{2}|)(\alpha_{-}+1/\tau_{p})^{2}}]-\frac{\lambda_{1}}{\lambda_{1}^{2}-\lambda_{2}^{2}}\frac{v^{4}\tau_{p}^{2}}{4}\},\nonumber \\ I_{2}(f)= & \frac{1}{2T}\{-\frac{v^{2}}{2}[\frac{(\lambda_{1}+\lambda_{2})K_{+}}{2(\lambda_{1}+|\lambda_{2}|)^{2}(\alpha_{+}+1/\tau_{p})}-\frac{\lambda_{1}}{\lambda_{1}+|\lambda_{2}|}\frac{(|\lambda_{2}|-\lambda_{2})K_{-}}{2(\lambda_{1}-|\lambda_{2}|)^{2}(\alpha_{-}+1/\tau_{p})}]\nonumber \\ & -\frac{v^{2}}{2}[\frac{(\lambda_{1}+\lambda_{2})^{2}K_{+}}{(\lambda_{1}+|\lambda_{2}|)^{2}(\alpha_{+}+1/\tau_{p})^{2}}-\frac{1}{\lambda_{1}+|\lambda_{2}|}\frac{(\lambda_{1}^{2}|\lambda_{2}|+\lambda_{2}^{2}|\lambda_{2}|-2\lambda_{1}^{2}\lambda_{2})K_{-}}{(\lambda_{1}-|\lambda_{2}|)^{2}(\alpha_{-}+1/\tau_{p})^{2}}]\nonumber \\ & -\frac{\lambda_{1}|\lambda_{2}|}{(\lambda_{1}^{2}-\lambda_{2}^{2})^{2}}\frac{v^{4}\tau_{p}}{2}(1+\lambda_{1}\tau_{p})+\frac{\lambda_{1}^{2}+\lambda_{2}^{2}}{(\lambda_{1}^{2}-\lambda_{2}^{2})^{2}}\frac{v^{4}\tau_{p}^{2}\lambda_{2}}{4}-\frac{2\lambda_{1}|\lambda_{2}|}{(\lambda_{1}^{2}-\lambda_{2}^{2})^{2}}2v^{2}D\nonumber \\ & +\frac{v^{2}}{2}[\frac{(\lambda_{1}+\lambda_{2})K_{+}+4\lambda_{1}R_{+}}{2(\lambda_{1}+|\lambda_{2}|)(\alpha_{+}+1/\tau_{p})^{2}}+\frac{\lambda_{2}K_{-}}{2(\lambda_{1}-|\lambda_{2}|)(\alpha_{-}+1/\tau_{p})^{2}}]-\frac{|\lambda_{2}|}{\lambda_{1}^{2}-\lambda_{2}^{2}}\frac{v^{4}\tau_{p}^{2}}{4}\},\label{seq:Ifreskern} \end{align}\tag{49}\] which determine the excess active energy input generated by a finite-time protocol, \(\Psi\equiv\int_{0}^{1}\boldsymbol{I}(f)\cdot[\boldsymbol{\lambda}']du\).
The quasistatic output work can likewise be evaluated analytically. Using Stokes’ theorem, it can be expressed as \(\Gamma=-\sum_{\mu}\oint_{\mathcal{C}}\langle\partial_{\mu}U\rangle_{s}d\lambda_{\mu}=\sum_{\mu\nu}\iint\mathcal{F}_{\mu\nu}d\lambda_{\mu}\wedge d\lambda_{\nu}\), where the only nonvanishing component of the thermodynamic curvature is \[\mathcal{F}_{12}=\frac{K_{+}[1-sgn(\lambda_{2})]}{4(\lambda_{1}+|\lambda_{2}|)^{2}}-\frac{\tau_{p}}{2v^{2}}K_{+}R_{+}[1-sgn(\lambda_{2})]+\frac{K_{-}sgn(\lambda_{2})}{4(\lambda_{1}-|\lambda_{2}|)^{2}}-\frac{\tau_{p}}{2v^{2}}K_{-}R_{-}sgn(\lambda_{2}).\label{seq:fcurveture}\tag{50}\] Finally, the steady-state energy injection rate due to self-propulsion is given by \(\Omega\equiv\int_{0}^{1}\langle f\rangle_{s}du\) with \[\begin{align} \langle f\rangle_{s} & =-\lambda_{1}(R_{+}+R_{-})-\lambda_{2}(R_{+}-R_{-})+v^{2}.\label{seq:averf} \end{align}\tag{51}\] Together with the dissipation metric and thermodynamic curvature derived above, these analytical results completely characterize the geometric thermodynamics of the harmonically confined active Brownian particle.