Impact of periodic thermal driving on heat fluctuations in a harmonic system

TL;DR

This paper investigates how sinusoidal thermal driving influences heat fluctuations in a harmonic oscillator, revealing resonance effects, non-Gaussian heat distributions, and the importance of inertia in non-equilibrium steady states.

Contribution

It provides analytical and numerical insights into the effects of periodic thermal driving on heat fluctuations and system dynamics in a harmonic system, highlighting inertia's role.

Findings

Resonant position-velocity correlations emerge under thermal driving.

Heat distribution becomes non-Gaussian and asymmetric, obeying Fluctuation Theorem.

Inertia dampens response and slows relaxation to steady state.

Abstract

The thermodynamics of mesoscopic systems driven by time-varying temperatures is crucial for understanding biological systems, designing nanoscale engines, and performing micro-particle cooling. In this work, we analyze an underdamped Brownian particle in a harmonic trap under a sinusoidal thermal protocol. Through analytical methods and numerical simulations, we analyze the system's dynamics and heat statistics. We report the emergence of resonant position-velocity correlations and a non-Gaussian, asymmetric heat distribution consistent with the Fluctuation Theorem. We demonstrate that inertia is a key parameter, damping the system's response and slowing its relaxation to a periodic non-equilibrium steady state. Our results show that oscillatory thermal driving is a powerful tool for controlling nanoscale energy flow.

Figures

Equations

$$\dot{x}(t)=v(t),m\dot{v}(t)=-\gamma v(t)-{\partial }_{x}V(x(t))+\sqrt{2\gamma k_{B}T(t)}\xi (t),$$

$$⟨x(t)⟩$$

$${\sigma }_x^2(t)$$

$${\sigma }_v^2(t)$$

$$⟨x(t)v(t)⟩$$

$$T(t)=T_0+T\sin (\Omega t),$$

$${\rho }_{xv}(t)=\frac{⟨x(t)v(t)⟩-⟨x(t)⟩⟨v(t)⟩}{{\sigma }_v(t){\sigma }_x(t)}.$$

$$Q[\vec{x}(t)]={\int }_{x_0}^x\left(-\gamma v(t)+\xi (t)\right)\circ dx(t)=\frac{1}{2}m[v^2(t)-v_0^2]+\frac{1}{2}k[x^2(t)-x_0^2].$$

$$P(Q)=⟨\delta (Q-Q[\vec{x}(t)])⟩=\frac{1}{2\pi }\int d\lambda e^{i\lambda Q}⟨e^{-i\lambda Q[\vec{x}(t)]}⟩.$$

$$Z(\lambda )=⟨e^{-i\lambda Q[\vec{x}(t)]}⟩=\int dx\int dv\int d{x}_0\int d{v}_{0}P(x,v,x_0,v_0)e^{-i\lambda Q[\vec{x}(t)]}.$$

$$Z(\lambda )=\frac{1}{\sqrt{\det (\mathbf{\Psi }+{\mathbf{I}}_4)}}$$

$$\displaystyle\mathbf{\Psi}=i\lambda\left[\begin{array}[]{cccc}k\,\sigma_{x}^{2}(t)&k\,\langle x(t)x_{0}\rangle&k\,\langle x(t)v(t)\rangle&k\,\langle x(t)v_{0}\rangle\\ -k\,\langle x(t)x_{0}\rangle&-k\,\sigma_{x_{0}}^{2}&-k\,\langle v(t)x_{0}\rangle&0\\ m\,\langle x(t)v(t)\rangle&m\,\langle v(t)x_{0}\rangle&m\,\sigma_{v}^{2}(t)&m\,\langle v(t)v_{0}\rangle\\ -m\,\langle x(t)v_{0}\rangle&0&-m\,\langle v(t)v_{0}\rangle&-m\,\sigma_{v_{0}}^{2}\end{array}\right],$$

$$⟨x_0{v}_0⟩$$

$${\sigma }_{x_0}^2$$

$${\sigma }_{v_0}^2$$

