Currents and fluctuations of quantum heat transport in harmonic chains

In the context of heat transfer, harmonic systems away from equilibrium have attracted a lot of attention over the last few years since the path from first-principle microscopic models to phenomenological results such as Fourier's law has turned out to be a formidable task [1–8]. Specifically challenging for quantum systems is the correct description of interactions between the system and environmental degrees of freedom. While in many situations, the assumption of weak system–bath interactions leads to powerful perturbative techniques such as master equations, the problem of quantum heat transfer seems to be particularly susceptible to the modeling of the correlations between the relevant system and environmental degrees of freedom. In fact, the underlying assumption of an unperturbed Gibbs state of the embedded system is only valid for vanishing coupling strengths, thus implying vanishing heat transfer. Weak system–bath interactions at low temperatures induce correlations between system and reservoir, which lead to non-negligible energy exchange [9] and produce unphysical results in conventional perturbative treatments [10, 11]. Moreover, the relaxation towards stationary states at low temperatures generally depends on the nature of the surrounding heat baths and, particularly, on non-Markovian dynamics [12–15]. Since harmonic systems allow for exact results at least in principle, they may on the one hand serve as non-trivial paradigmatic examples and, on the other hand, as starting points for more elaborate models.

Quantum heat transport through harmonic chains has been treated based on quantum Langevin equations which, for linear systems, leads to correct results for a large variety of couplings and temperatures [16–18]. However, while conceptually simple, practically these approaches focus on symmetrized bath-correlation functions to avoid the sampling of quantum time evolutions in the presence of complex-valued quantum noise. They are thus restricted to the study of symmetrized correlation functions, which reflect the dispersive part of the system and are connected with the fluctuation-dissipation theorem to the anti-symmetric correlation that constitutes the absorptive part [19]. These symmetrized correlations are sufficient to extract first moments of heat fluxes, but do not allow us to gain access to higher-order moments which are determined also by contributions from phase–space correlations that are not symmetrized [20, 21]. Further, these approaches represent the dynamics in terms of normal modes of a linear chain while the system–bath coupling remains local. Their numerical efficiency thus degrades substantially for sufficiently long chains. The aim of this paper is to provide a complete formalism which combines a non-perturbative treatment of the system–bath interaction with a high computational efficiency and allows us to study chains exposed to any damping or reservoir temperatures in the presence of external driving.

The formalism is based on a description of open quantum dynamics in terms of reduced-density operators as pioneered by Feynman and Vernon. The corresponding formally exact path integral expression allows for an exact mapping on a numerically more convenient stochastic Liouville–von Neumann equation (SLN) [22, 23]. The SLN is particularly beneficial for including external time-dependent forces and applies to any coupling strength and bath temperatures. Technically, system–bath interactions are exactly accounted for by introducing stochastic c-noises [24] the correlations of which reproduce the exact quantum bath–bath correlations. Here, we start from this very general formalism and adapt it to describe heat transfer through harmonic chains. Accordingly, the stochastic sampling can explicitly be performed to arrive at a deterministic equation of motion for the covariance matrix. The corresponding set of coupled ordinary differential equations can be solved in a straightforward manner even for a very large number of constituents in the chain.

Specifically, we apply the formalism to harmonic chains as shown in figure 1, where the system is out of equilibrium due to the coupling on reservoirs at different temperatures. The transient time evolution is investigated, as well as nonequilibrium steady states for heat fluxes and their correlations for a large variety of dampings, couplings, and temperatures. A high sensitivity of heat fluxes on damping strengths and intra-chain couplings are observed. Further insight is obtained by considering spatiotemporal heat-flux correlations which are compared with predictions from analytical Gibbs calculations. For stronger system–reservoir couplings, temperature gradients lead to substantial nonequilibrium effects in the correlation patterns. Beyond the scope of the present study is the extension to time-dependent driving to follow protocols to operate quantum heat engines also for stronger dissipation and in proximity to a possible quantum speed limit.

Zoom In Zoom Out Reset image size

2.1. Model

A safe starting point for treating open quantum systems is to consider the full Hamiltonian of the global system where energy is conserved. This Hamiltonian has three parts which reflect system, coupling, and reservoir, respectively:

$$\begin{eqnarray}&&H={H}_{{\rm{s}}}+{H}_{{\rm{I}}}+{H}_{{\rm{R}}}.\end{eqnarray} \tag{ 1 }$$

The system Hamiltonian is composed of N oscillators coupled by a quadratic nearest-neighbor, two-body potential with strength μ:

$$\begin{eqnarray}&&{H}_{{\rm{s}}}=\displaystyle \sum _{n=1}^{N}\displaystyle \frac{{p}_{n}^{2}}{2m}+\displaystyle \frac{1}{2}m{\omega }_{0}^{2}{q}_{n}^{2}+\displaystyle \frac{\mu }{2}\displaystyle \sum _{n=1}^{N-1}{({q}_{n}-{q}_{n+1})}^{2}.\end{eqnarray} \tag{ 2 }$$

The reservoir has the usual form that is used to model quantum Brownian motion of a particle coupled to a thermal heat bath: a collection of harmonic oscillators with positions ${x}_{\alpha }$, momenta ${p}_{\alpha }$, and masses ${m}_{\alpha }$ forming a quasi-continuum of reservoir modes with frequencies ${\omega }_{\alpha }$ (α, $\alpha ^{\prime}$ for reservoirs coupled to q₁ and q_N, respectively) [25]:

$$\begin{eqnarray}&&{H}_{{\rm{R}}}=\displaystyle \sum _{\alpha }\displaystyle \frac{{p}_{\alpha }^{2}}{2{m}_{\alpha }}+\displaystyle \frac{{m}_{\alpha }{\omega }_{\alpha }^{2}}{2}{x}_{\alpha }^{2}+\displaystyle \sum _{\alpha ^{\prime} }\displaystyle \frac{{p}_{\alpha ^{\prime} }^{2}}{2{m}_{\alpha ^{\prime} }}+\displaystyle \frac{{m}_{\alpha ^{\prime} }{\omega }_{\alpha ^{\prime} }^{2}}{2}{x}_{\alpha ^{\prime} }^{2}.\end{eqnarray} \tag{ 3 }$$

We assume that the first and last oscillators with index n = 1 and n = N are coupled to a thermal reservoir (see figure 1). The coupling to the baths with the strengths ${c}_{\alpha }$ is bilinear in the positions of the system oscillators and reservoir coordinates

$$\begin{eqnarray}&&{H}_{{\rm{I}}}={q}_{1}\displaystyle \sum _{\alpha }\,{c}_{\alpha }{x}_{\alpha }+{q}_{N}\displaystyle \sum _{\alpha ^{\prime} }\,{c}_{\alpha ^{\prime} }{x}_{\alpha ^{\prime} }+{q}_{1}^{2}\displaystyle \sum _{\alpha }\,\displaystyle \frac{{c}_{\alpha }^{2}}{2{m}_{\alpha }{\omega }_{\alpha }^{2}}+{q}_{N}^{2}\displaystyle \sum _{\alpha ^{\prime} }\,\displaystyle \frac{{c}_{\alpha ^{\prime} }^{2}}{2{m}_{\alpha ^{\prime} }{\omega }_{\alpha ^{\prime} }^{2}}\end{eqnarray} \tag{ 4 }$$

with two additional quadratic potential terms ('counterterms'). These ensure that the coupling to the reservoir induces no net potential when the reservoir is eliminated adiabatically. In the classical limit, this model corresponds to a purely velocity-dependent friction force.

One can show that the reduced system dynamics depends on the microscopic parameters in (3) only in an aggregate form, which is commonly denoted by a spectral density [26]

$$\begin{eqnarray}&&{J}_{r}(\omega )=\displaystyle \frac{\pi }{2}\displaystyle \sum _{\alpha }\displaystyle \frac{{c}_{\alpha }^{2}}{{m}_{\alpha }{\omega }_{\alpha }}\delta (\omega -{\omega }_{\alpha }).\end{eqnarray} \tag{ 5 }$$

