Rectification of heat currents across nonlinear quantum chains: A versatile approach beyond weak thermal contact

Within the emerging field of quantum thermodynamics the issues of heat transfer and heat rectification are basic ingredients for the understanding and design of heat engines or refrigerators at nanoscales. Here, a consistent and versatile approach for mesoscopic devices operating with continuous degrees of freedom is developed valid from low up to strong system-reservoir couplings and over the whole temperature range. It allows to cover weak to moderate nonlinearities and is applicable to various scenarios including the presence of disorder and external time-dependent fields. As a particular application coherent one-dimensional chains of anharmonic oscillators terminated by thermal reservoirs are analyzed with particular focus on rectification. The efficiency of the method opens a door to treat also rather long chains and extensions to higher dimensions and geometries.


Introduction
Thermodynamics emerged as a theory to describe the properties of systems consisting of macroscopically many degrees of freedom. While it served from the very beginning as a practical tool to quantify the performance of heat engines, it also initiated substantial efforts in fundamental research to pave the road for fields like statistical and nonequilibrium physics. In recent years, triggered by the on-going progress in miniaturizing devices down to the nanoscale, the crucial question if, and if yes, to what extent macroscopic thermodynamics is influenced by quantum mechanics has led to a flurry of literature and the appearance of the new field of quantum thermodynamics.
Of particular relevance as a pre-requisite for the design of actual devices is the understanding of heat transfer and heat rectification. In macroscopic structures heat transport is typically characterized by normal diffusion such that the heat conductivity is independent of the size of the probe (Fourier's law). The requirements for microscopic models that support this type of normal heat flow has been subject to controversial debates [1,2,3,4,5,6,7,8] with dimensionality, disorder, and nonlinearities as potential ingredients.
Thermal rectification describes a difference of the heat fluxes' absolute values dependent on the direction of the temperature gradient the system is exposed to. It has turned out that nonlinearities in combination with spatial symmetry breaking are pivotal conditions for the occurrence of rectification [9,10]. Rectification can be used in thermal diodes or heat valves which have been proposed in classical [11,12,13,14,15,16,17,18] as well as in quantum systems [9,10,19,20,21,22,23]. As an alternative to nonlinearities, a coupling of the system to self-consistent reservoirs which guide the constituents of a chain to local equilibrium have been studied [24,25,26,27,28,29,30]. Experimentally, thermal rectification has been realized for example in carbon nanotubes [31], trapped ions in optical lattices [32] and, recently, in superconducting quantum bits [33]. Particularly superconducting circuits allow for a well-controlled modulation of nonlinearities [33,34,35] and may thus serve as promising candidates for the realization of heat engines operating in the deep quantum regime. In fact, situations where heat transfer is beneficial for the performance of devices have been discussed for the fast initialization of quantum bits by cooling [36] as well as for properties of heat valves, thermal memories, switches, and transistors [23,37,38,39].
From the theory point of view, the description of quantum heat transfer is a challenging issue. As part of the broad field of open quantum systems, it became clear in the last years that the underlying processes are extremely sensitive to a consistent treatment of the coupling between system and thermal reservoirs. Conventional weak coupling approaches (master equations) based on local dissipators may even violate the second law [40,41,42]. Moreover, due to current experimental activities a systematic approach for quantum heat transfer valid from the weak up to the strong coupling regime and down to very low temperatures is urgently required. The goal of this work is to present such a formulation that, in addition, stands out for its efficiency. While here for specific applications we focus on the experimentally relevant case of one-dimensional chains of anharmonic oscillators, it thus allows for generalizations to higher dimensions, more complex geometries, external driving, or set-ups such as heat valves and heat engines.
This approach is based on a stochastic description of non-equilibrium quantum dynamics (stochastic Liouville-von Neumann equation or SLN) which is an exact representation of the Feynman-Vernon path integral formulation in terms of a stochastic equation of motion for the unravelled reduced density operator [43,44,45]. For reservoirs with nearly ohmic spectral densities this SLN can further be reduced to contain only one real-valued noise field (so-called SLED) instead of two complex-valued ones [46]. It has then proven to be particularly powerful to describe systems with continuous degrees of freedom and in presence of external time-dependent driving [47,48]. However, the main challenge to practically implement the SLED is the degrading signal to noise ratio for increasing simulation time. Here, we circumvent this problem by formulating the quantum dynamics in terms of hierarchies of cumulants which are truncated properly. As we explicitly demonstrate in comparison with benchmark results from full SLED simulations, this treatment provides a very versatile tool to analyze quantum heat transfer in steady state across single or chains of anharmonic oscillators with weak to moderate anharmonicities. It thus allows to cover within one formulation broad domains in parameters that are not accessible by perturbative approaches [19,20]. Other alternative schemes are either limited to closed systems [49,50,51,52,53] or to very short structures interfacing the reservoirs [21].
The paper is organized as follows: In Section 2, we introduce the SLN which represents the basis for the cumulant truncation scheme presented in Section 3. First applications and comparison with benchmark resulrs are discussed in Section 4. The main results for heat rectification are then presented in Sections 5 and 6 including an analysis of the physical mechanism that determine the occurrence of rectification. The impact of disorder is considered in Section 7 before a summary and outlook is given in 8.

Non-perturbative reduced Dynamics
The description of heat transfer is a delicate issue. It is known that standard weak coupling approaches for the reduced dynamics of an aggregate terminated by two thermal reservoirs such as master equations suffer from fundamental inconsistencies [41,54]. In addition they are restricted to weak aggregate reservoir couplings and elevated temperatures. Hence, we here develop an approach which is based on a formally exact formulation of open quantum dynamics derived within a system + bath model. The total Hamiltonian consists of a system Hamiltonian H s , a reservoir Hamiltonian H R (for the moment we consider a single reservoir), and a system reservoir coupling H I = −qX, where the latter captures a bilinear coupling of a system's coordinate q and a collective bath coordinate X. Due to the macroscopically many degrees of freedom of the reser-voir, the fluctuations of the latter in thermal equilibrium can assumed to be Gaussian.
Then, as shown previously, the formally exact path integral representation of the dynamics of the reduced density ρ(t) = Tr R {W (t)} with W (t) being the time evolved density operator in full Hilbert space, has an equivalent representation in terms of a stochastic Liouville-von Neumann equation (SLED) [46], i.e., Here, ξ(t) denotes a c-valued noisy driving force whose auto-correlation reconstructs the real part of the bath quantum correlation L(t − t ) where with inverse temperature β = 1/(k B T ) and the spectral density J(ω). The latter we assume to be ohmic with a large cut-off frequency ω c and coupling constant γ acting on a system with mass m, i.e., In this regime (large cut-off), the imaginary part of the correlation function ImL(t) collapses to a δ-function and is accounted for by the γ-dependent in (1). We note in passing that Gardiner has identified this equation as the adjoint of a quantum Langevin equation [55,56]. The random density ρ ξ propagated according to (1) by itself lacks of a physical interpretation, while only the expectation value of ρ = ρ ξ ξ represents the physical reduced density of the system. The SLED is particularly suited to capture the dynamics of open quantum systems with continuous degrees of freedom and has been explicitly applied so previously in various contexts [42,57]. In the sequel, we follow a somewhat different route though by exploiting that the multiplicative noise ξ in the SLED turns into additive noise if an adjoint equation governed by L † [58] for the dynamics of Heisenberg operators A is used, i.e., d dt Expectation values are then obtained from a quantum mechanical average · tr = tr(·ρ ξ ) and a subsequent noise average · ξ , i.e., Now, to avoid an explicit noise sampling with its inherent degradation of the signal to noise ratio for longer simulation times, we do not explicitly work with (5) but rather derive from it sets of equations of motion for expectation values [47,48]. This way, one arrives at a very efficient scheme to construct the open system dynamics for ensembles of anharmonic oscillators also in regimes which are not accessible by perturbative methods for open quantum systems. A central element of this formulation is the covariance of two operators A and B which can be transformed to [42,59] Cov(A, and thus provides a separation into two covariances, one with respect to the trace average and one with respect to the noise average. Hence, choosing A and B as elements of the operator valued vector σ = (q, p) t carrying position and momentum operators, the corresponding covariance matrix Σ takes the form Technically, this decomposition requires to carefully distinguish between trace-moments and trace-cumulants. Namely, the elements of Σ contain noise expectation values of trace-moments, while Σ (mtr) contain noise expectation values of trace-covariances. For general operators A, B we thus introduce the compact notation for noise expectation values of trace-covariances Cov tr (A, B) ξ : with AB tr,c = AB tr − A tr B tr . We also emphasize the equality of the first moments and cumulants A tr,c = A tr , AB c = AB tr,c ξ,c .
With these tools at hand, we can now proceed to develop the approach for nonlinear oscillators in detail.

Cumulant formulation of open quantum dynamics for nonlinear oscillators
We will start to consider a single anharmonic oscillator to arrive at a formulation that can then easily be generalized to chains of oscillators. As a paradigmatic model for a nonlinear oscillator of mass m we choose with fundamental frequency ω and anharmonicity parameter κ, see figure. 1, which can be both positive (stiffer mode) or negative (softer mode). Our main interest is in weakly anharmonic systems, i.e., even for negative κ, where the potential is not bounded from below, resonances are long-lived and may be treated approximately as eigenstates of the Hamiltonian. This is valid if thermalization through external reservoirs is fast compared to the resonance lifetime.
While for a purely linear system (κ = 0) a formulation in terms of cumulants leads to closed equations of motion for the first and second order cumulants [47,48], this is no longer the case for nonlinear oscillators. The challenge is thus to implement a systematic procedure for the truncation of higher order cumulants which on the one hand provides an accurate description for weak to moderate anharmonicity parameters and on the other hand leads to an efficient numerical scheme also for ensembles of those oscillators. We emphasize that we are primarily interested in quantum heat transfer in steady state situations and not in the full wealth of nonlinear dissipative quantum dynamics. Accordingly, in a nutshell, a formulation which satisfies both criteria factorizes all higher than second order moments (Wick theorem) and allows to capture anharmonicities effectively by state dependent frequencies which must be determined self-consistently. This perturbative treatment of nonlinearities for quantum oscillators has some similarities in common with the one for classical systems [60] but turns out to be much more involved due to the highly complex equation of motion (5) as we will show now.
For this purpose, for systems, where the relevant dynamics occurs around a single potential minimum, it is convenient to set initial values of operators A equal to zero and use henceforth. Then, starting from (5) one obtains two coupled equations for position and momentum according to d dt q tr,c = 1 m p tr,c d dt p tr,c = −mω 2 q tr,c − mκ q 3 tr − γ p tr,c + ξ(t) .
Here, the third trace-moment can be expressed as a linear combination of products of cumulants q 3 tr = q 3 tr,c + 3 q 2 tr,c q tr,c + q 3 tr,c .
This kind of transformation represents a systematic separation and summation over all possible subsets of trace-averaged products, for details see [61].

a) Truncation of trace-and noise-cumulants
A straightforward approach to treat higher than second order moments is to assume the approximate validity of Wick's theorem and neglect all higher than second order cumulants so that (17) reduces to q 3 tr ≈ 3 q 2 tr,c q tr,c + q 3 tr,c .
This procedure does not immediately lead to a closure of (16) although with parametric driving through q 2 tr,c , which couples to the equations for the second cumulants which we analyze later.
As shown for the linear system, the covariance matrix consists of two parts, where the elements of one part is given by the noise averaged product A tr,c B tr,c ξ [59]. To handle these terms, we derive the equations of motion for these products with L † and turn the noise-moments to noise-cumulants which is trivial for all linear contributions since A tr,c ξ,c = A = 0 and also for the system-bath correlation which contains the quantum noise whose expectation value is zero. The resulting equations of motion then read d dt q tr,c q tr,c ξ,c = 2 m q tr,c p tr,c ξ,c d dt p tr,c p tr,c ξ,c = − 2mω 2 q tr,c p tr,c ξ,c − 2mκ[3 q 2 tr,c q tr,c p tr,c ξ + q 3 tr,c p tr,c ξ ] − 2γ p tr,c p tr,c ξ,c + 2 ξ(t) p tr ξ,c d dt q tr,c p tr,c ξ,c = 1 m p tr,c p tr,c ξ,c − mω 2 q tr,c q tr,c ξ,c − mκ[3 q 2 tr,c q 2 tr,c ξ + q 4 tr,c ξ ] − γ q tr,c p tr,c ξ,c + ξ(t) q tr ξ,c .
In equivalence to the equations of the trace-cumulants, here the contributions which can not be turned immediately from moments to cumulants are the ones which enter from the nonlinearity of the system. The second and third equations in (19) contain higher order moments with respect to the noise average. The trace-cumulants in these moments constitute c-numbers, and hence, the order of the higher order noise-moments is determined by the sum of the exponents outside of the trace averages . tr,c . Hence, q 2 tr,c A tr,c B tr,c ξ (A, B=q, p) constitutes a third order moment with respect to the noise and q 3 tr,c A tr,c ξ a fourth order moment. These moments as linear combinations of cumulants are Setting the third-and fourth-order cumulant to zero gives and where the equalities in (22) and (23) follow directly from (13). With (22), products of the noise expectation values of the first and second trace-cumulants enter. The dynamics of the second trace-cumulants is shown in Appendix A, where an analysis of the steadystate shows that the fluctuations induced by the coupling to the first trace-moments are exponentially suppressed and, therefore, this second trace-cumulants have a stable fixed point at zero. This leads to a decoupling of the equations for the first-and second trace-cumulants and a reduction of the covariance matrix to Σ = Σ (msc) . This allows us to simplify the notation by using A tr B tr ξ = AB which is valid for steady-states. The elements of Σ are then determined by a system of differential equations d dt which represent the steady-state by algebraic equations if the left hand sides are set to zero and the system-bath correlations ξ(t) σ j tr ξ are in their respective steady-state value. We introduced the effective frequencỹ which shows a similarity of our formalism to the treatment of classical anharmonic oscillators where amplitude dependent frequencies are well known [60]. In the quantum model considered here, the expectation value of the amplitude is zero ( q = 0) but the width of the state q 2 enters into the effective frequency. The noise average leads to a temperature dependency of q 2 and therefore to an effective frequencyω which depends on the temperature of the heat bath. A direct calculation of the dynamics of the full covariance matrix Σ reveals a dependence on higher order moments which are given by q 4 and q 3 p+pq 3 2 . The presented formalism represents a transformation of the moments contained in the covariance matrix into cumulants and an expansion with subsequent truncation of higher order cumulants. This procedure is equivalent to an approximation of the higher order moments given by The first equation is obviously the content of Wick's theorem, just like the second equation if the system is in steady-state where qp+pq 2 = 0 holds. Nevertheless, the truncation scheme we present here provides a systematic approach that accounts for both, the quantum and the thermal fluctuations. Therefore it represents an extension of previous schemes based on truncations of higher order moments or cumulants for closed systems [49,50,51,52,53].

b) System-bath correlation function
For a compact treatment of the system-bath correlations and a subsequent extension to multipartite systems we treat the system as a linear one with effective frequency and use a matrix notation of the derived steady-state equations (24) via a Lyapunov equation as in previous work [54,59]: The matrices have dimension 2N × 2N with N being the number of oscillators. M contains details of the model and the damping and reads for one oscillator with the effective frequencyω 2 = ω 2 + 3κ q 2 which itself depends on an element of the covariance matrix Σ. Therefore, (27) is solved self-consistently with the solution for the linear system (κ = 0) as an initial guess for M and the system-bath correlations. These are shifted to the matrix y which reads In terms of a tensor product of the vector carrying the phase space variables σ = (q, p) t and a vector with the noise ξ = (0, ξ) t , the correlation matrix is If the effective frequency is supposed to be fixed, the equations (24) can considered as a linear set of equations with additional driving ξ(t) which are formally solved by the Green's function: Figure 2. Schematic diagram of a chain of anharmonic oscillators which is terminated by thermal reservoirs. The dotted lines illustrate the on-site potentials for κ > 0 (blue) and κ < 0 (red), while the solid black line represents the harmonic case with κ = 0. The width of the potential determines the quantum state (green) with its effective frequencyω n .
The noise expectation value of the product of σ tr with the noise vector ξ t (t) gives a form of y where all noises are shifted to the bath auto-correlation function: The elements of the matrix L (t − t ) = ξ(t) ξ t (t ) ξ are given by (3) which completes the presented formalism to a deterministic description based solely on model parameters. For consistent results with (27) one has to integrate to sufficiently large times until a constant steady-state value of y is reached. By usage of the time-translational symmetry of the system, this integral can be derived with respect to time and written as a system of two coupled differential equationṡ which allows a simple and computationally efficient implementation. Treating both, the system-bath correlation and the dissipative parts as a linear system with effective frequencyω provides a consistent formalism for the nonlinear system.

