Quantum Thermodynamic Uncertainty Relations, Generalized Current Fluctuations and Nonequilibrium Fluctuation-Dissipation Inequalities

Thermodynamic uncertainty relations (TURs) represent one of the few broad-based and fundamental relations in our toolbox for tackling the thermodynamics of nonequilibrium systems. One form of TUR quantifies the minimal energetic cost of achieving a certain precision in determining a nonequilibrium current. In this initial stage of our research program, our goal is to provide the quantum theoretical basis of TURs using microphysics models of linear open quantum systems where it is possible to obtain exact solutions. In paper [Dong \textit{et al.}, Entropy {\bf 24}, 870 (2022)], we show how TURs are rooted in the quantum uncertainty principles and the fluctuation-dissipation inequalities (FDI) under fully nonequilibrium conditions. In this paper, we shift our attention from the quantum basis to the thermal manifests. Using a microscopic model for the bath's spectral density in quantum Brownian motion studies, we formulate a ``thermal'' FDI in the quantum nonequilibrium dynamics which is valid at high temperatures. This brings the quantum TURs we derive here to the classical domain and can thus be compared with some popular forms of TURs. In the thermal-energy-dominated regimes, our FDIs provide better estimates on the uncertainty of thermodynamic quantities. Our treatment includes full back-action from the environment onto the system. As a concrete example of the generalized current, we examine the energy flux or power entering the Brownian particle and find an exact expression of the corresponding current-current correlations. In so doing, we show that the statistical properties of the bath and the causality of the system+bath interaction both enter into the TURs obeyed by the thermodynamic quantities.


I. INTRODUCTION
Perhaps one of the most interesting recent developments in the pursuit to understand nonequilibrium sciences are the family of thermodynamic uncertainty relations [1][2][3][4][5][6][7][8] (see also Refs. [9][10][11]). In one form, a TUR relates the fluctuations of a nonequilibrium current with the minimum of dissipation energy during the nonequilibrium process. Originally derived for classical, Markovian systems, the TUR has been quickly extended to over- [12] and underdamped [13] dynamics of a Brownian particle as well as Markovian quantum systems by means of large deviation methods [14] and the Cramér-Rao bound for the quantum Fisher information [6]. The latter technique has also been used to derive a connection between the TUR and fluctuation theorems for classical systems [15] and has just recently been used to derive an extension of the TUR for quite general open quantum systems [7]. We further mention that a TUR has been investigated for steady-state heat transfer [5,16] and a similar expression could be found for the relation between the entropy production and the time it takes to complete a nonequilibrium process [17].
Aiming at finding the quantum roots of the uncertainties of thermodynamic quantities in nonequilibrium sys-tems, we demonstrated [18] that for Gaussian open quantum systems, thermodynamic functions are functionals of the Robertson-Schrödinger uncertainty (RSU) function. Using recent results on the nonequilibrium free energy and nonequilibrium effective temperature [19], we showed that a fluctuation-dissipation inequality (FDI) exists at all times in the nonequilibrium dynamics of the open system. In this sequel paper, we continue these veins of investigation and show how a thermodynamic uncertainty relation (TUR) for macroscopic quantities can formally be derived. While these are often motivated by phenomenological considerations in the literature, we show how some such relations can be obtained in rigorous ways within a microscopic quantum, even quantum field theory, framework. The centerpiece is the generalized current fluctuations. In short, from the mathematical perspective, while our first paper [18] is concerned with physical quantities that are derived from secondorder correlations functions, we now direct our attention to higher-order correlation functions and their physical imprints on Gaussian open quantum systems. Let us examine the ingredients one by one, and then look at their synthesis which leads to the aforementioned inequalities and uncertainty relations.
themselves as a useful tool in controlling physical systems on the micro-and nanoscale. A precise understanding and manipulating of the fluctuations is hence crucial in the design of future technologies, especially in the strive for further miniaturization of devices (see, e.g., the recent reviews [20][21][22][23] and references therein). Some fluctuations are introduced by operational techniques or imperfections and can in principle be eradicated by optimizing the experimental protocol. Other fluctuations are more fundamental and can never be circumvented. Possibly the most famous of such fluctuations are the quantum fluctuations whose physical implications are summarized by the quantum uncertainty principles (QUP) [24] and, in particular, the Robertson-Schrödinger uncertainty principle [25,26]. In many realistic scenarios, the system of interest is also coupled to an environment. Thermal noise enters from a finite temperature bath, and a finite coupling strength between the open system and its environment can also be represented by noises. The nonequilibrium dynamics of the open quantum system can further modify the uncertainty relations [27][28][29].
In equilibrium, the corresponding variance of a quantum observable in an open system is accounted for by the fluctuation-dissipation theorem [30,31]. It states an exact relation between variance and dissipation due to a detailed balance between the average incoming and outgoing power [32].

Nonequilibrium Relations
Departing from equilibrium enriches the fluctuation spectrum with respect to the corresponding equilibrium situation. The search for the statistical description of both classical and quantum systems in nonequilibrium situations has instigated the development of a number of theorems. To name a few examples, we find the fluctuation theorems [33][34][35] (see, e.g., also the reviews of refs. [36][37][38][39]), the nonequilibrium fluctuation-dissipation relations for steady states [40][41][42][43][44] (see also refs. [45][46][47][48] for nonequilibrium fluctuation-dissipation relations of classical active systems), and the fluctuation-dissipation inequality [49,50]. The latter is a reflection of the nonequilibrium condition the system exists in for the full duration before it reaches equilibrium. In simple terms, it states that the energy stored in the fluctuations of the environment (statistical operator) is always equal to or exceeds the energy dissipated into the environment [49,51]. It can hence be understood as a generalization of fluctuation-dissipation relations which can often be shown to hold at late times after a system has settled down to equilibrium while interacting with its environment.

Fluctuations of Generalized Current
Due to their fundamental nature, it can be expected that the (often microscopically formulated) uncertainty principles are to some extent inherited by certain macroscopic observables. This line of thought draws attention to another perspective on uncertainties in nonequilibrium systems that has seen a surge of interest in the recent years, namely the thermodynamic uncertainty relation (TUR) [9]. Despite versatile contexts and examples (see beginning of current section), the particular relations can often be summarized in terms of a generalized currentĴ. May it be, for example, the position of a particle [1], certain measurements on atomic systems [6], or the exchanged energy in heat transfer [5], the thermodynamic uncertainty relation states that the fluctuations of the generalized current Ĵ 2 are always larger than or equal to the average current Ĵ , with a thermal proportionality factor depending on the temperature of the heat bath(s). If the average current can be connected to dissipation from the system of interest into the environment, the TUR provides a nonequilibrium measure for the thermodynamic cost of achieving a desired precision on measurements of the current [1].

B. This Work-Key Findings and Organization
In the present manuscript, we aim to connect the microscopic picture of interacting quantum systems (quantum uncertainty principles) with the macroscopic picture of nonequilibrium thermodynamics (thermodynamic uncertainty relation) from a conceptual point of view. To this end, based on the generic model for quantum Brownian motion, we first study the connection between fluctuations in the system and the bath spectral density. In this way, we extend previous work on the nonequilibrium fluctuation-dissipation inequality [18] to something we call "thermal fluctuation-dissipation inequality" which additionally takes the impact of the bath temperature into account. Secondly, we calculate the fluctuations of the current of energy entering the system. We then combine the exact result for the current fluctuations with the thermal fluctuation-dissipation inequality to derive a combined inequality. This inequality incidentally turns out to formally represent a thermodynamic uncertainty relation, i.e., it is proportional to the average current times a thermal factor. In contrast to parts of the related literature, our result -showing the emergence of a thermodynamic uncertainty relation from the nonequilibrium evolution -fully incorporates the non-Markovian features of the system+bath interaction (exact in all orders of coupling). Since our formalism can be easily extended to the problem of heat transfer, we confirm the soundness of our results by reproducing a known TUR for steady-state heat transfer [5,52]. It turns out that a limited number of assumptions is needed to connect the two worlds -quantum uncertainty principles and the TUR -which, in turn, establishes a clear hierarchy of different inequalities that are connected to their respective physical perspectives: from microscopic quantum mechanics, over thermodynamic quantities, and all the way to concrete physical realizations. This paper is structured as follows. We work in an open quantum system framework [53,54] and introduce the quantum stochastic dynamics of the system using the time-honored quantum Langevin equation of generic quantum Brownian motion in Section II A [55][56][57][58][59][60][61][62][63][64]. In Section II B, we deduce the fluctuation-dissipation inequality and derive the space-momentum uncertainty. We relate the two to the power flow between system and environment as well as the total entropy production and the entropic uncertainty relation during the equilibration process in Section III. We then specify the spectral density and the dissipation kernel characterizing the environment of the particle (see Section IV A) and derive a thermodynamic uncertainty relation (see Section IV). Specifying the fluctuation-dissipation inequality to finite temperatures to include thermal noise contributions (see Section IV A), we can highlight when thermodynamics enters the stage by building a concrete connection to the functional form of the thermodynamic uncertainty relation in Section V. In Section VI, we briefly show how our formalism can be applied to heat transfer and comment on the consistency with the respective (steady-state) TURs that have been derived earlier. We conclude our manuscript with a discussion in Section VII.