Indexing each reservoir by the chain element $r\in \{1,N\}$ it is coupled to, all dissipative effects can be described by the reservoir correlation functions

$$\begin{eqnarray}&&{L}_{r}(t)=\displaystyle \frac{{\hslash }}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega {J}_{r}(\omega )\left[\coth \left(\displaystyle \frac{{\beta }_{r}{\hslash }\omega }{2}\right)\cos (\omega t)-{\rm{i}}\,\sin (\omega t)\right]\end{eqnarray} \tag{ 6 }$$

that depend on the respective spectral density and inverse temperature ${\beta }_{r}=1/{k}_{{\rm{B}}}{T}_{r}$.

For the spectral density of our model, we choose ${J}_{r}{(\omega )=m\gamma \omega (1+{(\omega /{\omega }_{c})}^{2})}^{-2}$ with a UV cut-off frequency ${\omega }_{c}$. This constitutes an Ohmic reservoir with damping strength ${\gamma }_{1}={\gamma }_{N}=\gamma$ for frequencies far below the cut-off ${\omega }_{c}$. With the choice of the spectral density and temperature, the thermally averaged motion of the chain is determined as a reduced dynamics by formally tracing out the reservoir. The imaginary part in (6) arises from time ordering; Weyl ordering yields only the real part.

This is sufficient when expectation values of energy currents are computed as, e.g. in [18], but the imaginary part is essential when current correlations are considered, which are given by fourth-order correlations of positions and momenta. Since the decomposition of higher-order correlations using Wick's theorem relies on time-ordered two-time correlations, the full complex expression (6) will be needed when we consider current–current correlations in section 5.

2.2. Reduced dynamics and stochastic mapping

The global Hamiltonian (1) uniquely defines the dynamics of the global density operator in terms of the Liouville–von Neumann equation. Since this is a high-dimensional system, a reduced dynamics in terms of the system density operator is needed. However, defining a reduced dynamics fully consistent with the global dynamics is difficult, particularly when the damping is large enough to modify the reduced equilibrium state or when the system is too complex for dissipation to be analyzed in terms of the system's level structure.

Path integrals involving a Feynman–Vernon influence functional [27] provide an exact approach in these settings. However, they cannot easily be transformed into an equivalent equation of motion. The influence functional of a thermal oscillator bath is a Gaussian functional of the double path representing the propagation of the reduced density matrix. This property can be used to construct the functional as the stochastic average of an exponentiated action functional corresponding to external driving by colored Gaussian c-number noise. A corresponding unraveling procedure [22, 28] of the functional employs two Gaussian processes whose correlation matrix reproduces both the real and imaginary parts of the reservoir correlation function (6). This translates into a time-local stochastic equation of motion for the reduced density called the SLN. This equation is formally exact and has proven to be useful for damped systems exposed to strong driving [24, 29] and energy transfer in systems coupled to reservoirs with structured spectral density [30].

If the spectral density is Ohmic, the imaginary part in L_r(t) reduces to a derivative of a delta function and therefore is time-local. Then, only the real part in (6) needs to be reconstructed by a non-Markovian noise term. The SLN then reduces to the so-called SLN for dissipation (SLED) [23]:

$$\begin{eqnarray*}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\rho }_{\xi }={ \mathcal L }{\rho }_{\xi } & = & -\,\displaystyle \frac{{\rm{i}}}{{\hslash }}[{H}_{{\rm{s}}},{\rho }_{\xi }]+\displaystyle \frac{{\rm{i}}}{{\hslash }}{\xi }_{1}(t)[{q}_{1},{\rho }_{\xi }]+\displaystyle \frac{{\rm{i}}}{{\hslash }}{\xi }_{N}(t)[{q}_{N},{\rho }_{\xi }]\\ & & -\displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \frac{{\gamma }_{1}}{2}[{q}_{1},\{{p}_{1},{\rho }_{\xi }\}]-\displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \frac{{\gamma }_{N}}{2}[{q}_{N},\{{p}_{N},{\rho }_{\xi }\}],\end{eqnarray*}$$

with the damping constants γ providing the coupling of the chain's endpoints to the baths, and the quantum noise ${\xi }_{r}(t)$, whose correlation function is given by the real part of ${L}_{r}(t-t^{\prime} )$,

$$\begin{eqnarray}&&\langle {\xi }_{r}(t){\xi }_{r^{\prime} }(t^{\prime} )\rangle ={\delta }_{r,r^{\prime} }\mathrm{Re}{L}_{r}(t-t^{\prime} ),\end{eqnarray} \tag{ 7 }$$

where two distinct reservoirs indexed with r and $r^{\prime}$ act independently on the system. We emphasize that, even at ${k}_{{\rm{B}}}T=0$, quantum noise is still present since the coth function in (6) does not vanish in this limit. Gardiner has identified a similar equation as an adjoint equation [31] of the quantum Langevin [32]. In the present context, we consider the Schrödinger picture the primary formalism of the dynamics; the adjoint dynamics introduced later in the present work will propagate time-dependent observables. The physical reduced density matrix ρ of the system is given by the expectation value of the stochastic density $\rho =\langle {\rho }_{\xi }{\rangle }_{\mathrm{stoch}}$.

2.3. Parameterization through cumulants

For a harmonic system, SLED dynamics lead to Gaussian states for long-enough times, both for individual samples ${\rho }_{\xi }$ and the physical density matrix even when the initial state is not Gaussian. We restrict ourselves to Gaussian states in the following and characterize them through first and second cumulants (expectation values and covariances) of positions and momenta.

Moments of observables are obtained from ${\rho }_{\xi }$ by a double average, the combination of a trace operation $\langle \cdot {\rangle }_{\mathrm{tr}}=\mathrm{tr}(\cdot {\rho }_{\xi })$ and an expectation value with respect to the noise statistics,

$$\begin{eqnarray} &&\langle\kern-2pt\langle A \rangle\kern-2pt\rangle = \langle\kern-2pt\langle A{\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}.\end{eqnarray} \tag{ 8 }$$

Since we are dealing with a linear system, the twice-averaged first moments show effective classical behavior, i.e., exponentially damped oscillations. We eliminate these by choosing their initial values to be zero.

For the covariance of two arbitrary operators A and B, the double average allows the transformation

$$\begin{eqnarray}\mathrm{Cov}(A,B) & = & \tfrac{1}{2} \langle\kern-2pt\langle {AB}+{BA} \rangle\kern-2pt\rangle - \langle\kern-2pt\langle A \rangle\kern-2pt\rangle \langle\kern-2pt\langle B \rangle\kern-2pt\rangle \\ & = & \tfrac{1}{2} \langle\kern-2pt\langle {AB}+{BA} \rangle\kern-2pt\rangle - \langle\kern-2pt\langle A{\rangle }_{\mathrm{tr}}\langle B{\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}\\ & & + \langle\kern-2pt\langle A{\rangle }_{\mathrm{tr}}\langle B{\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}- \langle\kern-2pt\langle A \rangle\kern-2pt\rangle \langle\kern-2pt\langle B \rangle\kern-2pt\rangle \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&=\langle {\mathrm{Cov}}_{\mathrm{tr}}(A,B){\rangle }_{\mathrm{stoch}}+{\mathrm{Cov}}_{\mathrm{stoch}}(\langle A{\rangle }_{\mathrm{tr}},\langle B{\rangle }_{\mathrm{tr}})\end{eqnarray} \tag{ 10 }$$

which effectively splits the covariance into two terms which we will call mean trace covariance and stochastic covariance. It will be advantageous to treat these separately: when mean trace covariance (mtr) and stochastic covariance (msc) are considered with A and B among all the coordinates of the operator-valued vector $\vec{\sigma }={({q}_{1},{p}_{1},\ldots ,{q}_{N},{p}_{N})}^{t}$, the covariance matrix ${\boldsymbol{\Sigma }}$ of its components ${\sigma }_{j}$ thus can be split as

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}={{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}+{{\boldsymbol{\Sigma }}}^{(\mathrm{msc})},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{jk}}^{(\mathrm{mtr})}=\langle {\mathrm{Cov}}_{\mathrm{tr}}({\sigma }_{j},{\sigma }_{k}){\rangle }_{\mathrm{stoch}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{jk}}^{(\mathrm{msc})}={\mathrm{Cov}}_{\mathrm{stoch}}(\langle {\sigma }_{j}{\rangle }_{\mathrm{tr}},\langle {\sigma }_{k}{\rangle }_{\mathrm{tr}}).\end{eqnarray} \tag{ 13 }$$