c) Generalization to chains of oscillators
In order to step forward to chains of anharmonic oscillators, we choose a quadratic coupling between adjacent entities, i.e., and, for the sake of transparency, work with chains with homogeneous mass m, anharmonicitiy κ, and inter-oscillator coupling µ. The generalization to inhomogeneous structures is straighforward (see Sec. 7). Note that this Hamiltonian differs from the paradigmatic Fermi-Pasta-Ulam model [4] in that there nonlinearities appear in the inter-oscillator couplings while here they feature on-site pinning type of potentials. Now, to study quantum heat transfer through this extended structure, the formulation for a single oscillator is easily generalized along the lines presented in [42,59]. The chain is terminated at its left (l) and right (r) end by two independent reservoirs, see figure. 2, with bath correlations given by where the L ν (t) are specified in (3) and (4) with J(ω) → J ν (ω) according to γ → γ ν and β → β ν = 1/k B T ν . For heat transfer one considers situations with T l = T r corresponding to a 'cold' and a 'hot' reservoir with temperatures T h > T c . In absence of cross-correlations, the contributions of the two reservoirs are additive in the corresponding Caldeira-Leggett Hamiltonian [62] so that we have in generalization of (5) Accordingly, the cumulant truncation goes through as for a single oscillator and the dynamics of the phase space operators collected in the operator-valued vector Here ξ = (0, ξ l , 0, . . . , 0, ξ r ) t and the matrix M contains the dissipative parts and the effective frequenciesω n according to (25) to effectively account for the nonlinearities and to be determined self-consistently. The detailed structure of M for the Hamiltonian (35) is shown in Appendix B. Having M at hand, the multipartite system obeys the steady-state equations (27), (33), and (34), where the only non-zero elements in L(t−t ) are the second and last diagonal elements containing the auto-correlation functions L l (t) and L r (t), respectively.