$$⟨x(t)x_0⟩$$

$$⟨x(t)v_0⟩$$

$$⟨v(t)x_0⟩$$

$$⟨v(t)v_0⟩$$

$$\displaystyle P(Q)\propto\left\{\begin{array}[]{rcl}e^{-\beta_{+}Q}&\mbox{for}&Q>0\\ e^{-\beta_{-}|Q|}&\mbox{for}&Q<0.\end{array}\right.$$

$$\frac{P(Q)}{P(-Q)}=e^{\Delta \beta Q},$$

$$⟨Q(t)⟩$$

$$\displaystyle\langle Q^{n}\rangle=i^{n}\frac{\partial^{n}Z(\lambda)}{\partial\lambda^{n}}\Big{|}_{\lambda=0}.$$

$$\ddot{x}(t)+\beta \dot{x}(t)+{\omega }_0^{2}x(t)=A(t).$$

$$x(t)=A{e}^{{\alpha }_{1}t}+B{e}^{{\alpha }_{2}t},$$

$$\frac{dA}{dt}{e}^{{\alpha }_{1}t}+\frac{dB}{dt}{e}^{{\alpha }_{2}t}=0,$$

$$x(t)$$

$$v(t)$$

$$⟨x^2(t)⟩$$

$$⟨v^2(t)⟩$$

$$⟨x(t)v(t)⟩$$

$${\sigma }_x^2(t)$$

$${\sigma }_v^2(t)$$

$$⟨x(t)v(t)⟩$$

$${\partial }_{t}P(x,v,t)=\left[-v{\partial }_x+{\partial }_v\left(\gamma v+\frac{{\partial }_{x}V(x)}{m}\right)+\frac{\gamma k_{B}T(t)}{m}{\partial }_v^2\right]P(x,v,t).$$

$$P(x,v,t)$$

$$P(x_0,v_0,0)=\frac{1}{2\pi {\sigma }_{x_0}{\sigma }_{v_0}}\exp \left(-\frac{x_0^2}{{\sigma }_{x_0}}-\frac{v_0^2}{{\sigma }_{v_0}}\right),$$

Full text

Keywords: Brownian motion, underdamped, thermal protocol, heat distribution

1 Introduction

In mesoscopic systems, energy exchanges are inherently stochastic. Thermodynamic quantities, such as work and heat, exhibit fluctuations and are therefore described within the framework of stochastic thermodynamics [1, 2, 3]. In this framework, variable temperature processes emerge and play a fundamental role in many contexts, e.g., Brownian engines, micro-particle cooling, and biological systems [4, 5, 6].

To describe the motion of particles subjected to random forces in colloidal environments, a theoretical framework is provided by the Langevin and Fokker–Planck equations [7, 3]. These equations are employed in two distinct regimes, namely the underdamped and overdamped limits, each capturing different features of the stochastic behavior. The underdamped regime ($m\neq 0$) accounts for inertia and preserves correlations between the particle’s position and velocity in non-equilibrium situations, whereas the overdamped regime ($m=0$) neglects inertia, decouples position from velocity, and is appropriate for systems with fast relaxation times. While both regimes have been extensively studied from a dynamical perspective [8, 9, 10, 11, 12, 13], their thermodynamic implications — particularly the role of inertia in energy exchange, dissipation, and system performance — remain less understood [14, 15, 16, 17, 18, 19].

In stochastic thermodynamics, heat exchange between a Brownian particle and its thermal environment is always present, whereas work is not necessarily involved [20]. Work is the ordered energy exchange between a system and an external agent in a controlled or directed manner. Therefore, for work to be performed, the Brownian particle must interact with an external agent, an external system [3]. In contrast, heat is the spontaneous and disordered energy exchange with the thermal bath, and occurs inherently in such systems [3]. Despite their fundamental role, the energetics of Brownian systems cannot be directly accessed — typical experiments track the particle’s trajectory over time [21, 22]. To extract thermodynamic quantities, one must infer them from the statistical properties of the trajectories [23]. Consequently, the theoretical framework is essential for interpreting empirical data and characterizing the underlying thermodynamic processes.

Focusing on optically trapped Brownian particles, new experimental methods are already being implemented to mimic effective, and thereby controllable, thermal baths [24, 25, 26, 27, 28, 29]. Such controlled manipulation opens avenues for exploring a broad class of Brownian engines. Consequently, characterizing the thermodynamic properties of these non-isothermal processes is of considerable relevance for both fundamental studies in non-equilibrium systems and applications to nanoscale engines.

Motivated by these considerations, in this work we investigate how a temperature-time-dependent protocol affects the dynamics and thermodynamics of an underdamped Brownian particle in a harmonic potential. Specifically, we calculate analytical expressions for the system’s behavior under a generic thermal protocol and then evaluate these results for the case of a sinusoidal thermal protocol. Combining analytical methods with numerical simulations, we compute both the Pearson correlation coefficient between position and velocity and the heat distribution. For the Pearson correlation, we found a non-vanishing periodic correlation with resonant phenomena, while the heat distribution is highly sensitive to parameter control and is in accordance with the Fluctuation Theorem.

This paper is organized as follows: In Section 2, we describe the dynamics of the underdamped Brownian particle subjected to periodic thermal driving in a harmonic potential. In Section 3 we analyze the dynamics of the Pearson correlation coefficient. In Section 4 we analyze the probability distribution of the heat, derive the characteristic function, and calculate the moments of the distribution. We finish in Section 5 with the discussion of the results and conclusions.

2 Brownian dynamics with arbitrary temperature protocol

We consider an underdamped Brownian particle model, with a temperature-time-dependent protocol $T(t)$. The Langevin equations are

[TABLE]

where $m,\gamma,k_{B}$ and $V(x(t))$ are the mass, damping coefficient, Boltzmann constant, and potential energy, respectively. $\xi(t)$ is Gaussian white noise with a zero mean and delta correlated, i.e., $\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. Throughout this paper, we consider a harmonic potential $V(x)=kx^{2}/2$, where $k$ is an effective spring constant — for a Brownian particle trapped by an optical tweezer, $k$ is the stiffness of the trap and depends linearly on the optical power [30]. This potential implies that Eq.(1) is linear, and since the statistics of $x(t)$ and $v(t)$ depend linearly on the Gaussian noise, the statistics of $x(t)$ and $v(t)$ are also Gaussian; therefore, the quantities of greater interest — which dictate the probability distribution — are the first and second moments. For a generic temperature-time-dependent protocol, these moments are given by

[TABLE]

where $\beta=\gamma/m$, $\omega_{0}^{2}=k/m$ and $\omega_{1}=\sqrt{\omega_{0}^{2}-\beta^{2}/4}$ (See A for details). The foregoing equations assume that the system is initially at equilibrium with temperature $T_{0}$, i.e., the initial conditions obey $\langle x_{0}\rangle=\langle v_{0}\rangle=0$, $\langle x^{2}_{0}\rangle=k_{B}T_{0}/k$ and $\langle v^{2}_{0}\rangle=k_{B}T_{0}/m$. In addition, it is assumed $\omega_{1}>0$; if it is imaginary, the equations must be rewritten in terms of hyperbolic functions, using the relations $\cos(ix)=\cosh(x)$ and $\sin(ix)=i\sin(x)$.

Despite the calculations presented here being applicable to a generic temperature-time-dependent protocol, we focus the analysis on a specific case, a sinusoidal thermal protocol

[TABLE]

where $T$ defines the amplitude of the oscillations and $\Omega$ is the driving frequency. Such protocols have recently been suggested as a possible model for mechanical excitation and heat conduction in hydrocarbons [31].

3 Dynamics

We analyze the Pearson correlation coefficient between displacement and velocity

[TABLE]

An important property of the Pearson correlation coefficient is $|\rho_{xv}(t)|\leq 1$. Moreover, if its absolute value is equal to one, $|\rho_{xv}(t)|=1$, then there exist constants a and b such that v = ax + b with unit probability, i.e. v(t) depends linearly on x(t) [32].

For a system with constant temperature, i.e., in thermal equilibrium, it is well known that $\rho_{xv}(t)$ is null. However, with a sinusoidal temperature protocol, i.e., modulation of noise, a residual periodic correlation appears, as can be seen in Fig. 1.

Due to the time-dependent temperature $T(t)$, the naively overdamped approximation ($m/\gamma\rightarrow 0$) fails to describe the long-time behavior ($t\gg m/\gamma$). The oscillatory temperature protocol imposes a rapid, periodic driving at frequency $\Omega$ that prevents the velocity from fully relaxing before $T(t)$ changes again. Physically, the small but non-zero value of $\rho_{xv}(t)$ is a signature of the residual inertial effects in a driven steady state. As the particle heats and cools (i.e., $T(t)$ oscillates) energy is periodically injected into the system, sustaining inertial effects.

As displayed in Fig. 2 the set of mass $m$, effective spring constant $k$, and friction coefficient $\gamma$ has a substantial influence on the correlation $\rho_{xv}(t)$.

Numerical results in Fig. 2(a) reveal the existence of two distinct regimes, set by a critical mass $m_{c}=\gamma^{2}/4k$. For $m>m_{c}$, $\omega_{1}$ is real, inertia dominates over friction, and the dynamics exhibits an oscillatory relaxation. For $m<m_{c}$, $\omega_{1}$ is purely imaginary, friction dominates over inertia, and the dynamics relaxes exponentially. At $m=m_{c}$ the system reaches the critical damping point, marking the boundary between these two regimes. At this precise mass, the system’s response to any perturbation changes its fundamental character. This transition marks a fundamental change in the system’s dynamical behavior. Moreover, in both regimes a peak is observed, indicating an optimal response for a given mass $m$ — a clear manifestation of a resonant phenomenon. On the other hand, as $m\rightarrow\infty$, the position-velocity correlation is significantly reduced and tends to zero.

One might be induced to think that the drive frequency $\Omega$ only influences the oscillation frequency on the Pearson correlation $\rho_{xv}(t)$. However, as displayed in Fig. 2(b), the amplitude of the correlation oscillations (i.e., the contrast between the dark and bright bands) varies with $\Omega$. The strongest correlations occur for driving frequencies in the range $\Omega\in[1,2]$, indicating the presence of another resonance phenomenon.

4 Heat fluctuations

Since the parameters of the system do not change, there is no work done upon it. So, from the first law of stochastic thermodynamics, the stochastic heat — which is a functional of the trajectory — is equal to the variation in the internal energy $Q[\,\vec{x}(t)\,]=\Delta K+\Delta V$. For this reason, the stochastic thermodynamics of the system can be characterized solely by the heat, and therefore we focus on it.

The stochastic heat is the competition between the noise of the bath and the drift force [3]. By the underdamped dynamics in Eq.(1), we can construct the stochastic heat functional

[TABLE]

Stratonovich-type product $\circ$ is adopted in the foregoing equation — as in macroscopic thermodynamics, $dQ>0$ assigns heat received by the Brownian particle, and $dQ<0$ assigns heat released to the thermal bath. Despite having a time dependency in the temperature, the heat functional is just the sum of the difference in the kinetic energy and the potential energy.

4.1 Heat distribution

The random values of the heat $Q$ are given by the trajectory-dependent term $Q[\vec{x}(t)]$. Thereby, the probability density function (PDF) for the heat is obtained by averaging over the trajectories

[TABLE]

The foregoing equation can be rewritten in terms of the characteristic function $Z(\lambda)$  [33]

[TABLE]

Since the integrals in Eq.(10) are convergent for any generic smooth and continuous temperature-time-dependent protocol, $Z(\lambda)$ contains the same amount of information as $P(Q)$. On the other hand, the integral in Eq.(9) cannot be written in terms of conventional functions, i.e., $P(Q)$ cannot be calculated analytically. For this reason, the heat PDF is numerically obtained, and the characteristic function is calculated analytically, which yields

[TABLE]

where $\mathbf{I}_{4}$ is the identity matrix and $\mathbf{\Psi}$ is given by

[TABLE]

where its elements are

[TABLE]

(See B for details). As one can see in Fig. 3, our numerical results present strong agreement with what is expected from the theory.

Remarkably, we notice that the heat distribution is asymmetric around its peak. This asymmetry violates detailed balance, suggesting that the system is far from equilibrium. Moreover, Fig. 3 implies that the heat distribution curves are essentially linear in the semi-log scale, i.e., $P(Q)\propto e^{\beta_{i}Q}$. Since the slopes for $Q>0$ and $Q<0$ are different, there are two different decay constants, such that

[TABLE]

This asymmetric exponential behavior is a direct manifestation of the heat fluctuation relation for a Brownian particle and a nonstationary bath [34, 35, 36]

[TABLE]

where $\Delta\beta$ is related to the degree of entropy production of the system. Thereby, the asymmetry observed is a quantitative feature predicted by non-equilibrium statistical mechanics.

4.2 Statistical moments

The mean heat can be easily obtained by averaging Eq.(8) over the trajectories, which yields

[TABLE]

As one can see in Fig. 4, our numerical results present a nice agreement with what is expected from the theory.

The system is continuously exchanging heat, following the thermal modulation $T(t)$. Comparing the green ($\Omega=1$) and orange ($\Omega=5$) curves in both graphs, we can see the effect of the driving frequency. For instance, for both the light (Fig. 4(a)) and heavy (Fig. 4(b)) particles, the slower driving frequency results in a much larger amplitude of heat exchange. This happens because when the temperature changes slowly ($\Omega=1$), the system has more time to absorb or release heat before the cycle reverses. On the other hand, when the driving is fast ($\Omega=5$), the temperature changes so rapidly that the particle can’t exchange much heat before it reverses, leading to much smaller oscillations in $\langle Q(t)\rangle$. Furthermore, since the thermal protocol modulates noise strength and therefore the statistical energy exchange, it couples strongly into the particle’s energetic degrees of freedom. This produces a large mean heat response for $\Omega$ near $\omega_{0}$: the green $m=1$ curve exhibits the largest amplitude — a resonant-like amplification.

By comparing Fig. 4(a) ($m=1$) with Fig. 4(b) ($m=10$), the role of inertia becomes clear. The heavier particle has a significantly smaller amplitude of heat exchange compared to the lighter particle for the same driving frequency. This happens because inertia acts as a memory device and resists changes in motion: the heavier particle cannot respond as effectively to the changing thermal environment. Its kinetic energy lags behind the bath temperature, leading to less efficient and smaller heat exchange in each cycle.

The other statistical moments can be calculated by differentiating the characteristic function, Eq.(11), so that

[TABLE]

From the foregoing equation, we analyze the standard deviation, skewness, and excess kurtosis, as shown in Fig. 5.

Numerical results presented in Fig. 5(a) show that the light particle reaches the plateau much faster and has larger oscillations, whereas the heavy particle shows a slower rise and approaches more gradually. Fig. 5(b) confirms that the distribution remains asymmetric about its mean. In particular, the skewness $\mu_{3Q}$ changes sign periodically, and its positive peaks exceed the magnitude of its negative troughs. This bias toward positive skew reflects the fact that under the sinusoidal protocol, the particle absorbs heat more effectively than it releases. Moreover, an increase in the particle mass both reduces the amplitude of these skewness oscillations (damps the system’s response) and shifts the oscillations in time (creates a phase lag). Finally, the excess kurtosis displayed in Fig. 5(c) starts at $\kappa_{Q}=6$, indicating a strongly non-Gaussian, leptokurtic distribution and validating Eq.(27) at short times. As the system evolves, it relaxes into a periodic nonequilibrium steady state, set by the continuous driving. Furthermore, heavier particles, with greater inertia, exhibit longer relaxation times, so mass significantly prolongs the approach to that steady state.

5 Conclusion

In this work, we have studied a harmonically bound Brownian particle subjected to a sinusoidal thermal protocol. We employ the Langevin formalism to describe the dynamics and the stochastic energetics, which in this setup reduce to heat and internal energy. We analyze the Pearson correlation coefficient between displacement and velocity, and we derive analytical expressions for the characteristic function of the heat distribution from which its statistical moments are obtained. The heat distribution is evaluated numerically, and all theoretical predictions agree well with numerical simulations (see C for simulation details).

We find that the system dynamics and heat exchange depend sensitively on several parameters, including the driving frequency, the amplitude of the temperature oscillation, and the particle mass. The thermal frequency induces a time-scale competition between driving and relaxation: higher $\Omega$ would strengthen asymmetry (less time to equilibrate), while lower $\Omega$ would weaken it. This effect manifests in high order moments: both skewness and excess kurtosis are non-zero for $t>0$, confirming that the heat distribution $P(Q)$ is fundamentally non-Gaussian. Moreover, the heat flux is time-dependent, as expected, due to the thermal noise modulation. In the long-time limit, the driving periodicity is inherited by the probability distribution and by ensemble averages [37, 38, 39]; consequently, the system relaxes to a time-periodic non-equilibrium steady state.

Inertial effects — and thus the Pearson correlation coefficient — remain important when $T(t)$ changes on a time scale comparable to or faster than the velocity relaxation time. The particle’s mass is a critical parameter; it dampens the system’s response (lower amplitude of skewness oscillations), delays the response (phase lag in skewness) and slows down the relaxation to a steady state dramatically (longer transient time for standard deviation and excess kurtosis).

In conclusion, the results presented here demonstrate that dynamical and inertial contributions to the heat, induced by rapid thermal oscillations, produce complex hysteresis in energy flux. Furthermore, our findings indicate that time-periodic temperature modulations offer a viable route to control energy flow at the nanoscale. Exploiting these effects to design thermal devices constitutes a promising direction for future research.

Acknowledgments

FPA would like to acknowledge Dr. Leandro G. Rizzi for insightful discussion. The authors are especially grateful to Dr. Pedro V. Paraguassú, whose feedback and extensive discussions were fundamental to this study. Although not listed as co-author, his contributions were essential to the development of this manuscript.

This work is supported by the Brazilian agencies Coordenação de Aperfeiçoamento de Pessoal de Ensino Superior (CAPES) and Fundação Carlos Chagas Filho de Apoio à Pesquisa do Estado do Rio de Janeiro (FAPERJ). FPA acknowledges CAPES scholarship, finance code 001 and FAPERJ (Grants No. E-26/201.820/2025). WAMM acknowledges CNPq (Grant No. 308560/2022-1).

Appendix A Moments derivation

For a Brownian particle subject to a harmonic potential, we define the friction rate $\beta=\gamma/m$, natural frequency $\omega_{0}^{2}=k/m$, and assume that the stochastic force satisfies the fluctuation–dissipation relation $\langle A(t)\rangle=0$ $\langle A(t)A(t^{\prime})\rangle=2k_{\mathrm{B}}T(t)\beta\delta(t-t^{\prime})/m$. Under these definitions, the Langevin equation (Eq.(1) in the main text) reads

[TABLE]

From this equation, we seek to obtain the probability distribution $P(x,v,t;x_{0},v_{0},t_{0})$. To do so, we first obtain the formal solution of Eq.(30) regarded as an ordinary nonhomogeneous differential equation. The solution can be obtained in terms of the solution of the homogeneous equation given by

[TABLE]

where $\alpha_{1,2}=-\beta/2\pm i\omega_{1}$ and $\omega_{1}=\sqrt{\omega_{0}^{2}-\beta^{2}/4}$. To obtain the full inhomogeneous solution, we let the constants $A,B$ become time-dependent functions and impose the auxiliary condition [40]

[TABLE]

as well as the initial conditions $x(0)=x_{0}$ and $v(0)=v_{0}$. Therefore, one finds the solution

[TABLE]

Imposing an initial equilibrium distribution at temperature $T_{0}$, i.e., $\langle x_{0}\rangle=\langle v_{0}\rangle=0$, $\langle x_{0}^{2}\rangle=k_{B}T_{0}/k$, and $\langle v_{0}^{2}\rangle=k_{B}T_{0}/m$ and using the fluctuation–dissipation relation, one shows from equations Eq.(A) and Eq.(A) that $\langle x(t)\rangle=\langle v(t)\rangle=0$, and the second moments and cross-correlation, respectively, reads

[TABLE]

Initially, there was a Heaviside step function $\theta(t-s)$ in the integrand of the equations above, but the term can be suppressed since the integration limit is from [math] to $t$ and in this region $\theta(t-s)=1$. Since $\langle x(t)\rangle=\langle v(t)\rangle=0$, the variances are $\sigma_{x}^{2}(t)=\langle x^{2}(t)\rangle$ and $\sigma_{v}^{2}(t)=\langle v^{2}(t)\rangle$. After some further manipulation on the foregoing expressions, equations Eq.(2), Eq.(2) and Eq.(2) in the main text are obtained.

A.1 Moments for the sinusoidal thermal protocol

Considering the sinusoidal temperature protocol $T(t)=T_{0}+T\sin(\Omega t)$ and some further manipulation, equations Eq.(2), Eq.(2) and Eq.(2) in the main text reads

[TABLE]

Appendix B Characteristic function

The evolution of the probability distribution $P(x,v,t)$ from $t=0$ to $\tau$ is described by the Klein-Kramers equation [41]

[TABLE]

The Green’s function of Eq.(41) can be written as the bivariate Gaussian distribution

[TABLE]

If the initial probability distribution is Gaussian and is uncorrelated in position and velocity, i.e.,

[TABLE]

then the distribution remains Gaussian for all later times and is given by

[TABLE]

where $\mathbf{z}=[x,x_{0},v,v_{0}]^{T}$ is a $4$-dimensional column vector, and $\mathbf{\Sigma}$ is the covariance matrix

[TABLE]

Its elements are obtained using Eq.(A) and Eq.(A) and the initial equilibrium distribution, such that

[TABLE]

The characteristic function can be obtained by applying Eq.(44) and Eq.(8) in Eq.(10) and performing the Gaussian integrals, such that

[TABLE]

where $\mathbf{M}=\mathbf{Q}+\mathbf{\Sigma}^{-1}$ and

[TABLE]

Finally, after further simplification, the characteristic function becomes

[TABLE]

where $\mathbf{I}_{4}$ is the identity matrix and

[TABLE]

The probability normalization condition imposes that $Z(0)=1$.

Appendix C Numerical simulations

Considering a finite time interval $\Delta t=t_{i+1}-t_{i}$, Eq.(1) can be discretized and rewritten as

[TABLE]

$\tilde{F}_{a}(t)$ is a random force evaluated as [42]

[TABLE]

where $N(0,1)$ denotes a Gaussian variable with zero mean and variance equal to one. Thus, the stochastic heat can be obtained in terms of the trajectory, Eq.(69), such that

[TABLE]

This method is a physical approach provided by a simple discrete-state jump Markov process; it allows modeling the immediate effects of individual molecular collisions on the velocity of the Brownian particle. Despite this being a discrete and idealized model, the continuum limit duplicates the predictions of the Langevin equation [43].

For each simulation we consider $N_{t}=10^{6}$ independent trajectories, a time interval $\Delta t=10^{-4}$, effective spring constant $k=1$, temperatures $T_{0}=20$ and $T=5$, and the Boltzmann constant $k_{B}=1$. Also, we set initial conditions $x_{0}$ and $v_{0}$ in thermal equilibrium at temperature $T_{0}$. So, for each trajectory

[TABLE]