This split will allow the translation of the SLED dynamics into deterministic equations of motion for each of the two terms.

2.4. Time evolution of system-trace cumulants

As a first step toward equations of motion for ${{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}$ and ${{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}$, we consider the time dependence of observables which are only averaged through $\langle \cdot {\rangle }_{\mathrm{tr}}$, i.e., quantities which are still random variables in the probability space of the Gaussian noise $\xi (t)$.

Their time evolution can be obtained from the adjoint dynamics of observables associated with SLED. Quite generally, this 'Heisenberg picture' of a quantum master equation is a non-unitary equation of motion, governed by the adjoint generator ${{ \mathcal L }}^{\dagger }$ as described by Breuer and Petruccione [33]. It is not identical with the standard Heisenberg picture governed by the global Hamiltonian.

In the case of SLED, given by (7), the adjoint equation is a stochastic equation of motion. The random adjoint-propagated operators ${A}_{\xi }(t)$ follow

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{A}_{\xi } & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}[{H}_{{\rm{s}}},{A}_{\xi }]-\displaystyle \frac{{\rm{i}}}{{\hslash }}{\xi }_{1}(t)[{q}_{1},{A}_{\xi }]-\displaystyle \frac{{\rm{i}}}{{\hslash }}{\xi }_{N}(t)[{q}_{N},{A}_{\xi }]\\ & & +\displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \frac{{\gamma }_{1}}{2}\{{p}_{1},[{q}_{1},{A}_{\xi }]\}+\displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \frac{{\gamma }_{N}}{2}\{{p}_{N},[{q}_{N},{A}_{\xi }]\}.\end{eqnarray} \tag{ 14 }$$

When positions or momenta are inserted for A_ξ, this equation appears to be identical to the operator-valued quantum Langevin equation [31, 32] at first glance. The subtle difference between the two approaches is the meaning of ξ. In the case of the quantum Langevin equation, it is operator-valued; in our case, it is real-valued. This eliminates certain ordering problems, a feature which will be helpful when evaluating higher-order correlations. The price we pay for this lies in the fact that the operator algebra of canonical variables is lost in this propagation; the time evolution of products or functions of the canonical variables must be obtained separately by inserting them directly into (14).

It is helpful to note that the contributions in ${{ \mathcal L }}^{\dagger }$ which provide the noises ${\xi }_{r}(t)$ are linear in q_r while the others are all quadratic or bilinear in momentum and position. These contributions reduce to linear terms in the equations of motion for single coordinates, while the noise terms reduce to inhomogeneities ${\xi }_{r}(t)$. We therefore get a closed system of equations for the positions and momenta of the chain, which holds equally for first moments with respect to trace. These are shown in appendix A and can be summarized in the linear equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle \vec{\sigma }{\rangle }_{\mathrm{tr}}={\bf{M}}\langle \vec{\sigma }{\rangle }_{\mathrm{tr}}+\vec{\xi }(t),\end{eqnarray} \tag{ 15 }$$

where $\vec{\xi }={(0,{\xi }_{1},0,...,0,{\xi }_{N})}^{t}$ and where ${\bf{M}}$ is a $2N\times 2N$ matrix which is determined by the parameters of the system Hamiltonian.

The evolution of the second cumulants is accessible with the adjoint SLED in a similar manner. Inserting $A={\sigma }_{j}{\sigma }_{k}$ into (14) yields a right-hand side which is a linear combination of similar coordinate products and products of the form ${\sigma }_{l}{\xi }_{m}$. Combining the results for linear and quadratic terms yields the time derivatives of trace covariances, $\tfrac{{\rm{d}}}{{\rm{d}}t}\,{\mathrm{Cov}}_{\mathrm{tr}}({\sigma }_{j},{\sigma }_{k})=\tfrac{{\rm{d}}}{{\rm{d}}t}(\langle {\sigma }_{j}{\sigma }_{k}+{\sigma }_{k}{\sigma }_{j}{\rangle }_{\mathrm{tr}}/2-\langle {\sigma }_{j}{\rangle }_{\mathrm{tr}}\langle {\sigma }_{k}{\rangle }_{\mathrm{tr}})$. The covariance matrix ${{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}$ with elements ${\mathrm{Cov}}_{\mathrm{tr}}({\sigma }_{j},{\sigma }_{k})$ obeys the simple equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}={\bf{M}}{{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}+{{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}{{\bf{M}}}^{\dagger }.\end{eqnarray} \tag{ 16 }$$

This is independent of $\vec{\xi }(t)$, reflecting the fact that a spatially homogeneous force cannot induce squeezing in this model.

Due to the absence of noise in (16), ${{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}$ can immediately be identified with the mean trace covariance matrix ${{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}$, i.e., the mean trace covariance can be computed by propagating

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}={\bf{M}}{{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}+{{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}{{\bf{M}}}^{\dagger }.\end{eqnarray} \tag{ 17 }$$

The stochastic covariance part simplifies, under the initial conditions we have assumed, as

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{j,k}^{(\mathrm{msc})}= \langle\kern-2pt\langle {\sigma }_{j}{\rangle }_{\mathrm{tr}}\langle {\sigma }_{k}{\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}.\end{eqnarray} \tag{ 18 }$$

Inserting the formal solution of (15),

$$\begin{eqnarray}&&\langle \vec{\sigma }{\rangle }_{\mathrm{tr}}={\int }_{0}^{t}{\rm{d}}t^{\prime} {\bf{G}}(t-t^{\prime} )\vec{\xi }(t^{\prime} ),\end{eqnarray} \tag{ 19 }$$

and using Green's function, we obtain

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}={\int }_{0}^{t}{\rm{d}}t^{\prime} {\int }_{0}^{t}{\rm{d}}t^{\prime\prime} {\bf{G}}(t-t^{\prime} )\langle \vec{\xi }(t^{\prime} ){\vec{\xi }}^{t}(t^{\prime\prime} ){\rangle }_{\mathrm{stoch}}{{\bf{G}}}^{\dagger }(t-t^{\prime\prime} ),\end{eqnarray} \tag{ 20 }$$

which can be evaluated without sampling over noise. The correlation matrix $\langle \vec{\xi }(t^{\prime} ){\vec{\xi }}^{t}(t^{\prime\prime} ){\rangle }_{\mathrm{stoch}}$ can easily be computed and tabulated from (7) by numerical Fourier transform or summation over Matsubara frequencies. The only non-zero entries in this matrix are the second and last diagonal elements, which contain the autocorrelation functions of the two reservoirs. In the stationary limit, (20) arises from both the present method and the quantum Langevin approach [18, 34]. When also incorporating the contribution of ${{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}$, the present approach is exact at any timescale.

Instead of computing this double convolution integral through a normal mode analysis, we construct a formal dynamical system of modest size which contains the elements of ${{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}$ as dynamical variables. Performing the derivative with respect to time and introducing an auxiliary variable ${\rm{y}}$ yields a closed system of linear differential equations

$$\begin{eqnarray}&&{\dot{{\boldsymbol{\Sigma }}}}^{(\mathrm{msc})}(t)={{\rm{y}}}^{\dagger }(t)+{\rm{y}}(t)+{\bf{M}}{{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}(t)+{{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}(t){{\bf{M}}}^{\dagger },\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\dot{{\bf{y}}}(t)={\bf{G}}(t){\bf{L}}(t),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\dot{{\bf{G}}}(t)={\bf{MG}}(t),\end{eqnarray} \tag{ 23 }$$

with the matrix ${\bf{L}}(t)=\langle \vec{\xi }(t){\vec{\xi }}^{t}(0){\rangle }_{\mathrm{stoch}}$. As initial conditions we have ${\bf{G}}(0)={\mathbb{1}}$, and we choose throughout the whole paper ${{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}(0)=0$ as well as ${\bf{y}}(0)=0$ (no initial system–reservoir correlations). ${{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}(0)$ is determined by the unperturbed ground states of the chain's oscillators. The solutions from (16) and (21)–(23) enable one to exactly calculate the transient dynamics of ${\boldsymbol{\Sigma }}$, which provides all information about the state of the system. Since ${\bf{M}}$ is taken from a dissipative linear system, the propagation of these equations is numerically stable even for large chain length. For an efficient computation of steady-state behavior, we set the derivative with respect to time in (21) equal to zero and obtain a Lyapunov equation

$$\begin{eqnarray}&&{\bf{M}}{{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}+{{\boldsymbol{\Sigma }}}^{(\mathrm{msc})}{{\bf{M}}}^{\dagger }+{{\bf{y}}}^{\dagger }+{\bf{y}}=0\end{eqnarray} \tag{ 24 }$$

with an inhomogeneity that is determined by integrating (22) and (23) to times large enough that ${\bf{y}}$ reaches a constant steady-state value.

A physical interpretation of ${\bf{y}}$ is manifested by integration of (22)

$$\begin{eqnarray}&&{\bf{y}}(t)={\bf{y}}(0)+{\int }_{0}^{t}{\rm{d}}t^{\prime} {\bf{G}}(t-t^{\prime} )\langle \vec{\xi }(t){\vec{\xi }}^{t}(t^{\prime} ){\rangle }_{\mathrm{stoch}}.\end{eqnarray} \tag{ 25 }$$

By comparison with (19), this can be identified as the correlations of the bath fluctuations with the system's degrees of freedom given by first cumulants:

$$\begin{eqnarray}&&{\bf{y}}(t)= \langle\kern-2pt\langle \vec{\sigma }{\rangle }_{\mathrm{tr}}(t){\vec{\xi }}^{t}(t){\rangle }_{\mathrm{stoch}},\end{eqnarray} \tag{ 26 }$$

where ${\bf{y}}(0)=0$. With this interpretation, ${\xi }_{r}(t)$ appears as a 'stochastic substitute' for a reservoir operator. The time until ${\bf{y}}(t)$ is in equilibrium is provided by the time until the integrand in (25) decays to zero. Therefore, for moderate damping, the correlation time in ${\bf{L}}(t-t^{\prime} )$ is the pivotal time scale.

With the approach we presented in the previous section, we are able to determine the state of a harmonic quantum chain coupled to Ohmic reservoirs with arbitrary temperatures. This allows us to study nonequilibrium situations with finite heat fluxes causing an energy transport from a hot to a cold reservoir. The link topology of the chain shown in figure 1 suggests that we consider three cases. Two of them account for the coupling of the endpoints to the neighbor and hot or cold reservoirs, respectively, while the third covers the oscillators in the bulk which are coupled to next neighbors. By employing ${{ \mathcal L }}^{\dagger }$ on ${H}_{n}=\tfrac{{p}_{n}^{2}}{2m}$ $+\,\tfrac{1}{2}m{\omega }_{0}^{2}{q}_{n}^{2}$ $+\,\tfrac{\mu }{2}[{({q}_{n-1}-{q}_{n})}^{2}$ $+\,{({q}_{n}-{q}_{n+1})}^{2}]$, we obtain three adjoint dynamical equations for the energy operators:

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{H}_{1} & = & \displaystyle \frac{1}{m}{\xi }_{1}(t){p}_{1}-{\gamma }_{1}\displaystyle \frac{{p}_{1}^{2}}{m}+\displaystyle \frac{\mu }{m}{q}_{2}{p}_{1},\\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{H}_{n} & = & \displaystyle \frac{\mu }{m}{q}_{n-1}{p}_{n}+\displaystyle \frac{\mu }{m}{q}_{n+1}{p}_{n},\quad 2\leqslant n\leqslant N-1,\\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{H}_{N} & = & \displaystyle \frac{\mu }{m}{q}_{N-1}{p}_{N}+\displaystyle \frac{1}{m}{\xi }_{N}(t){p}_{N}-{\gamma }_{N}\displaystyle \frac{{p}_{N}^{2}}{m}.\end{eqnarray} \tag{ 27 }$$

For simplicity, the index ξ has been omitted from the adjoint-propagated observables. The assumption of locality of energy transfer between nearest neighbors allows the identification of individual terms on the rhs of $\tfrac{{\rm{d}}}{{\rm{d}}t}{H}_{n}={j}_{n-1,n}-{j}_{n,n+1}$, leading to

$$\begin{eqnarray}{j}_{n-1,n} & = & \displaystyle \frac{\mu }{m}{q}_{n-1}{p}_{n},\quad 2\leqslant n\leqslant N,\\ {j}_{n,n+1} & = & -\,\displaystyle \frac{\mu }{m}{q}_{n+1}{p}_{n},\quad 1\leqslant n\leqslant N-1.\end{eqnarray} \tag{ 28 }$$

To get from the energy flux operators to the first moments of the heat flux, we have to apply the stochastic average on the trace-averaged operators. Since we use initial conditions where all expectation values of phase-space operators are zero, the resulting currents are given by the elements of the covariance matrix:

$$\begin{eqnarray} &&\langle\kern-2pt\langle {j}_{n-1,n} \rangle\kern-2pt\rangle =\displaystyle \frac{\mu }{m} \langle\kern-2pt\langle {q}_{n-1}{p}_{n} \rangle\kern-2pt\rangle ,\quad 2\leqslant n\leqslant N,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&\langle\kern-2pt\langle {j}_{n,n+1} \rangle\kern-2pt\rangle =-\displaystyle \frac{\mu }{m} \langle\kern-2pt\langle {q}_{n+1}{p}_{n} \rangle\kern-2pt\rangle ,\quad 1\leqslant n\leqslant N-1,\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &&\langle\kern-2pt\langle {j}_{r,1} \rangle\kern-2pt\rangle =\displaystyle \frac{1}{m}\langle {\xi }_{1}(t)\langle {p}_{1}{\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}-{\gamma }_{1}\langle \hspace{-3.5pt}\langle \displaystyle \frac{{p}_{1}^{2}}{m}\rangle \hspace{-3.5pt}\rangle ,\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} &&\langle\kern-2pt\langle {j}_{r,N} \rangle\kern-2pt\rangle =-\displaystyle \frac{1}{m}\langle {\xi }_{N}(t)\langle {p}_{N}{\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}+{\gamma }_{N}\langle \hspace{-3.5pt}\langle \displaystyle \frac{{p}_{N}^{2}}{m}\rangle \hspace{-3.5pt}\rangle ,\end{eqnarray} \tag{ 32 }$$

where the second term in (31) and (32) is classical while the first term reflects correlations of the bath fluctuations with the system's degrees of freedom. According to (26), this term is an element of ${\bf{y}}$ which needs already to be computed for ${\boldsymbol{\Sigma }}$ and does not demand any further effort. In a steady state, the current in the bulk is constant over the chain and equal to the current between either the bath or the system. Therefore, for the analysis of steady-state currents, we write $\langle\kern-2pt\langle {j}_{r,1} \rangle\kern-2pt\rangle = \langle\kern-2pt\langle {j}_{n,n+1} \rangle\kern-2pt\rangle = \langle\kern-2pt\langle j \rangle\kern-2pt\rangle$.

The heat current provides us intuitive consistency checks such as a vanishing heat flux in the absence of a temperature gradient. This seemingly trivial result cannot be reproduced by naively constructed local Lindblad operators as was shown for harmonic oscillators and two level systems [10, 11]. In this respect, our approach delivers consistent results for various combinations of couplings and temperatures. In particular, low temperatures and unequal damping strengths of the endpoints do not lead to a violation of the second law.

An interesting feature of the steady-state flux that is shown in figure 2 is the interplay of the damping strengths γ and the couplings μ within the oscillators. For small damping, the plots show a linear increase of $\langle\kern-2pt\langle j \rangle\kern-2pt\rangle$ with the damping γ up to a maximum followed by an algebraic decay according to $\propto {\gamma }^{-1}$ for large damping strengths. This behavior was studied by Rieder et al in [35] for a system of classical harmonic oscillators without local frequencies ${\omega }_{0}$ and by Gaul and Büttner [18] for a similar quantum chain to the one used here. We could verify the observation of [18] that the position of the current's maximum increases linearly with the next-neighbor coupling μ for small couplings and falls back for large μ. The distinct sensitivity of the current on the damping will be relevant when we study heat-flux correlations under equilibrium and nonequilibrium conditions in the following sections.

Zoom In Zoom Out Reset image size

5.1. Stochastic calculations

The relative simplicity of the formal dynamics for the covariance matrix invites the study of space–time correlations of heat flux expressed through the covariance. As the fluxes themselves are given by a correlation of phase-space variables (see (29) and (30)), the heat-flux correlations represent an average over four operators which are simplified by an application of Wick's theorem reducing the Gaussian heat-flux correlations to products of phase-space correlations. The covariance of two currents ${j}_{n,n+1}(t)={j}_{n}(t)$ and ${j}_{l,l+1}(t^{\prime} )={j}_{l}(t^{\prime} )$ reads

$$\begin{eqnarray} &&\langle\kern-2pt\langle {j}_{n}(t){j}_{l}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c}= \langle\kern-2pt\langle {j}_{n}(t){j}_{l}(t^{\prime} ) \rangle\kern-2pt\rangle - \langle\kern-2pt\langle {j}_{n}(t) \rangle\kern-2pt\rangle \langle\kern-2pt\langle {j}_{l}(t^{\prime} ) \rangle\kern-2pt\rangle ,\end{eqnarray} \tag{ 33 }$$

with the last term being finite for nonequilibrium settings. Inserting (30) gives

$$\begin{eqnarray} \langle\kern-2pt\langle {j}_{n}(t){j}_{l}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c} & = & {\left(\displaystyle \frac{\mu }{m}\right)}^{2}[ \langle\kern-2pt\langle {q}_{n+1}(t){p}_{n}(t){q}_{l+1}(t^{\prime} ){p}_{l}(t^{\prime} ) \rangle\kern-2pt\rangle \\ & & - \langle\kern-2pt\langle {q}_{n+1}(t){p}_{n}(t) \rangle\kern-2pt\rangle \langle\kern-2pt\langle {q}_{l+1}(t^{\prime} ){p}_{l}(t^{\prime} ) \rangle\kern-2pt\rangle ],\end{eqnarray} \tag{ 34 }$$

where Wick's theorem decomposes the first term to products of phase-space variables

$$\begin{eqnarray} \langle\kern-2pt\langle {j}_{n}(t){j}_{l}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c} & = & {\left(\displaystyle \frac{\mu }{m}\right)}^{2}[ \langle\kern-2pt\langle {q}_{n+1}(t){q}_{l+1}(t^{\prime} ) \rangle\kern-2pt\rangle \langle\kern-2pt\langle {p}_{n}(t){p}_{l}(t^{\prime} ) \rangle\kern-2pt\rangle \\ & & + \langle\kern-2pt\langle {q}_{n+1}(t){p}_{l}(t^{\prime} ) \rangle\kern-2pt\rangle \langle\kern-2pt\langle {p}_{n}(t){q}_{l+1}(t^{\prime} ) \rangle\kern-2pt\rangle ].\end{eqnarray} \tag{ 35 }$$

Here we consider time-ordered products; operators appearing at equal time in (35) commute since they bear different site indices. For definiteness, we assume $t\gt t^{\prime}$. Since we wish to consider these correlations in the stationary state, the initial preparation will not be at time ${t}_{0}=0$, but at time ${t}_{0}=-\infty$ in the following. The time-ordered correlation matrix ${{\boldsymbol{\Sigma }}}_{\gt }(t-t^{\prime} )$ is then only a function of the time difference ${\rm{\Delta }}=t-t^{\prime}$; its value at ${\rm{\Delta }}=0$ is related to the results of section 3,

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\gt }(0)=\mathop{\mathrm{lim}}\limits_{\tau \to \infty }{\boldsymbol{\Sigma }}(\tau )+{{\boldsymbol{\Sigma }}}_{0}^{\prime\prime }=\mathop{\mathrm{lim}}\limits_{\tau \to \infty }{{\boldsymbol{\Sigma }}}^{(\mathrm{mtr})}(\tau )+{{\boldsymbol{\Sigma }}}_{0}^{\prime\prime }\end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray}{{\boldsymbol{\Sigma }}}_{0}^{\prime\prime }=\left(\begin{array}{ccccc}0 & \displaystyle \frac{{\rm{i}}}{2} & & & \\ -\,\displaystyle \frac{{\rm{i}}}{2} & 0 & \ddots & & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & 0 & \displaystyle \frac{{\rm{i}}}{2}\\ & & & -\,\displaystyle \frac{{\rm{i}}}{2} & 0\end{array}\right)\end{eqnarray} \tag{ 37 }$$

reflects the usual commutation relations, with a sign convention which determines a convention for the ordering of positions and momenta at equal time.

For arbitrary ${\rm{\Delta }}\gt 0$, the elements of ${{\boldsymbol{\Sigma }}}_{\gt }({\rm{\Delta }})$ are correlation functions of the form $\langle\kern-2pt\langle A{ \mathcal V }({\rm{\Delta }},0)B{ \mathcal V }(0,{t}_{0}) \rangle\kern-2pt\rangle$ with

$$\begin{eqnarray}&&{ \mathcal V }({t}_{2},{t}_{1})={{ \mathcal T }}_{\leftarrow }\exp \left({\int }_{{t}_{1}}^{{t}_{2}}{ \mathcal L }(s){\rm{d}}s\right).\end{eqnarray} \tag{ 38 }$$

The derivative of these correlation functions with respect to Δ can be analyzed in terms of the adjoint propagation and cast into matrix-valued equations of motion as in section 3; the time derivative transforms as

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial {\rm{\Delta }}} \langle\kern-2pt\langle A{ \mathcal V }({\rm{\Delta }},0)B{ \mathcal V }(0,{t}_{0}) \rangle\kern-2pt\rangle & = & \langle\kern-2pt\langle A{ \mathcal L }{ \mathcal V }({\rm{\Delta }},0)B{ \mathcal V }(0,{t}_{0}) \rangle\kern-2pt\rangle \\ & = & \langle\kern-2pt\langle {{ \mathcal L }}^{\dagger }A{ \mathcal V }({\rm{\Delta }},0)B{ \mathcal V }(0,{t}_{0}) \rangle\kern-2pt\rangle .\end{eqnarray} \tag{ 39 }$$

In appendix B, we show that this translates into a system of equations of motion which can be summarized in matrix form as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {\rm{\Delta }}}{{\boldsymbol{\Sigma }}}_{\gt }({\rm{\Delta }})={\bf{M}}{{\boldsymbol{\Sigma }}}_{\gt }({\rm{\Delta }})+{\bf{z}}({\rm{\Delta }})\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {\rm{\Delta }}}{\bf{z}}({\rm{\Delta }})=-{\bf{L}}({\rm{\Delta }})-{\bf{z}}({\rm{\Delta }}){{\bf{M}}}^{\dagger }\end{eqnarray} \tag{ 41 }$$

with initial conditions given by (36) and by ${\bf{z}}(0)={\mathrm{lim}}_{\tau \to \infty }{{\bf{y}}}^{\dagger }(\tau )$.

As with the results for the states shown in the previous sections, the correlations presented here are valid for arbitrary damping and temperatures. Therefore, our results overcome the restriction of weak damping that has to be respected when applying the quantum regression hypothesis which is commonly used in models for quantum optics [36]. It was pointed out that the hypothesis, which is based on a Born–Markov approximation, violates fundamental consistency criteria such as the fluctuation-dissipation theorem [19, 37, 38]. The formal reason for the validity of (40) and (41) lies in the fact that the propagation of the inhomogeneity ${\bf{z}}({\rm{\Delta }})$, as well as the propagations leading to its initial value, take into account fluctuations up to arbitrarily high order, which also reflects the system–reservoir correlations for strong damping and ensures the accordance with the fluctuation-dissipation theorem.

With the real and imaginary parts of the phase-space correlations given, one can construct the real and imaginary parts of the heat-flux correlations according to (35)

$$\begin{eqnarray}\mathrm{Re} \langle\kern-2pt\langle {j}_{n}(t){j}_{l}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c} & = & {\left(\displaystyle \frac{\mu }{m}\right)}^{2}[\mathrm{Re} \langle\kern-2pt\langle {q}_{n+1}\,{q}_{l+1} \rangle\kern-2pt\rangle \mathrm{Re} \langle\kern-2pt\langle {p}_{n}{p}_{l} \rangle\kern-2pt\rangle +\mathrm{Re} \langle\kern-2pt\langle {q}_{n+1}\,{p}_{l} \rangle\kern-2pt\rangle \mathrm{Re} \langle\kern-2pt\langle {p}_{n}\,{q}_{l+1} \rangle\kern-2pt\rangle \\ & & -\mathrm{Im} \langle\kern-2pt\langle {q}_{n+1}\,{q}_{l+1} \rangle\kern-2pt\rangle \mathrm{Im} \langle\kern-2pt\langle {p}_{n}\,{p}_{l} \rangle\kern-2pt\rangle -\mathrm{Im} \langle\kern-2pt\langle {q}_{n+1}\,{p}_{l} \rangle\kern-2pt\rangle \mathrm{Im} \langle\kern-2pt\langle {p}_{n}\,{q}_{l+1} \rangle\kern-2pt\rangle ],\end{eqnarray} \tag{ 42 }$$

where we skipped the dependence of the correlations on $(t,{\rm{\Delta }})$ for better legibility. For the imaginary part, one needs to compute

$$\begin{eqnarray}\mathrm{Im} \langle\kern-2pt\langle {j}_{n}(t){j}_{l}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c} & = & {\left(\displaystyle \frac{\mu }{m}\right)}^{2}[\mathrm{Re} \langle\kern-2pt\langle {q}_{n+1}\,{q}_{l+1} \rangle\kern-2pt\rangle \mathrm{Im} \langle\kern-2pt\langle {p}_{n}\,{p}_{l} \rangle\kern-2pt\rangle +\mathrm{Re} \langle\kern-2pt\langle {q}_{n+1}\,{p}_{l} \rangle\kern-2pt\rangle \mathrm{Im} \langle\kern-2pt\langle {p}_{n}\,{q}_{l+1} \rangle\kern-2pt\rangle \\ & & +\mathrm{Im} \langle\kern-2pt\langle {q}_{n+1}\,{q}_{l+1} \rangle\kern-2pt\rangle \mathrm{Re} \langle\kern-2pt\langle {p}_{n}\,{p}_{l} \rangle\kern-2pt\rangle +\mathrm{Im} \langle\kern-2pt\langle {q}_{n+1}\,{p}_{l} \rangle\kern-2pt\rangle \mathrm{Re} \langle\kern-2pt\langle {p}_{n}\,{q}_{l+1} \rangle\kern-2pt\rangle ].\end{eqnarray} \tag{ 43 }$$

Again, one should note that the pairs of imaginary parts in (42) provide contributions to the real parts of the heat-flux correlations.

5.2. Current fluctuations in a Gibbs state

Finally, we perform analytic calculations which resemble the limit of zero damping ($\gamma =0$) and allow us to verify the validity of our results obtained by the previous calculations of heat-flux correlations. Therefore, we use a Hamiltonian with fixed ends

$$\begin{eqnarray}&&{H}_{{\rm{s}}}=\displaystyle \sum _{n=1}^{N}\displaystyle \frac{{p}_{n}^{2}}{2m}+\displaystyle \frac{1}{2}m{\omega }_{0}^{2}{q}_{n}^{2}+\displaystyle \frac{\mu }{2}\displaystyle \sum _{n=1}^{N-1}{({q}_{n}-{q}_{n+1})}^{2}+\displaystyle \frac{\mu }{2}({q}_{1}^{2}+{q}_{N}^{2}),\end{eqnarray} \tag{ 44 }$$

which is an extension of (2) to boundary conditions that easily can be incorporated in the numerical SLED dynamics for comparison. The Hamiltonian in (44) can be interpreted as a chain with $N+2$ oscillators which yield ${q}_{0}={q}_{N+1}=0$. To achieve a normal mode representation of H_s, we introduce the transformations

$$\begin{eqnarray}&&{q}_{n}=\sqrt{\displaystyle \frac{2}{N+1}}\displaystyle \sum _{\nu =1}^{N}\sin ({k}_{\nu }n){{ \mathcal Q }}_{\nu }\qquad {p}_{n}=\sqrt{\displaystyle \frac{2}{N+1}}\displaystyle \sum _{\nu =1}^{N}\sin ({k}_{\nu }n){{ \mathcal P }}_{\nu }\end{eqnarray} \tag{ 45 }$$

with the wave vector of the νth normal mode ${k}_{\nu }=\pi \nu /(N+1),\ \nu ={{\mathbb{N}}}_{0}$. Applying these transformations and exploiting the fixed boundary conditions, we get

$$\begin{eqnarray}&&{H}_{{\rm{s}}}=\displaystyle \sum _{\nu }\left[\displaystyle \frac{{{ \mathcal P }}_{\nu }^{2}}{2m}+\displaystyle \frac{1}{2}m{\omega }^{2}({k}_{\nu }){{ \mathcal Q }}_{\nu }^{2}\right].\end{eqnarray} \tag{ 46 }$$

This Hamiltonian constitutes a sum of uncoupled oscillators, each with frequency $\omega ({k}_{\nu })\,=\,$ $\sqrt{{\omega }_{0}^{2}+(4\mu /m){\sin }^{2}({k}_{\nu }/2)}$. Applying the Hamilton equations gives us the standard differential equation of the harmonic oscillator for the νth mode and the classical dynamics of ${{ \mathcal Q }}_{\nu }({\rm{\Delta }})$ and ${{ \mathcal P }}_{\nu }({\rm{\Delta }})$, whose quantum nature is provided by the non-commutativity of the initial values ${{ \mathcal Q }}_{\nu }(0)$ and ${{ \mathcal P }}_{\nu }(0)$. Therefore, thermal averages over pairs of the initial values have to be calculated and put together with the classical dynamics plus the transformations from (45) to arrive at

$$\begin{eqnarray}\langle {q}_{n}({\rm{\Delta }}){q}_{l}(0){\rangle }_{\mathrm{tr}} & = & \displaystyle \frac{2}{N+1}\displaystyle \frac{1}{2m}\displaystyle \sum _{\nu }\sin ({k}_{\nu }n)\sin ({k}_{\nu }l)\\ & & \times \displaystyle \frac{1}{\omega ({k}_{\nu })}\left[\coth \left(\displaystyle \frac{\beta \omega ({k}_{\nu })}{2}\right)\cos (\omega ({k}_{\nu }){\rm{\Delta }})-{\rm{i}}\,\sin (\omega ({k}_{\nu }){\rm{\Delta }})\right],\end{eqnarray} \tag{ 47 }$$

where we only have to take the trace average of a canonical Gibbs state into account and can skip the stochastic average that occurs in the numerical method. From the correlation of positions, by considering ${p}_{n}({\rm{\Delta }})=m{\dot{q}}_{n}({\rm{\Delta }})$, one can calculate all other elements of the correlation matrix ${\boldsymbol{\Sigma }}({\rm{\Delta }})$ from (35) and the sought-after higher-order correlations.

Without the particular focus on correlations of heat currents, the present canonical calculation has been presented in [32] under allusion to the extendability to any higher-order correlation function with an even number of operators. In that paper and, for example, in [16], the thermodynamic Gibbs calculation is compared to dynamical models treated with Langevin equations. Multi-time correlations are particularly challenging in this framework, as the operator-valued quantum noise demands careful operator ordering.

In order to check the validity of the dynamical description based on the SLED and our numerical implementation, we compare its results with the ones from a thermal Gibbs state. We consider the Hamiltonian in (44) and compute the covariance of heat fluxes for varying time difference ${\rm{\Delta }}=t-t^{\prime}$ and lead index n. Figure 3 plots the real and imaginary parts of $\langle\kern-2pt\langle {j}_{n}(t){j}_{N/2}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c}$ in density plots (top) for a chain with N = 30 oscillators where the ends are weakly damped with $\gamma ={10}^{-4}$ and the reservoir temperatures are equal with ${k}_{{\rm{B}}}{T}_{r}=0$. For such a weak damping, we expect a good agreement with the numerical method based on the SLED and the results of the Gibbs calculations, which represent the limit of zero damping. Indeed, a quantitative agreement is shown in the lower plots of figure 3, where the covariance versus the index n is shown for two different values of Δ. For a small time difference ${\rm{\Delta }}=5$ (c), the results of SLED and Gibbs agree with high accuracy, as well as for a larger ${\rm{\Delta }}=50$ (d). In particular, the agreement of the two methods at the endpoints of the chain indicate that we are in a limit where damping is negligible over the considered time scales since one would expect that effects of damping first start to emerge at the oscillators close to the reservoirs.