d) Energy fluxes
Temperature gradients along the chain induced by the reservoirs at both ends give rise to an energy flux and thus heat transfer. Locally, for the change of energy at an individual site n one must distinguish between the oscillators at the boundaries (n = 1, N ) and the oscillators in the bulk. For the latter, the change in energy follows from and includes also the nearest neighbor coupling, while for the oscillators at the boundaries only one coupling term appears and the other one is replaced by the coupling to the respective reservoir. Calculating the time evolution for each of these parts according to d dt H n tr = L † H n tr leads, after performing the stochastic average, to the dynamics of the on-site energies. This then directly leads to the energy fluxes between adjacent sites and between the oscillators at the ends and the respective reservoirs, namely, d dt where respective heat fluxes follow from One should note that the above averages are trace-and noise-moments which means that the respective currents in the bulk in (41) and (42) are determined by elements of the covariance matrix. The first terms in (43) and (44) represent system-bath correlations and have to be computed as elements of y (30) to obtain the covariance matrix. Therefore, once the covariance matrix is known, the energy currents are implicitly known as well.

e) Rectification of heat transfer
An important measure for the asymmetry in heat transfer is the rectification coefficient. It quantifies the difference of the absolute values of the heat currents for opposite directions of the thermal gradient. The configuration with the left reservoir being hot at temperature T l = T h and the right one being cold at temperature T r = T c < T h is denoted by h → c with the corresponding heat current in steady state j hc . The opposite situation, where only the temperatures are exchanged, i.e. T l = T c < T r = T h , is termed c ← h with the respective heat current j ch . Based on this notation the rectification coefficient reads Since heat always flows from hot to cold, α > 0 (< 0) means that the heat flow T l → T r (T r → T l ) exceeds that from T r → T l (T l → T r ). The basic ingredients for finite rectification are both nonlinearity and spatial symmetry breaking [9,10]. A purely harmonic system does not show any rectification even in presence of spatial asymmetry which implies that non-equidistant energy level spacings induced by the nonlinearities are a crucial pre-requisite. In the present case of anharmonic oscillators, nonlinearities are effectively accounted for by effective state dependent frequenciesω n that depend on temperature and damping.