II. STOCHASTIC DYNAMICS OF GAUSSIAN SYSTEMS: 2ND AND HIGHER-ORDER CORRELATIONS
We begin with a brief review of stochastic dynamics of open quantum systems using the ubiquitous quantum Brownian motion model. We show from the Langevin equation how the open system's dissipative dynamics is linked to the fluctuations in its environment, registered in the dissipation (response) and noise (correlation) kernels. From this, we show how to obtain the second and fourthorder correlations of currents. This prepares us to tackle the main tasks we set forth in our goals.

A. Fluctuations and Stochastic Dynamics of Open Quantum Systems
In order to introduce our system, we loosely follow refs. [49,63,[65][66][67][68] and refer readers to the monograph [69] for a more comprehensive discussions.
We consider the dynamics of a single bosonic quantum degree of freedomq under the influence of a quadratic potential V (q) = ω 2 0q 2 /2 with bare frequency ω 0 [70]. The system is coupled to a bosonic bath at finite temperature T described by the generic environmental operatorÊ . We assume a simple bilinear coupling [68] where we absorbed the coupling constant in the definition ofÊ. We note that this choice of coupling implicitly demands that the Hamiltonian is appropriately renormalized in order to account for the (Lamb-)shift of the potential's minimum frequency ω 0 due to the interaction with the environment [65,68] [see also Equation (7)]. Under the assumption of a linear (Gaussian) environment, the bath operator can be expressed self-consistently as [49,63,65].Ê whereξ is the unperturbed stochastic operator of the environment in absence of the system, and we set the initial time of the experiment to t i = 0. The integral in Equation (2) connects the impact of the system on the environment which, in turn, back-acts on the particle by means of the real-valued and causal response kernel Γ(t, τ ). The noise spectrum is intimately connected to the response kernel.
where we average over the density matrices of the system ρ S and the bathρ E , which we assume to factorize at initial time, i.e.,ρ(0) =ρ S (0) ⊗ρ E (0). The real-valued correlation function ν(t, t ) and the kernel Γ(t, t ) are given by [68].
which, for a given set of initial operators {q(0),q(0)}, fully determines the dynamics of the system. For simplicity, we from now on consider stationary functions ν(t, t ) = ν(τ = t − t ) and Γ(t, t ) = Γ(τ ). Equation (5) allows for an intuitive interpretation of the separate terms: due to the resemblance between Equation (5) and the equation of motion for a randomly moving particle in a thermal viscous environment, it is usually referred to as quantum Brownian motion [57][58][59].ξ acts as noise "driving" the system. The kernel Γ, on the other hand, encodes the dissipative processes, i.e., energy losses of the system to the environment, as well as the rescaling of the bare frequency ω 0 . The latter can be made explicit by defining the dissipation kernel.
−∂ τ γ(τ ) = Γ(τ ) (6) which yields after a partial integration where we have redefined the resonance energy asω 2 0 = ω 2 0 − γ(0). In order to keep notation simple, we do not print the extra tilde in the following. Furthermore, in the remainder of the manuscript, we assume that q(0) = q (0) = 0. We note that, at the price of restricting to linear systems, the previous equation of motion is non-Markovian and general to all orders of coupling between system and environment, i.e., we are not limited by the Born-Markov or rotating wave approximation [71][72][73]. It is well-known that the quantum Langevin equation [Equation (7)] can be solved for the canonical pair Q(t) := [q(t)q(t)] T by means of the response function (see Appendix A for details).
as the linear susceptibility and γ θ (t) = θ(t)γ(t) and θ(t) as the Heaviside step function. We note that we distinguish between a function and its Fourier transform only by the different argument. The corresponding second-order fluctuations are given by with σ 0 = ∆Q 2 (0) s as the covariance matrix at t = 0, ν(τ ) = diag[0, ν(τ )], X(t) as a matrix comprised of (timederivates of) the response function (see Equation (A2)), and the superscript T indicates the transpose of a matrix. It is interesting to note that, at arbitrary times, the autocorrelation ofq is not a stationary function, i.e., even though the response kernel χ and the correlator of the environment degrees of freedom ν are. Indeed, Equation (10) is given by a complex convolution of the system's self-consistent interaction with the environment. On top of that, the initial conditions evolve by construction, not necessarily in a stationary way. True stationarity can only be achieved at late times, where the impact of the initial conditions on the dynamics abates and the convolution can be expanded in Fourier modes of the form exp(−iω[t − t ]) (see, e.g., Section III.B of ref. [74] for details). Since the noise operator assumes a general Gaussian form, the statistics of the system are fully determined by the two-point function [65]. We hence have that any even-order correlation reduces to the sum of twopoint functions of all possible pairings of operators with preserved order, while any odd-order correlation vanishes [75]. For example, we obtain for the fourth-order correlation function (see Appendix B for a detailed proof).
We note that the previous relation is not time-ordered and will come in handy when we compute the current fluctuations (see, e.g., Section V).

B. Fluctuation-Dissipation Inequality and Robertson-Schrödinger Relation
Due to its strong formal resemblance to classical equations of motion, the quantum Langevin equation perhaps evokes the impression that the dynamics of the particle progresses deterministically: the particle absorbs energy from its environment via the "force"ξ, processes it following its harmonic constraints, and emits it back into the environment; quantitatively described by the dissipation kernel γ. However, such a descriptions lets slide the quantum stochastic properties of the system. During the nonequilibrium evolution, even though the self-consistency of our approach ensures thermodynamic stability at all times, we are in lieu of an exact and transparent relation between the fluctuations of the environment [ν(τ )], the fluctuations of the reduced system of interest [ q(t)q(t ) s ], and the corresponding dissipation [γ(τ )]. Instead, even for systems which do not possess a fluctuation-dissipation relation, a corresponding fluctuation-dissipation inequality (FDI) [49,50] can be found, and can serve as useful bounds on the physical quantities of the system evolving under dynamical conditions. In a recent work [18], we showed that a fluctuation-dissipation inequality exists at all times in the nonequilibrium dynamics of open quantum systems. We traced back the uncertainties of thermodynamic quantities in nonequilibrium systems to their quantum origins. We summarize the main points of the FDI in Appendix A 1. For our purposes, it is sufficient to recall its main statement as [18,49].
for any complex function f . In particular, in a frequency domain, it is possible to show (see Appendix A 1) With this, it is straightforward to translate the results from the fluctuation-dissipation inequality for the environment to the fluctuations of the system of interest. Indeed, any property of ν and γ is inherited by Equations (8) and (9) such that also the fluctuations of the reduced systemq are a positive semidefinite function. However, there is no such simple frequency-domain relation [Equation (13)] as in the case for the environment, since q 2 is not stationary [Equation (10)] at arbitrary times t. Alternatively, we can compare the position-momentum uncertainty prescribed by the fluctuation-dissipation inequality to the Robertson-Schrödinger inequality [25,26]. In fact, it turns out that the former (FDI) can provide a stronger, more stringent bound to the uncertainties in the system at late times and only reduces to the Robertson-Schrödinger inequality by an additional application of the Cauchy-Schwarz inequality [18]. The same situation is achieved for vanishingly small system+bath coupling, as one would expect from the limiting case of conventional thermodynamics [74]. For details, we refer to Appendix A 2.

III. ENERGY FLOW AND ENTROPY PRODUCTION
We are now interested in the direct consequences of the fluctuation-dissipation inequality on the system's thermodynamic properties.
During the course of the interaction with the environment, the system's instantaneous mechanical Hamiltonian Ĥ (t) s = [q 2 (t) + ω 2 0q 2 (t)] s /(2) varies with time starting from initially ω 0 /2, eventually increasing to its equilibrium value at late times. From the perspective of the system (q), this is due to the (stochastic) interaction with the environment at any instance t and can be associated to the power entering P in and exiting P out the system, respectively. Upon multiplying from both sides of the quantum Langevin equation (Equation (5)) witḣ q, we obtain We remark that the incoming power is unaffected by the initial conditions. The outgoing power, on the other hand, is determined by the cross-correlation of position and momentum of the system, and hence by both the initial conditions as well as the response to the noisy environment (see Equation (A3) for the exact expression). Due to the linearity of our system+bath coupling and the fact that q(0)ξ s = q (0)ξ s = 0, the corresponding dynamics decouple, and we can write P out (t) = P init out (t) + P fluc out (t) (see also Equation (9)). As one could have expected, the information and energy initially stored in the system dissipates into the environment over time. Instead, over time, additional excitations are transferred to the system from the fluctuating environment. These environment-induced excitations again become dissipated. The system and environment are self-consistently back-acting onto one another. At late times, the system can reach a nonequilibrium steadystate (equilibrium in our case [74,76]). Concretely, upon defining the part of the system's dynamics, denoted bŷ Q(t), due to the noise operator ξ (Equation (A3)).
The change in the system's instantaneous mechanical energy due to the fluctuating interaction with the environment obeys the relation (see Appendix C for details).
At late times, we can further invoke the fluctuationdissipation inequality to derive a lower bound on the fluctuating part of the, say, outgoing power. Indeed, upon partially integrating with respect to τ in Equation (14b), we can write Performing the limit t → ∞, we obtain while P init out → 0. Since they equal at late times, the same relation holds for P in . Clearly, the net flow of energy into the system due to the interaction with the environment modifies its mean energy, and hence leads to a change in the particle's von Neumann entropy over the course of equilibration.
Due to the linearity of our system, the density matrix can be calculated exactly (see Appendix D and references therein). Quite remarkably, the result shows a formal resemblance to the density matrix of a quantum harmonic oscillator coupled to a bath with time-dependent temperature [19,52]. In particular, the entropy of the system at time t is given by [19].
One can readily check that this is equivalent to the von Neumann entropy S = S vN = −Tr[ρ logρ], where the trace is performed over the system's degrees of freedom, see, e.g., refs. [29,77,78]. From the perspective of the system, it is further interesting to note that the von Neumann entropy production rate is bounded from below by the time-derivative of its covariance matrix.
Starting with the contact between system and environment at t = 0 to the equilibration at t → ∞, the von Neumann entropy varies. The total production of entropy in the system is given by We note that this is not the total entropy production of the combined system, which would also need to consider the heat flow between system and environment (see, e.g., refs. [78,79] for a comparison of different definitions for the total entropy). At the beginning of the experiment, the system and environment are decoupled. If we further assume vanishing means and an unsqueezed initial state with σ(0) = diag[1/2, 1/2], the system is in its welldefined ground state, hence S(0) = 0. At late times, for comparison, the interaction with the environment has increased the uncertainty in the state of the system with respect to its initial (isolated) value, and we have that the entropy reaches its equilibrium value. Equation (19) is a monotonic function of u, but not necessarily of time t, as u is not a monotonic function of time. For a discussion of the von Neumann entropy over the course of the nonequilibrium evolution, we refer to ref. [18]. At late times, however, it becomes a constant depending on the system parameters as well as the coupling strength.
For instance, in the case of strong system+bath coupling, where u 1/4, we have that ∆S ∼ 1 + log[u] varying logarithmically with u. Hence, the uncertainty relations derived in Equation (A14) (see also Equation (18)) are inherited also by ∆S in exactly the same hierarchy. In this way, the fluctuation-dissipation inequality provides a lower bound to the entropy production in the system. Incidentally, that bound exceeds the bound one would expect on the basis of the Robertson-Schrödinger inequality. By implication, similar arguments can be drawn for the effective temperature β(t) −1 at late times.
Lastly, we comment on the consequences of the fluctuation-dissipation inequality on the Shannon entropy in the phase-space (Wigner) representation, i.e., the entropic uncertainty relation [80,81]. For nonvanishing second-order cross-correlation and Gaussian systems with positive Wigner function, the entropic uncertainty relation can be formulated on the basis of the system's Wigner function [82,83].
where, in the second line, we inserted the expression for the Wigner function in Equation (D6) and performed the Gaussian integrals. Equation (22) gives the Shannon entropy of the Wigner distribution and can be particularly useful from the information theoretical point of view in multipartite Gaussian systems [84]. Even though the difference is rather subtle in our single-oscillator case, we remark that h(W r ) is not primarily connected to the von Neumann entropy (i.e., a Rényi-1 entropy S) but rather the generalized concept of the Rényi-2 entropy (compare also to Equation (20)). We refer to [84] and references therein for further details. In Figure 1, we report a numerical evaluation of the uncertainty and the (positive) Wigner entropy over the course of the nonequilibrium evolution as well as study the distinct impact of quantum and thermal fluctuations. Clearly, from the Robertson-Schrödinger inequality, , which is the traditional statement of the entropic uncertainty relation. For instance, in our case, at t = 0, where u(0) = 1/2 corresponds to a pure Gaussian state, the inequality is saturated, h(W r (0)) = log[eπ]. At finite times, however, even though the system's Wigner function remains Gaussian, the system+bath coupling introduces additional fluctuations into the system so that the Wigner entropy h(W r ) exceeds its lower bound. Here, for comparison, the fluctuation-dissipation inequality in Equation (A13) can provide a more precise statement. Indeed, since it is equal to or exceeds the conventional Robertson-Schrödinger relation [Equation (A14)], the fluctuation-dissipation inequality can be used to extract information on the system+bath coupling from the entropic uncertainty relation, even when no exact solution (as is the case for Gaussian systems) is available.  (22), normalized to the minimal uncertainty log[πe]) as a function of time (measured in multiples of the dissipation rate). We employ the bath spectral density of Equation (E3) and use parameters ω0/γ0 = 10, Λ/ω0 = 10, and dimensions where = 1 as well as choose σ 0 = diag[1/2, 1/2] for the initial conditions. The lower dashed line gives the lower bound prescribed by the fluctuation-dissipation inequality (FDI) at late times (Equations (12), (13) and (22)), which can be connected to quantum fluctuations in the coupled sys-tem+bath system. The upper dashed line gives the exact late time limit of the Gaussian evolution (Equations (22) and (A12)), which also includes thermal fluctuations (see Section IV A). (bottom) Late-time quantum uncertainty [Equations (9) and (A13)] as a function of βω0, i.e., a measure of the respective impact of quantum or thermal fluctuations. For finite system-bath coupling (black, solid line; ω0/γ0 = 1), the uncertainty always exceeds the minimal bound of 1/4 given by the Robertson-Schrödinger equation and saturates the fluctuation-dissipation inequality for βω0 1 (gray, horizontal, dashed line). This discrepancy fades for smaller coupling (gray, solid line; ω0/γ0 = 100). Additionally, for βω0 1, thermal fluctuations start to prevail over the quantum fluctuations, and the more accurate bound (comparing to the FDI) can be provided by the thermal fluctuationdissipation inequality (TFDI; gray, dashed, nonhorizontal line; see Section IV A).

A. Thermal Fluctuation-Dissipation Inequality
So far, our results are derived from mathematical and conceptual considerations, i.e. the hermicity of the noise operator and the self-consistency of the interaction. We now turn our full attention to the specific case of quantum Brownian motion and derive a thermodynamic uncertainty relation for the nonequilibrium energy current between system and environment. To this end, we write the dissipation and the noise kernel in the convenient form (see, e.g., refs. [85,86]).
where I(ω) is the bath spectral density and β −1 = k B T the temperature of the environment with k B the Boltzmann constant (not to be confused with the effective parameter used in Equation (19)). The bath spectral density is an odd function in frequency. Moreover, using the relation ωγ(ω) = 2Im[Γ θ (ω)], we have that Γ θ (ω) is an analytic function in the upper complex frequency plane with singularities, i.e., the physical resonances of the system, that are symmetrically distributed with respect to the complex frequency axis. For each singularity at ω n with Im[ω n ] < 0, solving Γ −1 θ (ω n ) = 0, Γ θ (ω) also features a singularity at −ω * n (usually referred to as crossing relation). We note that this allows for (degenerate) purely complex solutions. Hence, ωγ(ω) is in general not analytic in the complete complex frequency plane. Instead, due to the Im[·] in Equation (24), for each singularity at ω n , the spectral density I(ω) shows the pair {ω n , ω * n } as poles in the complex frequency plane. In Appendix E, we summarize the properties of dissipation kernels for some common bath spectrum models.
Due to the simple form of Equation (23), a Fourier transform is well-defined at all times, and we obtain the exact relation.
This is the celebrated (local) fluctuation-dissipation theorem [30,64] and means that for the bath degrees of freedom in quantum Brownian motion, detailed balance is fulfilled at all times. We note that the latter equality in Equation (25) follows as a special case for zero temperature since coth βω/2 → sgn[ω] for β → ∞. Indeed, for vanishing temperatures, the equality holds even in time domain.
The difference per frequency between ν and γ stems from the thermal occupation of the field. Moreover, using that coth[xω] ≥ [xω] −1 (ω > 0), we can relate in frequency domain ν(ω) ≥ (2/β)γ(ω). The corresponding relation in time domain does not generally hold since the Fourier transform [Equation (23)] does not necessarily preserve order. Instead, we can specify the notion of the fluctuation-dissipation inequality in Equation (12) by using our knowledge of the detailed balance of the bath degrees of freedom in quantum Brownian motion, i.e., For brevity, and in order to discern Equations (12) (second line of the previous equation) and (27), we denote the first line of the previous equation as the thermal fluctuation-dissipation inequality. We emphasize that the first line of Equation (27) is not just a hightemperature approximation of the second line. When comparing with Equation (12), the thermal version of Equation (27) is valid for temperatures larger than the resonances of the system+environment composite. It provides a more accurate estimate than (12) does in the classical, high-temperature regime. In this regime, it has the advantage of providing the tighter lower bound on the system's fluctuations. The quantitative advantage, however, comes at the price of weakened generality. Indeed, at lower temperature, it needs to be replaced by the inequality in Equation (12) which gives the absolute quantum limit. We finish this paragraph by considering some concrete examples for the thermal fluctuationdissipation inequality as well as deduce its consequences on the fluctuations of the system degrees of freedom.
While the fluctuation-dissipation inequality [Equation (25)] gives the absolute quantum limit of the fluctuation's correlations and provides the most accurate bound for large frequencies (small time delays τ ), the thermal inequality (Equation (27)) puts the focus on the environment's temperature and becomes more accurate at small frequencies (large time delays τ ). As expected, the thermal fluctuation-dissipation inequality is hence connected to the classical thermodynamic limit. For analytical details on how the noise kernel is affected by different types of bath spectral densities, we refer to Appendix F.

B. Numerical Results and Quantifying Error
We report the numerical results for the autocorrelation of the momentum operator using the Λ-model for the bath spectral density (see Equation (E3) and Appendix F for details on the numerics), i.e., Q (t)Q(t ) s , in Figure 2. We take the momentum correlation to illustrate our findings as they play a role in the next section dealing with the energy flow between system and environment. As we anticipated from the expressions for the noise correlations ν and the response function χ, the autocorrelation is maximal if t = t and decays exponentially for increasing time delay |t − t | (top of Figure  2). Per construction (see Equation (15)), the momentum correlator vanishes if either t or t is smaller than zero, and the asymmetric shape with respect to the value at t = t stems from the fact that we have set t i = 0 as the initial time of our experiment. For comparison, if we were to set t i → −∞, the shape would be fully symmetric. The momentum correlation at equal times increases over time and reaches its late-time limit for times t 1/γ 0 (see Equation (A12)) (bottom of Figure 2). As expected from Equation (27), at all times, the momentum correlation exceeds the value prescribed by the thermal fluctuation-dissipation inequality, where we replace ν → (2/β)γ in the evaluation of the time integrals (dashed line, bottom of Figure 2).
Starting from the microscopic model in Equation (23), the difference between a regime that is dominated by quantum fluctuations and a regime that is dominated by thermal fluctuations becomes immediately clear from the coth function weighted by the bath spectral density. From this perspective, the thermal fluctuationdissipation inequality in Equation (27) may appear deceptively simple. However, when concrete examples for realistic situations are considered, where the coth function and the bath spectral density are buried in many layers of modeling and experimental circumstances, the situation can be less clear. Sometimes, it may seem legitimate to neglect certain types of fluctuations in the interaction, and hence it becomes interesting to quantify the resulting error. To this end, let us consider the difference between the exact momentum fluctuations and their lower limit in terms of the thermal fluctuationdissipation inequality, i.e., which is strictly positive by means of Equation (27). The previous relation reveals quite intuitively the physical meaning of the thermal fluctuation-dissipation inequality. The fluctuations of the system can be connected to their thermal (∼β −1 ) and quantum (∼ ω/2) nature. Both are interwoven (∼coth[ βω/2]) into the non-Markovian and nonequilibrium dynamics of the evolving system by means of the time and frequency integrals in Equation (28). Furthermore, they are weighted by the properties of the system and the environment (χ and γ) which set the relevant energy scales for the quantum fluctuations, i.e., ω 0 and Λ, γ 0 or Ω, κ, ρ(r a ) for the model in Equation (E3) or (E5) (see Appendix E), respectively. The thermal fluctuation-dissipation inequality now puts a special emphasis on the thermal (classical) fluctuations of the interaction, i.e., approximates the hyperbolic cotangent in the frequency domain (Equation (23b)) by its value for small arguments.
Since there is no way to circumvent quantum mechanics, the corresponding fluctuations always come on top, and hence the thermal fluctuation-dissipation inequal-ity provides a lower bound for sufficiently high temperatures (see discussion after Equation (27)). As a consequence, in the extreme limit of vanishing temperature, the thermal fluctuation-dissipation inequality becomes trivial and the general fluctuation-dissipation inequality or the Robertson-Schrödinger inequality (see, e.g., Equation (A14)) provide a sharper and physically nontrivial bound, since their focus lies on the quantumness of the system. Even though such considerations are well-known from equilibrium physics or nonequilibrium steady states, our considerations show that such concepts can to some extent be translated into the full nonequilibrium evolution of the system. In other words, Equation (28) could be understood as a measurement of the "quantumness" of the nonequilibrium dynamics of a quantum Brownian particle that also takes the particular bath spectral density into account. Indeed, for illustration, let us consider the Markovian and late-time (equilibrium) limit of Equation (28). Using the Λ model (Equation (E3) and Appendix E) in the limit Λ → ∞ (formally yields the Markovian limit) and additionally assuming that γ 0 ω 0 , the polarizability |α(ω)| 2 ∼ πδ(ω 2 − ω 2 0 )/(2γ 0 ω) and we obtain lim t→∞ ∆(t) = ([ βω/2] coth[ βω/2] − 1)/β → ω 0 /2 as β → ∞. The fluctuations of the (momentum) operator are purely quantum. For a numerical evaluation of the role of the thermal fluctuation-dissipation inequality on the uncertainty function of the system, we refer to Figure 1.
For a more advanced example, we can extend the previous discussion to the interaction of an atom with the material-modified electromagnetic field. Considering alkali-metal atoms and conducting macroscopic bodies, thermal fluctuations are usually much smaller than any of the system's resonances [87]. If the atom moves with constant velocity in the vicinity of a macroscopic body, it experiences a decelerating force, the so-called quantum friction [88]. Even at zero temperature, the motion-induced Doppler shift then instigates additional low-frequency fluctuations into the system that, for simplicity, can be understood as mimicking certain aspects of thermal fluctuations [51,89]. The important point is that, if one assumes that equilibrium can be established locally, even though the system is in a nonequilibrium state [90], the power incoming into the system can be written in a form very similar to Equation (28) (see Equation (11) in ref. [50]). In the particular example, a nonvanishing ∆(t) hinted towards the negligence of important fluctuations in the power spectrum [51] and could be shown to be a serious defect in the underlying statistical modeling of the interaction [50]. One is hence often well-advised to minimize ∆(t) in (quantum) fluctuation-induced systems. and normalized to its equal-time correlation at t = γ −1 0 . At t = γ −1 0 , the apparent kink is owed to the numerical resolution in time. The curve is smooth. (bottom) Equal-time correlation normalized to its late-time limit (solid). Lower bound prescribed by thermal fluctuation-dissipation inequality (dashed). The difference between the two is given by ∆(t) (Equation (28)).

V. NONEQUILIBRIUM CURRENT, ENERGY FLOW, AND CURRENT FLUCTUATIONS
In the literature on thermodynamic uncertainty relations, it has become customary to use a nonequilibrium current operator to quantify the degree of nonequilibrium, namely, the thermodynamic uncertainty relations provide a lower bound on the fluctuations of that current operator [9]. Quite generally, we can define the current operator aŝ with an, in principle, arbitrary function f that, for simplicity, is not explicitly time-dependent. We use the anticommutator here in order to ensure the real expectation values in our later examples. Various examples of such can be found in the literature (see Section I). As the particle is constantly exchanging energy with its surroundings, we could for example choose to consider the net power entering and exiting the particle, i.e., where the first term corresponds to the incoming power and the second term to the outgoing power. For simplicity, we focus on the contributions to the power connected to the noise operator (P in and P fluc out ), as the contributions connected to the initial conditions simply decay into the environment over time (see Equation (14b) and discussion below). Given the specific form off , due to the selfconsistency of our system, any autocorrelation ofĴ reduces to calculating the moments of the (Gaussian) noise operatorξ and solving the subsequent time integral [see Equation (11)]. For the net power the system exchanges with the environment, we obtain We remark that the clear splitting in different contributions of energy flows is connected to the linearity of the system. By using the symmetric average, we can render the individual parts as physically meaningful. This has been extensively discussed in the literature, and we refer the reader to refs. [32,[91][92][93] for further details.
In Figure 3, we report the outgoing power connected to the system's fluctuating dynamics using the bath spectral density of Equation (E3) (second part of Equations (31), see also Equation (17)). At late times, where the system can equilibrate, it balances the ingoing power (first part of Equation (31), see also Equation (14a)), and there is no net transfer of energy between system and its environment on average [94,95]. During the full course of the nonequilibrium evolution, the corresponding expression using the thermal fluctuation-dissipation inequality [Equation (27)] provides a lower bound to the outgoing power. Since the properties of the momentum correlation are inherited by the power, this behavior can be understood in full similarity to our discussion of the momentum correlations (see Equation (28) and discussion below).

A. Generalized Current Fluctuations
In the context of the thermodynamic uncertainty relation, it is more interesting for us to consider the fluctuations of the generalized current (see Section I). As an illustration and to keep the expressions transparent, we focus on the incoming power only and consider the total input power, i.e., On average, as the integral kernel of the previous equation approaches a constant at late times (see Section III), the mean current is equal to the incoming instantaneous power, i.e., Ĵ s /t → P in for t → ∞ (we use the time average for convergence following Refs. [16,96]). Indeed,  (17) and (31)) as a function of time in multiples of the dissipation rate γ −1 0 . Parameters are chosen as in Figure 2. We normalize to the expression for the ingoing power at late times Pin(∞) [Equation (14a)] in order to indicate the balancing of the two at equilibrium (solid line). We further report the corresponding expression using the thermal fluctuation-dissipation inequality, i.e., replacing ν → (2/β)γ in the numerical evaluation, which is always smaller than the full expression (dashed). Lastly, we give the upper estimate of Equation (42) (dotted).
For finite times, on the other hand, using Equation (11), we obtain This result is exact for all times of the nonequilibrium evolution and we note that we do not take the symmetric average on the right-hand side of the previous line. In order to employ our results on the (thermal) fluctuationdissipation inequality, we need to reorder Equation (34) using the commutators.
[ξ(x),Q(y)] = 2i which yields for the variance of the generalized current operator The previous line might appear more complicated than Equation (34), but it enables us to dissect the underlying physics. Firstly, the foremost term of Equation (36) is related to the incoming power P in (t) in the frequency domain.
Secondly, the last two lines of Equation (36) feature two terms that are given by average commutators only (last term in brackets in the last two lines). These are temperature-independent (see Equation (35) in combination with Equation (12) and discussion below) and are hence solely dependent on quantum fluctuations. In the high-temperature (classical) limit, where thermal fluctuations dominate, they can be ignored.
Thirdly, Equation (36) contains a number of terms that are given by a combination of a symmetric average and an average of a commutator. For times larger than the typical dissipation time scale of the system (∼γ), these can be expected to become exponentially small due to symmetry reasons under the integral.
Lastly, the remaining term (first term of the second line) of Equation (36) is even more interesting from our perspective, as the noise kernel is evaluated over the two-time correlations of the momentum operator, and we can determine a lower bound by means of the thermal fluctuation-dissipation inequality [Equation (27) where we used that the integral kernel is invariant with respect to replacing x ↔ y. Again, at later times, where the initial jolt has been damped and the system is approaching the steady-state, we can employ the thermal fluctuation-dissipation inequality (Equation (27)) and recall the expression for the outgoing power P fluc (17)], in order to write Collecting the results from the previous points, we can formulate the relation which is a new form of thermodynamic uncertainty relations: the fluctuations of the generalized current are greater than or equal to the (time-averaged) outgoing power times a thermal factor of 2k B T . From the perspective of the Brownian particle, the latter can be interpreted as dissipation/loss of energy to the environment. Increasing the accuracy of the incoming power comes at least at the price of the time-averaged dissipated energy (compare with the original formulation of the TUR in ref. [1]). Establishing an unambiguous connection to entropy production in the complete system-although rather commonly invoked in conventional formulations of TURs-is nontrivial for general non-Markovian quantum systems under fully nonequilibrium conditions. We refer the reader, e.g., to our first work [18] and references therein for more details. Instead, we outline in Section VI how our formalism can be applied to steadystate heat transfer where an unambiguous connection to entropy production is possible. To derive Equation (41), we had to assume that (i) the transitional dynamics is starting to settle, i.e., we work at times larger than the characteristic damping scales of the system+environment composite approaching the steady state (equilibrium in the current case), and (ii) that quantum fluctuations can be ignored. In any other situation, one needs to exercise particular care and is probably better advised to use Equation (36).

B. Non-Markovianity of the Damping Kernel
Remarkably, our thermodynamic uncertainty relation connects the fluctuations of the nonequilibrium current connected to the incoming power operatorĴ in to the power leaving the system P fluc out . It thereby establishes a statistical statement on the energetic interaction between system and environment in the course of their equilibration process. Only at late times, where P out → P in , the incoming power itself can be related to the lower thermodynamic bound of its corresponding fluctuations. Since we could provide an exact relation [Equation (36)], the inequality can be refined by (i) including the additional terms from Equation (36) or (ii) by including higherorder corrections from the noise correlation (Equation (F2)).
Furthermore, we would like to comment on the Markovian and high-temperature limit of Equation (40). Starting from the Λ-model for the bath spectral density (Equation (E4) and Appendix E), the Markovian case is achieved in the limit Λ → ∞, where γ Λ (τ ) → 2γ 0 δ(τ ). In this case, the fluctuating part of the outgoing power reduces to P fluc out (t) → 2γ 0 Q 2 (t) s and is fully determined by the momentum fluctuations of the system. For comparison, we have seen earlier ( Figure 2) that it is generally true that Q 2 (t) s ≥ Q (t)Q(τ ) s , i.e., the momentum correlations become maximal at equal times. Hence, a general upper bound for the incoming power is given by Choosing once again γ = γ Λ in the limit Λ → ∞, we would reproduce the Markovian result for P fluc out . However, for the non-Markovian case with finite Λ, we instead obtain P fluc . This is less than or equal to the Markovian limit, and the equality is only achieved at late times (t → ∞) (see Figure  3). In other words, our non-Markovian thermodynamic uncertainty relation (Equation (40)) can be expected to provide a lower bound than its Markovian counterparts. Furthermore, we remark that Equation (42) is sensitive to the particular regularization scheme [18,74].
Lastly, let us once again emphasize that it is not the particular example of the nonequilibrium current we are interested in, but Equation (32) was rather chosen for its simple form. The important insight is that we could deduce a relation from a chain of arguments starting from the hermicity of the noise operators and the selfconsistency of our system+bath dynamics over the uncertainty relations and the fluctuation-dissipation inequality, all the way to what we coined the thermal thermodynamic uncertainty relation. In this way, the thermodynamic uncertainty relation finds its clear footing in the microscopic stochastic properties of open quantum systems.
We conclude our discussion in the next section by discussing one more example and connecting our formalism to the related problem of steady-state heat transfer.

VI. NONEQUILIBRIUM STEADY STATE AND CONNECTION TO HEAT TRANSFER
At late times, the energy lost to the environment and the energy inflow into the system balance on average such that the system can equilibrate [74]. The net current (Equation (31)) vanishes in the mean, i.e., Ĵ s = 0. This, however, is not true for its fluctuations that can remain finite even in equilibrium. To see this more clearly, let us again for simplicity consider the current connected to the incoming power only (Equation (32)). From Equation (40), we immediately find the expected late-time limit (2/β) Ĵ fluc out s → (2/β) lim t→∞ P out (t) as found in Equation (18) The formalism we used above can be readily extended to a situation of multiple particles connected to multiple heat baths at different temperatures [94,95,98]. At late times, the system then does not equilibrate but rather reaches a nonequilibrium steady state, where it mediates a constant average heat transfer between the heat baths [99,100] (see also refs. [43,94,95,101] for a modern take on the topic). The fluctuations of such a nonequilibrium steady-state current have been calculated quite generally in refs. [102][103][104], and a particular focus on the thermodynamic uncertainty relation was put just recently in refs. [5,16]. Comparing with our full nonequilibrium form in Equation (34), it can be instructive to explore how the thermal fluctuation-dissipation inequality affects the thermodynamic uncertainty in steady-state heat transfer. It again turns out to be the hidden fundamental principle determining the macroscopic statistical behavior.
For simplicity, we consider two interacting particles connected to two heat baths individually. Furthermore, we consider the simplified example of a delta-functionshaped dissipation kernel γ(t) ∼ γ 0 δ(t). Following the approach and notation of ref. [94], we then obtain a system of two coupled quantum Langevin equations (both particles feature equal mass m = 1), i.e., where Ω = is the coupling matrix,χ = (χ 1 ,χ 2 ) T the vector operator for the two quantum degrees of free-dom, andξ = (ξ 1 ,ξ 2 ) T the corresponding noise operator of the respective quantum baths with ξ Here, δ ij is the Kronecker-delta and β 1|2 the inverse temperature of the respective quantum baths. The heat flow can be measured at various different points of the system. However, in the nonequilibrium steady state, the particular choice does not influence the result anymore. For convenience, we define the nonequilibrium current as the energy transferred from particle 2 to particle 1, i.e., We note that earlier references, e.g., refs. [5,52,102], considered the energy transferred from one reservoir into the coupled system, which requires a somewhat different calculation. In Appendix G, we explicitly show that the same result can be obtained with the point of measurement in Equation (44). After a lengthy but straightforward calculation, we find with the transmission function T (ω) = (γ 0 ωσ) 2 Π ± [ω 2 0 − ω 2 ±σ±2iγ 0 ω] −1 (the product runs over all possible combination of ±1) This result, in different contexts leading to various transmission coefficients, has been found by many authors before, see, e.g., [52,96,102,103,105] and references therein.
Recently, it was found [5,16] that steady-state heat transfer is in accordance with the thermodynamic uncertainty relation by means of an expansion of the thermal functions in Equation (45), independent of the concrete form of the transmission coefficient (see Equations (13)- (15) in [5]). To this end, we use that coth[x/2] = 1 + 2N (x) with N (x) = [e x − 1] −1 as the bosonic occupation number. If we further neglect the strictly positive first term in Equation (45), we follow ref. [16] and readily find Although the previous result could have been similarly obtained from purely mathematical arguments, it becomes clear from our present discussion that the thermodynamic uncertainty relation in steady-state heat transfer is deeply rooted in the fluctuation-dissipation inequality, and hence the causality and the statistical properties of the underlying Hamiltonian. Indeed, in our simple case of heat transfer with time-local dissipation kernel, the thermal fluctuation-dissipation inequality is practically a statement on the properties of the hyperbolic cotangent in the second line of Equation (45) (see discussion below Equation (27)), which was the necessary step in order to derive the inequality in Equation (46). Lastly, we comment on the connection between the thermodynamic uncertainty relation and entropy production in the system. In the nonequilibrium steady state, the density matrix of the system of interest becomes stationary per definition, as does the von Neumann entropy of the system. However, due to the constant energy flow mediated by the system, the total entropy of the sys-tem+bath increases at a constant rate [78,106] and the average entropy production rate in the total system was identified via ∂ t S s = (β 2 −β 1 ) −1 Ĵ s [5,16]. In the case of quantum Brownian motion, where the system can equilibrate at late times, the entropy production vanishes. It is rather the total change in entropy of the system over the course of the equilibration process (Equation (19)) that quantifies the uncertainties. From the perspective of the thermodynamic uncertainty relation, it can then be more convenient to use the power flow between system and environment in order to quantify dissipation (Equation (40)).

VII. CONCLUSIONS
We systematically explored the emergence of macroscopic thermodynamic uncertainty relations [1] for generalized currents from microscopic uncertainty principles of interacting open quantum systems. Using the generic example of quantum Brownian motion, we set the microscopic framework by adopting the fluctuationdissipation inequality [49] for the system's microscopic quantum degrees of freedom (canonical operators) [18]. The fluctuation-dissipation inequality is solely based on reasonable physical assumptions, such as the hermicity of the involved operators and a causal interaction between system and bath. It provides the inexpugnable lower bound for the system's microscopic uncertaintyincluding the Robertson-Schrödinger inequality as the nonequilibrium standard quantum limit. However, when thermal fluctuations are dominating, the fluctuationdissipation inequality may not provide the most accurate lower bound. Moreover, although the fluctuationdissipation inequality is a fairly fundamental statement, it is not sufficient to fully explain the emergence of a thermodynamic uncertainty relation for macroscopic thermodynamic quantities. Instead, we needed to specify the statistical properties of the bath in order to extract information on the temperature dependence of the system's fluctuations. Based on a microscopic model for the bath spectral density, this led us to formulating a thermal fluctuation-dissipation relation that is valid at high temperatures, which provides a tighter bound than the fluctuation-dissipation inequality. This means that it can provide a better estimate in the thermal-energydominated regimes (see Equation (27)). This thermal fluctuation-dissipation inequality which applies to the classical domain enables us to compare with popular TURs in the literature [1,9,13,107], where the hightemperature limit is a requirement rather than a particular limit. When applied to the generalized current, we show how it formally leads to a thermodynamic uncertainty relation. In essence, it is the combination of the statistical properties of the bath and the causality of the system+bath interaction that is inherited by thermodynamic quantities (e.g. generalized currents) and can be seen as the microscopic origin of thermodynamic uncertainty relations in linear open quantum systems.
As an instructive example, we examine the power entering the Brownian particle as the generalized current and find an exact expression of the corresponding current-current correlations. At high temperatures, applying the thermal fluctuation-dissipation inequality, in Equation (41) we showed that the corresponding fluctuations are always equal to or greater than the average outgoing power weighted by a thermal factor 2k B T . For temperatures that are low with respect to the dominating energy scales in the system, quantum corrections need to be taken into account. Our result fully includes back-action from the environment onto the system and respects the particular anatomy of realistic bath spectral densities, albeit limited to the linear regime. It is straight-forward to extend our formalism to multiple particles connected to several heat baths. This enabled us to check the consistency of our result by rederiving a thermodynamic uncertainty relation already known from the literature on steady-state heat transfer.
Our analysis sheds some light on how thermodynamic uncertainty relations are deeply rooted in the quantum uncertainty relations and how they impact macroscopic observables. It further provides some basic tests and requirements in the form of inequalities that any trustworthy physical result within the realm of our assumptions must comply with. In the following, we provide details on the properties of the Langevin equation and its solutions discussed in Section II.
The response function χ(t) obeys the homogeneous equation with the initial values χ(0) = 0,χ(0) = 1. We note that the response function is causal χ(τ ) = 0 for τ < 0. Furthermore, we can transform the convolution in Equation (A1) to an algebraic expression in Fourier domain and obtain for the response function the expression in Equation (8). Defining the matrix The general time-dependent solution for the quantum degree of freedom reads From this, we can derive the second-order correlations in Equation (9).

Fluctuation-Dissipation Inequality
We begin our discussion with the relation between ν and γ characterizing the environment. For any quantum mechanically well-defined operatorÔ(t) and scalar product · , we can state that This implies that the quantum average ofÔ † (t)Ô(t) is a semipositive definite function of time in the sense of [49,65].
We note that, ifÔ is hermitian and f a real function, the r.h.s. of the previous line is antisymmetric with respect to interchanging τ ↔ τ and hence becomes trivial. For nonhermitian operators, that is not necessarily the case. Equation (12) is a slightly generalized form of the fluctuation-dissipation inequality put forward in ref. [49]. If we now setÔ =ξ, Equation (12) simply states that the hermicity ofξ translates into a correlation function ν that is a positive semidefinite function in the sense of Equation (A5). Alternatively, setting f (τ ) = 1, substituting x = τ − τ and using We find that the fluctuations integrated over the full span of time always are equal to or exceed the corresponding first moment of the time delay, i.e., Additionally, we can make a statement on the relation between ν and γ in frequency domain. To this end, without loss of generality, we extend the integral boundaries in our definition of the scalar product in Equation (A5) to ±∞. It is then possible to show that for stationary noise and any instance of time, arbitrarily far from equilibrium, the noise spectrum of the environment is always equal to or exceeds the spectrum of the dissipation kernel weighted by the factor ω (see ref. [49] for details) which is stated in Equation (13). We note that the generalization to multivariate noise is straight-forward [49,50]. Equation (A8) explicitly shows that detailed balance is not necessarily fulfilled in nonequilibrium.

Robertson-Schrödinger Inequality
In terms of the covariance matrix ∆Q 2 s , the Robertson-Schrödinger inequality can be stated as We can specify the previous line to our specific situation. First, we note that the cross correlations can be written as where the subscript denote the corresponding components of the matrix and we used that, at initial time t = 0, the system and environment fulfill the usual canonical commutation relations [q(0),q(0)/ω 2 0 ] = i . Here, the extra ω −2 0 is due to the proper definition of the canonical momentum. Combining Equations (A9) and (A10), we obtain the time-dependent Robertson-Schrödinger function for strongly coupled system-bath dynamics It is quite interesting to consider the late-time limit of the previous inequality. To this end, we note that almost all components of the first three terms of Equation (A11) are given by elements of the linear response tensor X(t).
Since the latter obeys the damped Equation (A1), it can generally be expected to vanish in the limit t → ∞. The remaining term lim t→∞ t 0 dτ t 0 dτ χ(τ )χ(τ )ν(τ −τ ) = 0 also vanishes since we integrate an odd function over a symmetric integral. Only the last term in Equation (A11) prevails at late times, such that the Robertson-Schrödinger equation evolves according to lim t→∞ ∆q 2 (t) s ∆q 2 (t) s (A12) For comparison, at late times, the fluctuation-dissipation inequality applied to the variance of position and momentum operators directly yields The connection between the fluctuation-dissipation inequality (FDI) and the Robertson-Schrödinger (RS) uncertainty relation is provided by the Cauchy-Schwarz (CS) inequality for integrals [108], which further bounds Equation (A13) as In fact, it can be shown that the last line of the previous equation, i.e., the lower bound of the Robertson-Schrödinger equation, reproduces the conventional Heisenberg bound, since dω 2π ω 2 |α(ω)| 2 |γ(ω)| = 1/2 [18]. We kept the frequency representation as it makes the connection to the fluctuation-dissipation inequality more transparent. In this way, at late times, the fluctuation-dissipation inequality not only reproduces the Robertson-Schrödinger equality as a special case, but actually provides a stronger bound.
From the previous uncertainty relations, especially the FDI in the second line of Equation (A14), it becomes clear that the strong coupling between the system and its environment leads to fluctuations that can exceed the usual Heisenberg bound. Indeed, the Heisenberg uncertainty relation should be restored in the limit of vanishingly small system+bath coupling which is implicitly encoded in the magnitude of γ. To show this explicitly, we focus on the late-time limit for simplicity and note that γ(τ ) is a real function that is even in its argument such that γ(ω) = 2Re[γ θ (ω)]. In the limit of small coupling, we then have that Here, 2ωIm[γ θ (ω)] can be seen as system's resonance frequency due to the coupling with the environment. In leading order coupling, this is subleading with respect to ω 0 , so that we obtain lim t→∞ γ→0 This is nothing but the Heisenberg uncertainty relation, and we refer to, e.g., refs. [27,28] for further details as well as ref. [18] for a finite-time version of Equation (A14).

Appendix B: Wick's Theorem for Thermal States
In the perturbative treatment of in-out quantum field theory, the series is often expanded in terms of the vacuum n-point correlation functions of free fields, that is, the vacuum expectation values of the time-ordered product of n free-field operators. The Wick's theorem states [109,110] that the time-ordered product of the free-field operators can be expressed as the sum of the normal-ordered product and the permutation sums of the products of pair-wise contractions. When the state is the vacuum state of the free field, the normal-ordered term vanishes and each contraction gives a Feynman propagator. Then, the Wick's theorem implies that the vacuum n-point correlation functions of the free field can be expanded by the permutation sums of the products of the corresponding two-point Feynman propagators. Thus, the perturbative quantum field theory reduces to the knowledge of the two-point Feynman propagator.
However, the thermal expectation value of the normalordered term in general does not vanish. So, we may wonder whether the thermal expectation values of the timeordered product of n free-field operators can be likewise expanded by the finite-temperature Feynman propagator. Various generalizations of the original Wick' theorem are formulated for finite-temperature field theory [111], in-in nonequilibrium field theory [112], and arbitrary initial state of the field operators [113].
for all permutations. Here φ (x 1 )φ(x 2 ) β , for example, is the thermal Feynman propagator, constructed from the field operator evaluated at the spacetime points x 1 and x 2 . In particular, the four-point function Note that the ordering of the operators matters. Equation (B13) implies that we can write This allows us to compute the energy current fluctuations.

Appendix C: Fluctuating Energy Flow
We provide the necessary steps to derive Equation (16), i.e., connecting the derivative of the mechanical en-ergy stores in the operatorQ explicitly with the incoming and the outgoing power. The nontrivial part of the relation is to connect the expression for P fluc out in Equation (17) to the form used in Equation (16). We start from Equation (17) and note that the threefold integral can be partially decoupled using for arbitrary functions f, g. Combining the previous lines yields The latter reproduces Equation (16).

Appendix D: Master Equation, Density Matrix, and Von Neumann Entropy
Following refs. [63,114] (see also Refs. [52,61,62,67,68,114] for details), it is possible to connect the solution to the quantum Langevin Equation (A3) to a master equation for the corresponding Wigner function W r (q,q, t) (q,q are now to be understood as coordinates and not operators). To this end, we start from the statistical relation where we average over the initial conditions Q 0 = Q(0) (see Equation (A3)) as well as the impact of the environment by means of the Feynman-Vernon functional integral (see refs. [56,63,114] for details). Equation (D1) yields for the characteristic function (Wigner function in Fourier domain Q ↔ k) where we use the matrix X defined in Equation (A2). Upon differentiation with respect to time and performing a Fourier transform, we obtain the master equatioṅ where we introduced the matrix Φ = −ẊX for convenience of notation and the gradient with respect to the coordinates, Q is to be understood as ∇ Q ≡ (∂ q , ∂q) T . Equation (D3) is the Hu-Paz-Zhang master equation with time-dependent coefficients [59]. We note that a direct connection to our initial quantum Langevin equation [Equation (7)] can be seen from recasting the latter in its time-local form [67,115], i.e., For given χ(t) and ν(t) (specifics of the system), the dynamics of the system in quasi-phase space is determined by the initial conditions as well as the environment-induced Gaussian covariance matrix. Choosing a simple Gaussian form for the initial Wigner function We can perform the Fourier transform exactly and obtain Note that d 2 Q W r (Q, t) = 1. Finally, for the density matrix in position representation, we need to perform yet another series of Gaussian integrals and obtain where we have definedσ(t) ≡ σ 0 (t) + σ(t) to shorten notation, and the subscript "qq" denotes the corresponding entry of the matrix. For vanishing mean values, i.e. q(t) = q (t) = 0, the previous line coincides with Equation (3.23) of ref. [19]. Such a result is quite remarkable and was extensively discussed in refs. [19,52] because the density matrixρ can now be rewritten in a Gibbs formρ with effective temperature Quantifying the coupling between system and environment and the effective time-dependent Hamiltonian Note that the Robertson-Schrödinger inequality is directly encoded in the temperature as well as in the corresponding partition function Combining the previous results leads to Equation (19) used in the main text [19,78].

Appendix E: Bath Spectral Density and Dissipation Kernel
We review the common bath spectral densities and their analytic properties (see Equation (24) and discussion below). First, we note that we can write (ω n = −ω * −n , ω n = ω R n − iω I n ) [116].
with complex constants c −n = c * n . We can then evaluate the dissipation kernel by means of the residue theorem We note that, for simplicity, we have implicitly assumed that there are no connected branch cuts in the bath spectral density. Let us introduce two relevant examples for the bath spectral density. One convenient example that is often used is to have I(ω) either featuring exponential [28] or Lorentzian damping [114,117] at the order of a cut-off frequency Λ. In the latter case, choosing so-called Ohmic damping (linear in frequency), we can write with γ 0 > 0 a phenomenological damping rate. We can identify ω R n = 0, ω I n = Λ, and c n = 2γ 0 Λ, and hence find Another, more physical example, occurs when the system is coupled to a bath that features collective oscillations at frequency Ω and width κ, where I Ω (ω) = ω π 2ρ(r a )κΩ 2 (Ω 2 − ω 2 ) 2 + κ 2 ω 2 . (E5) Here, ρ(r a ) is a geometric factor that gives the strength of the interaction and could, in principle, depend on the position of the particle. For instance, if the system is an atom in the vicinity of a macroscopic object, ρ(r a ) would be connected to the small frequency value of the local density of states of the material-modified electromagnetic vacuum [118]. For positions r a close to an interface, Ω would then be related to the surface plasmonor phonon-polariton frequency. See, e.g., refs. [119][120][121][122] for some applications of such a description. Furthermore, we note that Equation (E5) is substantially different from Equation (E3), as it features its very own damping mechanism described by the phenomenological dissipation rate κ that is connected to the width of Ω. Although it is possible to solve the frequency integral for γ(τ ) with the spectral bath density of Equation (E5) by means of an exact partial fraction decomposition, it is simpler to consider the limit κ Ω, which holds for many realistic cases, i.e., I Ω (ω) ∼ ωρ(r a ) π κ/2 ( κ 2 ) 2 + (ω − Ω) 2 + κ/2 ( κ 2 ) 2 + (ω + Ω) 2 .
The spectral features of the environment asymptotically show as an overall oscillatory behavior in the dissipation kernel. Lastly, it is relevant to mention that the bath spectral densities lead to an intrinsically non-Markovian dissipation kernel γ. Indeed, only in the limit of an infinitely large cut-off frequency Λ, the Ohmic damping turns to its (classical) Markovian value γ Λ (τ ) → 2γ 0 δ(τ ). When the bath features at least one finite collective resonance, the situation is even more intricate, as the Given the bath spectral densities in Appendix E, we compute analytical expressions for the noise kernel. We further give details on the numerical integration used for calculating correlation functions for the system degrees of freedom.
From the expressions in Equation (23), we can derive corrections to the inequality in Equation (27) up to arbitrary order. To this end, we employ the series expansion of the hyperbolic cotangent [126][127][128].
and integrate the frequency integral for the noise correlation in every order n, i.e., ν(τ ) = 2 β γ(τ ) + ∞ n=1 dω I(ω) ωe −iωτ ω 2 + (nη) 2 , (F2) where η = 2π/ β is the first Matsubara frequency [129,130]. We can solve the frequency integral quite generally by means of the residue theorem. For simplicity, we focus on isolated, simple poles of the bath spectral density in the complex frequency plane such that the spectral density takes the form where F is a positive constant and P (ω) an analytic function for complex frequencies (e.g., polynomial). We have seen earlier that, given a residue ω j , ±ω * j are also residues of the bath spectral density. For purely imaginary residues ω j , the set of additional solutions reduces toω * j = −ω j . In this way, P (ω) becomes an even function of frequency as does j (ω − ω j ). We note that I(ω) does not preserve causality under frequency integration (see also Equation (24)). Causality, as we see in the following paragraph, is mathematically established by the convergence criterion of the Fourier transform in Equation (23).
For illustration, let us consider the two bath spectral densities in Equations (E3) and (E5). Applying the residue theorem to the overdamped model in Equation (E3), we obtain which was, e.g., already reported in ref. [117]. where we assumed κ/2 Ω in the final expression. Together with the expression for the system's response function, an exact integration of the system's fluctuations Q (t)Q(t ) s becomes possible for our choices of the spectral bath density. To this end, we first note that the polarizability α(ω) features poles in the lower complex frequency plane only, in order to preserve causality. We can therefore evaluate the response function by means of the residue theorem and obtain where 1 ≤ j ≤ N is the number of roots ω j solving [ω 2 0 − ω 2 j ]P 2 (ω j ) − 2iω j P 1 (ω j ) = 0 (F7) and we have written the causal dissipation kernel as the fraction of two complex polynomials γ θ (ω) = P 1 (ω)/P 2 (ω) (see, e.g., Equations (E3) and (E6) in combination with Equation (23)). χ(τ ) is a real function due to the crossing relation which translates into ω j = −ω * k for k = j [97]. In other words, for purely complex resonances, as it is, e.g., the case in the model of Equation (E3), the response function becomes exponentially damped. In the other case, i.e., the model in Equation (E6), the system's response function generally behaves as the linear combination of damped oscillations with coupling constants (we follow the notation of ref. [117]).
We iterate that the generalization to include branch cuts in the bath spectral density is possible but requires more care in the limiting procedures. The properties of the field fluctuations ν and the system's response function χ are inherited by the fluctuations of the system degree of freedom, i.e., Q (t)Q(t ) s = t 0 dx t 0 dy χ(x)χ(y)ν(y − x + t − t ) and similarly for momentum and cross-correlations. For instance, in the case of the Λ-model, due to the relatively simple form of the response function χ and the noise correlations ν as a linear combination of exponentially damped functions, the corresponding time integrals can be solved analytically (see also ref. [117] for the case of t = t ). The only complexity arises from properly accounting for the time ordering between t and t . To this end, we note that only the modulus square of the time argument τ of the noise kernel (see Equations (F4) and (F5)) is relevant. For t > t , we use that the argument of ν becomes negative (positive) for y ≤ (≥)x − (t − t ) and split the y-integral into t 0 dy = x−(t−t ) 0 dy + t x−(t−t ) dy. For t < t , we use that the noise kernel ν is symmetric in its argument and we reduce the problem to the previous case t > t by exchanging t ↔ t . The subsequent integration is then straight-forward and we employ a computer algebra system to carry out the calculation. The Mathematica code [131] and the calculated data is available upon reasonable request. Lastly, we numerically compute the sum over the coth's complex thermal resonances (nη in Equation (F4)), which is converging rapidly (at least as n −2 ). Taken together, the previous reasoning is the basis for the numerical evaluations performed for the present manuscript.