Zoom In Zoom Out Reset image size

If larger damping like $\gamma =0.5$ is considered in the SLED approach, as done in figure 4, one observes deviations of the covariances resulting from the stochastic approach and the assumption of a canonical Gibbs state. The upper plots show the real part of the spatiotemporal correlations from analytical Gibbs calculations (a) and numerical simulations (b). While the correlations are undamped for the canonical ensemble, the dynamics from the SLED are weakly damped as apparent for time differences Δ, which are large enough to find finite correlations at the chain's endpoints. This effect of damping is more pronounced in plots (c) and (d), showing the covariance over the index n for two different times Δ. While the left plot which shows $\mathrm{Re} \langle\kern-2pt\langle {j}_{n}(t){j}_{N/2}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c}$ for ${\rm{\Delta }}=5$ reveals a good agreement for the whole chain, deviations of Gibbs and SLED become apparent in the vicinity of the endpoints if the time difference is larger, like ${\rm{\Delta }}=40$ in the right plot. With the previous validity checks, we show that our numerical approach reflects the limit of zero damping of the canonical Gibbs ensemble while deviations for finite damping occur, as to be expected, at the ends of the chain which are coupled to reservoirs. We have also found agreement between the present methods for small dampings and finite but equal temperatures at the attached reservoirs. Nevertheless, nonequilibrium situations caused by different temperatures remain withheld to the numerical SLED method.