Test of the approach
a) The classical limit As a first and simple illustration of the presented method, we consider an anharmonic oscillator coupled to a single classical thermal reservoir. In thermal equilibrium, (24) lead to algebraic equations which can be rearranged as Now, only the system-bath correlations have to be calculated to arrive at a quadratic equation for q 2 .

b) Quantum oscillator in thermal equilibrium
We now proceed to analyze the performance of the approach in the quantum regime by considering the properties in thermal equilibrium of a single anharmonic oscillator embedded in a single thermal reservoir. In case of κ > 0 corresponding to a globally stable potential surface, this particularly allows to compare with numerically exact results from a SLED simulation with the full anharmonicity taken into account. Figure 3 displays the variances q 2 and p 2 obtained from both the deterministic method presented here and from the SLED calculation which serves as a benchmark for the cumulant truncation. Apparently, the latter performs very accurately in the considered window of values for κ for both phase space operators even for stronger dissipation. It is also seen (cf. the insets) that the perturbative treatment provides accurate results for higher order momements such as κ q 4 and κ qp 3 +p 3 q 2 , see (26). Note that in the considered range for κ one also clearly sees the impact of friction on the variances: finite dissipation tends to suppress fluctuations in position and to enhance those in momentum. The cumulant formulation nicely captures this feature.  In terms of the effective frequencyω, we find that the approach provides accurate results forω 2 − ω 2 = 3κ q 2 < δ × ω 2 with δ ≈ 1.5 and thus also for stronger anharmonicities. Note that in this domain the harmonic relations apply with effective frequencies though, i.e., q 2 = (h/2mω)coth(ωhβ/2) and p 2 = (mωh/2)coth(ωhβ/2).
In case of a softer mode (κ < 0) a direct comparison with full numerical findings is no longer possible as the anharmonic potential is only locally stable. The perturbative treatment is thus physically sensible only as long as all relevant processes remain sufficiently localized around the minimum of the potential. In particular, the approach does not capture quantum tunneling through the potential barriers from the well into the continuum. However, it provides the nonlinearity required for rectification of heat transfer with the finite dissipation promoting localized states in the well. Figure 4 shows the phase space variances q 2 and p 2 together with the uncertainty product q 2 p 2 . For very weak damping the latter remains at the minimal values 1/2, while position fluctuations increase and momentum fluctuations decrease with growing |κ|. This is in contrast to stronger friction, where already for κ = 0 (harmonic limit) the uncertainty product exceeds 1/2 with strongly suppressed (enhanced) momentum (position) fluctuations for larger |κ|. Note that the range, where the approach is expected to be applicable, is restricted to substantially smaller values of |κ| compared to the situation with κ > 0. We close this analysis by presenting Wigner functions W (q, p) for various values of κ at k B T = 0, see figure. 5. The squeezing of phase space distributions is clearly seen with momentum squeezing for κ < 0 and squeezing in position for κ > 0.

c) Towards quantum heat transfer: anharmonic chains
One advantage of the developed methodology is its high computational efficiency also for large systems. Here, we analyze heat transfer through a chain of anharmonic oscillators of length N = 10 according to (35) terminated at both ends by thermal reservoirs at temperatures T l = T h and T r = T c < T h , respectively. We put all intra- Density plots in figure 6 show the covariance matrices Σ for two anharmonicity parameters κ = −0.10 (top) and κ = 0.50 (bottom). One clearly sees the impact of a temperature gradient as the diagonal elements corresponding to the mode coupled to the hot bath are substantially larger than those attached to the cold one. The temperature dependence of the effective frequencies leads to a more distinct impact of the nonlinearity on the oscillators coupled to the hot bath. Apparently, for these sites the negative κ leads to a broading of the state in position as compared to the case for positive κ.
It is instructive to display the effective frequenciesω n along the chain (cf. figure 7 left) for a long chain N = 50. Apparently, the distribution is rather flat in the bulk and shows deviations only near the interfaces to the reservoirs with the frequencies being smaller (larger) for κ < 0 (κ > 0). We also show in figure 7 (right) the distribution of effective temperaturesT n reconstructed from p 2 n = (mω nh /2) coth[hω n /(2k BTn )]. For further discussions of thermometry for dissipative systems we refer to very resent literature such as in Ref. [64]. Figure 7 reveals a similar profile as that for the effective frequencies, thus indicating ballistic transport. This is confirmed by results for the heat currents (not shown), which do not reveal any dependence on the chain length. The fact that the temperature gradient drops only within narrow interfaces between phonon chain and reservoir is reminiscent of the voltage drop in molecular chains in contact to electronic leads [65].
The question to what extent anharmonicities alone may lead to normal heat flow has led to conflicting results recently depending on the model under consideration [1,11]. For the weak to moderate anharmonicities considered here, a combination of nonlinearity and disorder such as studied for example in [8,7] may provide a mechanism to induce normal heat flow. This is, however, beyond the scope of the current study, where we focus on heat rectification.