Zoom In Zoom Out Reset image size

For the model with reservoirs only being attached at the endpoints of a sufficiently long chain, one has to admit that the damping of the reservoir only affects correlations in the vicinity of the end points. This might call into question the need of a numerical method valid for strong damping, as the Gibbs calculation presented here delivers correct results for the bulk. But a pitfall that arises particularly for large chains is the decreasing level spacing of the Hamiltonian with increasing number of modes. This induces rigid upper bounds for the damping to ensure that the reservoir-induced level broadening is much smaller than the level spacing, which is important for all formalisms where the system is supposed to be in a Gibbs state. Besides the damping strength, we also have the freedom to vary the temperatures leading to a nonequilibrium situation. This cannot be calculated by a canonical Gibbs state and provides us effects on the entire chain.

In figure 5, we study nonequilibrium effects achieved by the coupling of a system with N = 20 oscillators to two reservoirs with different temperatures ${k}_{{\rm{B}}}{T}_{1}=2.0$ and ${k}_{{\rm{B}}}{T}_{N}=0$. While plot (a) shows the real part of the covariance $\langle\kern-2pt\langle {j}_{n}(t){j}_{10}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c}$ for very weak damping $\gamma ={10}^{-3}$, panel (b) shows the covariance achieved with $\gamma =0.3$ and (d) shows $\mathrm{Re} \langle\kern-2pt\langle {j}_{n}(t){j}_{10}(t^{\prime} ){ \rangle\kern-2pt\rangle }_{c}$ for $\gamma =1.5$. When comparing the covariances plotted over n, as is done for ${\rm{\Delta }}=20$ in (c), one immediately sees that the impact of the temperature gradient applied to the system is most distinct for $\gamma =0.3$. This is the damping γ providing the best 'match' with the couplings within the chain μ, and corresponds to the vicinity of the maximum of the steady-state heat flux $\langle\kern-2pt\langle j \rangle\kern-2pt\rangle$ shown in figure 2. Even the steady-state flux and the covariances of the heat fluxes are different observables and both are sensitive with respect to the coupling strengths to the heat baths. This shows that a method with validity beyond the weak coupling regime is crucial in order to study nonequilibrium effects on heat-flux correlations.