Rectification of quantum heat transport: single oscillator
One of the simplest models to realize rectification consists of an anharmonic oscillator coupled to two reservoirs at different temperatures T l and T r . The only way to induce a spacial asymmetry in such a model is to choose different coupling strengths γ l and γ r which remain constant under an exchange of reservoir temperatures. In steady-state, we are free to choose any of the two currents between system and reservoirs as they have identical absolute values. The energy current between the left reservoir and the oscillator reads according to (43) with y (l) p being the momentum part of the system-bath correlation function corresponding to the left reservoir.
In the classical limit, this is y (l) p = mγ l /β l , while p 2 for a system coupled to two baths is determined by additive fluctuations and damping strengths, respectively, i.e., The resulting heat current reads with ∆T = T r − T l . This immediately shows that no rectification can be observed for a single classical mode since an exchange of temperatures would result in a current with identical absolute value j hc = − j ch . The previous argument, however, does not apply to the quantum case, where the system-reservoir correlations depend on details of the Green's function caused by the non-Markovianity of the reservoirs. This Green's function is calculated with the effective frequencyω which varies with an exchange of the reservoir temperatures. Therefore, rectification is possible with a single nonlinear system degree of freedom as shown in refs. [19,21]. Figure 8 shows the rectification coefficient α versus ∆γ = γ r −γ l where γ l = 0.1 is kept constant. Results for positive and negative values of κ over the whole range of ∆γ reveal that κ < 0 leads to α < 0 and κ > 0 to α > 0. This difference in the signs can be attributed to the different level structure of the system, where for the softer mode the energy level spacings are narrower for higher lying states while the opposite is true for the stiffer mode, see figure 1.

Rectification of quantum heat transport: chains of oscillators
Going from a single oscillator to chains of oscillators is not just a quantitative modification of the setting. It is particularly interesting as it allows, by tuning the asymmetry either in the coupling to the reservoirs or in the on-site frequencies, to control the rectification of being positive or negative. This may be of relevance for experimental realizations of heat valves as they have very recently been explored in [33]. Underlying physical principles can be revealed already for a set-up consisting of two oscillators, where we keep frequency and coupling strength of the oscillator at site n = 1 fixed and vary those at n = 2, i.e. ω 2 = ω 1 + ∆ω and γ r = γ l + ∆γ.
Before we discuss specific results below, let us already here elucidate the main mechanism that governs the rectification process. For the Hilbert space of the oscillator chain alone (no reservoirs) two sets of basis functions are distinguished, namely, the one consisting of the localized eigenstates of the individual oscillators and the one consisting of the normal modes of the coupled system. Then, roughly speaking, if the inter-oscillator coupling µ dominates against the oscillator-reservoir couplings γ l , γ r , heat transfer occurs along the delocalized normal modes, while in the opposite situation localized states rule the game. It turns out that the changeover from one regime to the other one by tuning either ∆ω or ∆γ is associated with a sign change in the rectification coefficient α. This also implies that starting from the symmetric situation a growing asymmetry first leads to an extremum of α before one reaches the point α = 0. This discussion can be made more quantitatively by considering the effective spectral density J eff,r (ω) of the right reservoir as seen from the first oscillator (fixed parameters) through the second one (varying parameters). For this purpose, we employ the formalism developed in [66] and obtain for a purely ohmic bath [limit in (4) for ω c → ∞] the expression with the dimensionless inter-oscillator couplingμ = µ/mω 2 1 and effective frequencies ω 2 1/2 = ω 2 1/2 + 3κ q 2 1/2 . It has the expected Lorentzian form with a maximum around ω =ω 2 and reduces in the low frequency regime to an ohmic type of density J eff,r (ω → 0) ≈μ 2 mωγ r . Note that the effective spectral density (51) contains all three relevant parameters: the inter-oscillator coupling µ as well as γ r and ω 2 , which we use to induce spatial asymmetry to the system. Now, let us consider the case with ∆ω = 0 and varying ∆γ such that in the symmetric situation, ∆γ = 0, one has µ > γ l , γ r . As a consequence, for growing ∆γ first the delocalized mode picture applies and the rectification coefficient approaches an extremum around that asymmetry point, where heat transfer at resonance is optimal due to a matching of reservoir coupling constants to oscillator 1, i.e. γ eff,r ≈ γ l with the effective coupling from oscillator 1 to the right reservoir γ eff,r ≡ J eff,r (ω =ω 2 )/mω 2 ≡ µ 2ω2 2 /4γ r , i.e.,μ 2ω2 2 ≈ 4γ l γ r .
With further increasing ∆γ the rectification tends to zero, changes sign, the asymmetry dominates against the inter-oscillator coupling, and a picture based on localized modes captures the heat transfer. We will see below that the sign of the rectification coefficient α also depends on the sign of the anharmonicity parameter κ.