Zoom In Zoom Out Reset image size

The ability to study heat-flux correlations in nonequilibrium qualifies our method to be used for the quantification of transport in two respects. In order to decide if a system shows ballistic, diffusive, or localized behavior, one can consider a nonequilibrium situation and calculate the heat flux versus the chain length. This reveals an exponential or algebraic dependence and allows us to extract localization length scales. Second, one can use second-order correlations in thermal equilibrium, as done for anharmonic classical systems, in order to classify transport properties [39–41]. With the formalism presented here, we provide a generalization to the quantum regime and to situations out of equilibrium. According to [39–41] of specific relevance are momentum and energy correlations, with momentum correlations being inherent to fluctuations in the heat flux (see (42) and (43)), while energy fluctuations can be derived in a similar manner as the correlations of heat we targeted here. The impact of out-of-equilibrium situations on transport parameters such as localization length and diffusion exponent obtained in equilibrium will be studied in future work.

Although the damped harmonic oscillator has been studied extensively, a correct description of system–reservoir interactions (particularly for strong damping) is challenging but essential when computing thermodynamic quantities such as heat fluxes. Based on the stochastic Liouville–von Neumann equation, we presented an efficient computational scheme for damped systems which proved to be reliable over the full range of damping strengths and temperatures. We exploited this flexibility to also study heat currents for stronger coupling to thermal reservoirs and observed a similar behavior as is found in approaches based on Langevin equations—namely, a distinct maximum for the currents when the reservoir coupling strength is tuned; its position depends on the intra-chain couplings.

As a main result, the scheme also allows us to gain spatiotemporal correlations of heat fluxes. For weak damping, excellent agreement is found with analytical calculations in a canonical Gibbs ensemble. However, beyond this domain, substantial deviations occur. The parameters for which we found the maximum in the currents led to the most distinct nonequilibrium effects in the fluctuations. The most distinct nonequilibrium effects in the correlations appear in ranges of parameter space where mean currents exhibit maxima.

The present approach can now easily be extended to disordered systems. The computational efficiency allows one to study heat fluxes in large chains with randomization of any system parameter. The spatiotemporal correlation patterns might enable one to investigate localization and diffusion solely in the time domain, which would represent an alternative to methods based on a transformation to normal modes. Moreover, an extension to driven systems where the external driving can either couple linearly or quadratically to the system's position is possible. The latter represents the important case of parametric driving and is of particular relevance for nano heat engines. For example, the spatial motion of an ion in a linear Paul trap was induced in [42, 43] by a periodic narrowing and widening of the system's ground-state frequency. In these situations, ${\bf{M}}$ from (15) becomes time dependent and the Green's function in (23) does not remain time translational invariant. This demands that we decompose the Green's function into a forward and backward propagator ${\bf{G}}(t,t^{\prime} )={{\bf{G}}}_{{\rm{f}}}(t,0){{\bf{G}}}_{{\rm{b}}}(0,t^{\prime} )$ with the equations of motion given by ${\dot{{\bf{G}}}}_{{\rm{f}}}(t,0)={\bf{M}}(t){{\bf{G}}}_{{\rm{f}}}(t,0)$ and ${\dot{{\bf{G}}}}_{{\rm{b}}}(0,t^{\prime} )=-{{\bf{G}}}_{{\rm{b}}}(0,t^{\prime} ){\bf{M}}(t^{\prime} )$ for the forward and backward propagation, respectively. This formalism was used to study entanglement generation through local driving in a bipartite system in analogy to [29], and we consider it to open ways to study different settings for quantum heat engines beyond the limitations of weak couplings and adiabatic driving modes.