a) Rectification by variation of the damping
Detailed results are now first discussed for the situation, where asymmetry is induced by a varying damping, see Figure 9. As anticipated, one indeed observes a change in the sign of the rectification coefficient together with an extremum for moderate asymmetries. When comparing the location of the extrema with the prediction according to the matching condition (52), one finds an excellent agreement. For the chosen parameters µ > γ l , γ r in the symmetric case, the previous discussion directly applies. The sign change in α can also be understood from the impact of friction onto q 2 1/2 and thus onto the effective frequenciesω 1/2 : For strong coupling to the right reservoir (γ l γ r ) one always has (for identical harmonic frequencies)ω 2 <ω 1 for κ > 0 and ω 2 >ω 1 for κ < 0 due to the strong squeezing in position of oscillator 2. Quantum mechanically, the different energy level spacings then lead to | j ch | > | j hc | and thus α < 0 for κ > 0 and | j ch | < | j hc | with α > 0 for κ < 0. For weak asymmetry (γ l ≤ γ r ) the normal modes couple slightly more efficient to the right reservoir and the bottleneck is the coupling to the left reservoir. Consequently, the relation between the respective heat currents and the sign of α interchanges compared to the strongly asymmetric situation.
It is now also instructive to look for the rectification when the inter-oscillator coupling is tuned, see figure 10: Following our above reasoning the larger µ, the stronger needs the asymmetry to be in order to induce a sign change in the rectification, while for very weak coupling no sign change occurs at all. The mechanism developed above, namely, delocalized modes versus localized modes depending on the strength of the asymmetry, suggests that for purely classical reservoirs a qualitative similar behavior should be present. This is indeed the case as figure 10 (left) reveals.
The impact of the number of oscillators N in the chain on the rectification is displayed in figure 11. With growing N the asymmetry must increase as well in order to induce a sign change in α. We ascribe this to the feature already addressed above (cf. figure 7), namely, an effective screening of the asymmetry such that the bulk remains robust against variations of the coupling to the right reservoir. However, the absolute value of α increases with increasing chain length.

b) Rectification by modulation of the frequencies
As an alternative to a variation of the chain-reservoir coupling, spatial asymmetry can also be induced by varying on-site frequencies ω n . Figure 12 shows α for a system of two coupled anharmonic oscillators with the asymmetry being quantified by ∆ω = ω 2 − ω 1 with ω 1 = 1.0 put constant, while ω 2 is tuned. Instead of analyzing different values for the anharmonicity parameter as above, we here consider a constant temperature gradient at different individual temperatures from the quantum up to the classical regime. One again observes a distinct maximum (κ > 0) which becomes more pronounced towards the high temperature domain. What is more striking though is the behaviour for larger asymmetries: While the classical rectification approaches zero from above, in the low temperature quantum regime we see again a change in sign of α and then a saturation with further increasing ∆ω. This is a genuine quantum phenomenon as an inspection of the covariance matrix reveals, see figure 13. Classically, the phase space entries for the oscillator coupled to the reservoir at T c = 0 are absent, while this is not the case quantum mechanically. In the latter case, ground state fluctuations survive, where those of the oscillator with larger frequencyω 2 exceed those of the softer one. This then gives rise to the observed finite rectification also for large ∆ω. Why does the rectification saturate in this limit? To understand this we again exploit the effective spectral density specified in (51): For ω 2 ω 1 due toω 2 ω 1 the first oscillator effectively probes   only the ohmic-type low frequency portion of the distribution J eff,r (ω ω 2 ) ≈μ 2 mωγ r which is independent of ω 2 . Heat transfer is thus governed by low frequency quantum fluctuations, a process that may also be interpreted as quantum tunneling through the second oscillator.
The detuning of the oscillator's frequency as a resource for rectification suggests an extension to a chain consisting of an even number of oscillators where the frequencies with even index are set to ω high and are varied, while the ones with odd index are kept constant at ω low = 1.0. Such a setting can be considered as a chain of heat valves, where one element consists of two coupled oscillators with ω low and ω high , respectively. Figure 14 shows the rectification for such a series connection of heat valves. Apparently, α < 0 is only obtained for one single valve with N = 2, while a series of multiple valves leads to α > 0, also for large ∆ω and independently of the number of oscillators N . To understand the reason for this quite different behavoir, we consider the situation with two valve elements (N = 4): Then, we have for h → c a significant increase ofω 1 (strong fluctuations in position) which supports heat transport through the whole chain since then alsoω 2 is affected by the hot bath. Instead, for c ← h even a small ∆ω is sufficient to screen the effect of the hot bath onω 3 so that this respective oscillator remains in the ground state. Therefore, the detuning between the oscillators n = 4 and n = 3 is larger than for h → c which implies α > 0 for chains with N ≥ 4. In other words, the length dependence for N ≥ 4 is determined by the presence of detuned modes in the bulk which are not attached to a reservoir, while the total number of those modes is not important.

Rectification in disordered systems
The rather delicate interplay of nonlinearity and spatial asymmetry which the rectification is based on rises the question of the stability of these effects in presence of disorder. In the sequel, we restrict ourselves to a randomization of the on-site frequencies ω n and obtain the rectification as an average over several thousand of realizations, where the heat currents j hc and j ch are calculated with identical realizations. In the left panel in figure 15 we present α dis. for varying ∆γ = γ r −γ l and use normally distributed on-site frequencies with standard deviation σ ω and mean ω n dis. = 1.0, where we reject negative values of ω n . We emphasize that ω n dis. is an average with respect to random chains and not over the individual sites n of one chain. The right panel shows α dis. for varying ∆ω = ω N dis. − ω low dis. with normally distributed on-site frequencies and equal dampings γ l = γ r = 0.1 on both ends.
Obviously, for varying ∆γ (left), the disorder frustrates the zero-crossing of the rectification observed at ∆γ ∼ 2.0 for the ordered case shown in figure 11. The profile of α dis. for stronger disorder even suggests a convergence of α dis. → 0 for large ∆γ. For a system of two oscillators, we found that α < 0 arises from a slightly stronger detuning of the two frequencies for h → c than for c ← h. For the disordered case, this small effect seems to be washed out. Concerning the case of a varying ∆ω shown in the right panel, it is interesting to observe that here disorder is less effective: one has α < 0 for large ∆ω also for stronger disorder. Since in this regime ∆ω is very large, it is clear that the system is more robust against small variations of the ω n . This manifests as well in the relatively small error bars for the random α. Random couplings µ n were also investigated and found to cause a similar behavior as random frequencies (not shown).

Summary and Outlook
In this paper we developed a general framework to describe heat transfer in open quantum systems which is not plagued by inconsistencies of weak coupling approaches, allows to cover a broad range of parameter space (weak -strong coupling, highlow temperatures) and various realizations (static, external driving, disorder), and can be generalized to higher dimensions. It thus may serve as platform to quantify the performance of heat valves, heat engines, or cooling devices as they are currently under study in atomic and mesoscopic physics. More specifically, the heat transfer and rectfication across chains of anharmonic oscillators has been explored. Tuning the level of symmetry breaking by either changing the asymmetry in the chain-reservoir coupling or in the frequency distribution of the oscillators we find the rectification coefficient to pass from positive to negative or vice versa by running through an extremum. A deeper analysis reveals a mechanism, where heat is predominantly carried by non-local modes (weak asymmetry) or localized modes (strong asymmetry) with a smooth turnover between the two scenarios. Similar findings have been reported recently in an experimental realization of a heat valve based on superconducting circuits [33]. While for an asymmetry induced by the chain-reservoir coupling this mechanism also applies to the classical regime, a genuine quantum effect is found for an asymmetry induced by frequency mismatch between adjacent oscillators.
In this situation, a strong asymmetry gives rise to a finite rectification (and thus a finite heat transfer) in contrast to the classical prediction. This finding may have relevance for recent proposals for the fast initializations of quantum bits (cooling) by frequency tuning [36]. + 3 q tr,c q 2 p + pq 2 2 tr,c + 3 q 2 tr,c q tr,c p tr,c + 3 q 2 tr,c qp + pq 2 tr,c + q 3 tr,c p tr,c q 4 tr = q 4 tr,c + 4 q 3 tr,c q tr,c + 3 q 2 2 tr,c + 6 q 2 tr,c q 2 tr,c + q 4 tr,c .
Setting all trace-cumulants of order larger than two equal to zero and respecting the cancellation of some terms in the squared brackets in (A. It can be shown that this system of equations has a stable solution at zero for all trace-cumulants, although q 2 tr,c is a fluctuating quantity. To do so, we linearize (A.3) by setting all products of second trace-cumulants equal to zero and summarize the equations in a matrix-valued equation d dt σ 2 = M 2 σ 2 with σ 2 = ( q 2 tr,c , p 2 tr,c , qp+pq For positive κ, all eigenvalues have negative real part (stable fixed point). For negative κ and weak anharmonicity, this condition is violated only with exponentially small probability. This may result in an extremely small drift, which will be neglected in a similar manner as tunneling out of the metastable well defined by negative κ. With this analysis, we managed to decouple the dynamics of the first-and second trace-cumulants and can describe the steady-state covariance solely in terms of Σ (msc) .
for n = 1, N : withω 2 n = ω 2 + 3κ q 2 n , being the effective frequency where we used q n tr,c q n tr,c ξ = q 2 n which holds in the steady-state.