We thank Udo Seifert for valuable discussions. This work was supported by Deutsche Forschungsgemeinschaft through grant AN336/6-1.

Based on (2) and (14), we calculate the equations of motion for the first and second cumulants. We assume a general  model, where each oscillator n can be coupled to a reservoir leading to a nonzero damping ${\gamma }_{n}$. For the site indexed with n, we get (${a}_{n}=m{\omega }_{0}^{2}+{\mu }_{n}+{\mu }_{n-1}$)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {q}_{n}{\rangle }_{c}=\displaystyle \frac{\langle {p}_{n}{\rangle }_{c}}{m},\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {p}_{n}{\rangle }_{c} & = & +\,{\mu }_{n-1}\langle {q}_{n-1}{\rangle }_{c}-{a}_{n}\langle {q}_{n}{\rangle }_{c}\\ & & -{\gamma }_{n}\langle {p}_{n}{\rangle }_{c}+{\mu }_{n}\langle {q}_{n+1}{\rangle }_{c}+{\xi }_{n}(t),\end{eqnarray} \tag{ A.2 }$$

where ${\gamma }_{n}$ and ${\xi }_{n}(t)$ are only non-zero if $n=1,N$ for the system displayed in figure 1. The couplings are indexed since the topology of a chain where ${\mu }_{n}$ connects the site n with its neighbor $n+1$ is guaranteed by the boundary condition ${\mu }_{0}={\mu }_{N}=0$. For N modes, we have a matrix equation $\tfrac{{\rm{d}}}{{\rm{d}}t}\langle \vec{\sigma }{\rangle }_{c}={\bf{M}}\langle \vec{\sigma }{\rangle }_{c}+\vec{\xi }(t)$, with ${\bf{M}}$ being a $2N\times 2N$ matrix

$$\begin{eqnarray}{\bf{M}}=\left(\begin{array}{cccccc}0 & 1/m & 0 & ... & ... & 0\\ -\,{a}_{n} & -\,{\gamma }_{1} & {\mu }_{1} & ... & ... & 0\\ \vdots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ 0 & ... & 0 & 0 & 0 & 1/m\\ 0 & ... & {\mu }_{N-1} & 0 & -\,{a}_{n} & -\,{\gamma }_{N}\end{array}\right)\end{eqnarray} \tag{ A.3 }$$

where the coefficients of the nth phase-space variables form a 2 × 2 block structure and the couplings appear in the off-diagonals.

The second cumulants are equivalently calculated to the first and read for the site with index n:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {q}_{n}^{2}{\rangle }_{c}=\displaystyle \frac{2}{m}\langle {q}_{n}{p}_{n}{\rangle }_{c},\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {p}_{n}^{2}{\rangle }_{c} & = & -\,2{\gamma }_{n}\langle {p}_{n}^{2}{\rangle }_{c}-2{a}_{n}\langle {q}_{n}{p}_{n}{\rangle }_{c}\\ & & +2{\mu }_{n-1}\langle {q}_{n-1}{p}_{n}{\rangle }_{c}+2{\mu }_{n}\langle {q}_{n+1}{p}_{n}{\rangle }_{c},\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {q}_{n}{q}_{l}{\rangle }_{c}=\displaystyle \frac{\langle {q}_{l}{p}_{n}{\rangle }_{c}}{m}+\displaystyle \frac{\langle {q}_{n}{p}_{l}{\rangle }_{c}}{m},\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {p}_{n}{p}_{l}{\rangle }_{c} & = & -\,{a}_{n}\langle {q}_{n}{p}_{l}{\rangle }_{c}-{a}_{n}\langle {q}_{l}{p}_{n}{\rangle }_{c}-({\gamma }_{n}+{\gamma }_{l})\langle {p}_{n}{p}_{l}{\rangle }_{c}\\ & & +{\mu }_{n}\langle {q}_{n+1}{p}_{l}{\rangle }_{c}+{\mu }_{n-1}\langle {q}_{n-1}{p}_{l}{\rangle }_{c}\\ & & +{\mu }_{l}\langle {q}_{l+1}{p}_{n}{\rangle }_{c}+{\mu }_{l-1}\langle {q}_{l-1}{p}_{n}{\rangle }_{c}\end{eqnarray} \tag{ A.7 }$$

,

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {q}_{n}{p}_{l}{\rangle }_{c} & = & -\,{a}_{n}\langle {q}_{n}{q}_{l}{\rangle }_{c}+\displaystyle \frac{\langle {p}_{n}{p}_{l}{\rangle }_{c}}{m}-{\gamma }_{l}\langle {q}_{n}{p}_{l}{\rangle }_{c}\\ & & +{\mu }_{l-1}\langle {q}_{l-1}{q}_{n}{\rangle }_{c}+{\mu }_{l}\langle {q}_{n}{q}_{l+1}{\rangle }_{c}.\end{eqnarray} \tag{ A.8 }$$

The only entries in ${\bf{M}}$ which do not represent unitary evolution are the non-zero diagonal entries $-{\gamma }_{l}$. This allows us to rewite the preceding equations in the compact form of (16),

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}={\bf{M}}{{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}+{{\boldsymbol{\Sigma }}}^{(\mathrm{tr})}{{\bf{M}}}^{\dagger }.\end{eqnarray*}$$

We derive the dynamics of the elements in ${{\boldsymbol{\Sigma }}}_{\gt }({\rm{\Delta }})$ by considering (39) for all combinations of phase-space variables. The result is in close analogy to the equations of appendix A. However, it is simpler than the propagation leading to the covariance matrix since the adjoint Liouvillian ${{ \mathcal L }}^{\dagger }$ is applied only to one of the operators in the product. Abbreviating

$$\begin{eqnarray} &&\langle\kern-2pt\langle A({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle = \langle\kern-2pt\langle A{ \mathcal V }({\rm{\Delta }},0)B{ \mathcal V }(0,{t}_{0}) \rangle\kern-2pt\rangle \end{eqnarray} \tag{ B.1 }$$

and considering $A\in \{{q}_{n},{p}_{n}\}$, we find from (39) that

$$\begin{eqnarray}&&{\partial }_{{\rm{\Delta }}} \langle\kern-2pt\langle {q}_{n}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle =\displaystyle \frac{1}{m} \langle\kern-2pt\langle {p}_{n}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle ,\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}{\partial }_{{\rm{\Delta }}} \langle\kern-2pt\langle {p}_{n}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle & = & {\mu }_{n-1} \langle\kern-2pt\langle {q}_{n-1}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle -{a}_{n} \langle\kern-2pt\langle {q}_{n}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle \\ & & -{\gamma }_{n} \langle\kern-2pt\langle {p}_{n}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle +{\mu }_{n} \langle\kern-2pt\langle {q}_{n+1}({\rm{\Delta }})B(0) \rangle\kern-2pt\rangle +\langle {\xi }_{n}({\rm{\Delta }})\langle B(0){\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}.\end{eqnarray} \tag{ B.3 }$$

The last term in (B.3) represents elements of ${\bf{z}}({\rm{\Delta }})$. It requires further treatment in order to arrive at closed equations of motion. Observing that

$$\begin{eqnarray}&&\langle {\xi }_{n}({\rm{\Delta }})\langle B(0){\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}=\langle {\xi }_{n}(0)\langle B(-{\rm{\Delta }}){\rangle }_{\mathrm{tr}}{\rangle }_{\mathrm{stoch}}\end{eqnarray} \tag{ B.4 }$$

and applying the adjoint (14) to $B(-{\rm{\Delta }})$ leads to (41).