Nonequilibrium quantum heat transport between structured environments

We apply the hierarchical equations of motion technique to analyzing nonequilibrium heat transport in a spin-boson type model, whereby heat transfer through a central spin is mediated by an intermediate pair of coupled harmonic oscillators. The coupling between each pair of oscillators is shown to introduce a localized gap into the effective spectral densities characterizing the system-oscillator-reservoir interactions. Compared to the case of a single mediating oscillator, we find the heat current to be drastically modified at weak system-bath coupling. In particular, a second-order treatment fails to capture the correct steady-state behavior in this regime, which stems from the $\lambda^4$-scaling of the energy transfer rate to lowest order in the coupling strength $\lambda$. This leads naturally to a strong suppression in the steady-state current in the asymptotically weak coupling limit. On the other hand, the current noise follows the same scaling as in the single oscillator case in accordance with the fluctuation-dissipation theorem. Additionally, we find the heat current to be consistent with Fourier's law even at large temperature bias. Our analysis highlights a novel mechanism for controlling heat transport in nanoscale systems based on tailoring the spectral properties of thermal environments.


Introduction
Understanding and controlling energy transport mechanisms at the microscopic level is of critical importance to the development of nanoscale quantum technologies.Nanoscale devices that exploit the steady-state transfer of energy or particles (electrons, photons, phonons) through a quantum system-including thermal rectifiers [1,2,3], heat engines [4,5,6,7,8,9] and refrigerators [10,11,12,13]-have become an increasingly active topic of research due to their potential technological applications [14].From a more fundamental perspective, the study of such systems also plays an important role in advancing our understanding of quantum thermodynamics [15,16].Despite most activity in these fields remaining largely theoretical, there has been significant progress in recent years towards the realization of quantum thermal machines and, in particular, circuit quantum thermodynamics [17,18], with the thermal diode effect having recently been demonstrated in superconducting quantum circuits [19].The growing relevance of these technologies thus calls for deeper insights into ways of controlling energy transport at the nanoscale.
For boundary-driven systems, i.e. those in which a heat current is driven by local interactions with thermal baths [20], the nonequilibrium spin-boson model (NESB) provides one of the simplest yet most widely applied models of energy transport exhibiting a rich phenomenology.The model consists of a two-level system (spin-1 2 ) coupled to two Ohmic baths at different temperatures, and has found widespread applications to anharmonic molecular junctions [21,22], cold atoms, and biological systems [23].To date, numerous analytical and numerical techniques have been employed to study the NESB, each providing varying levels of insight into its heat transport characteristics.These include perturbative master equations [1,21], nonequilibrium Green's function methods (NEGF) [24,25], path integral techniques [26,27,28,29], the reaction-coordinate method [30,31], the hierarchical equations of motion (HEOM) [32,33,34,35,36,37], and the time-evolving matrix product operator (TEMPO) technique [38], among others.Perturbative master equations constitute a standard tool of open quantum systems theory [39] and have been shown to capture resonant heat transfer in the NESB up to second-order in the systembath coupling.Alternatively, the (nonequilibrium) non-interacting blip approximation (NIBA) [40,41,42,43], nonequilibrium polaron transformed Redfield equation (PT-RE) [44,45,46], and HEOM [34,37] have been used to treat the full counting statistics [47] of the NESB beyond second-order, demonstrating turnover behavior in the steady-state currents under strong system-bath interactions.
Recently, an extension to the NESB was proposed and studied in [31,37,48,49], where a unitary rotation is applied to one of the interaction terms so as to make the system-bath coupling operators noncommutative.This feature introduces a source of asymmetry into the model that facilitates two competing transport mechanisms: (i) one in which heat transfer occurs between baths through direct excitation of the system, termed system currents, and (ii) another pathway independent of the system in which a heat flow occurs directly between the baths.Notably, since the inter-bath current scales quartically with respect to the system-bath coupling parameter λ as λ → 0, it was established that a second-order approach is insufficient to capture the correct steadystate behavior in cases where the system current is zero [31,37].At the same time, other extensions to the NESB were examined in [50,51,52] to study heat transport through superconducting circuits, consisting of a transmon qubit and pair of resonators in a resonator-qubit-resonator configuration (quantum Rabi model) embedded between two thermal reservoirs [53,18].A similar setup whereby the reservoirs are 'filtered' through an intermediate oscillator has also been considered in the context of autonomous thermal machines [54].
Motivated by these studies, in this paper we investigate the transport properties of a different extended NESB-the gapped nonequilibrium spin-boson model (g-NESB)where the usual Ohmic baths are augmented by a pair of coupled harmonic oscillators (see figure 1).To our knowledge, this model has not previously been proposed or studied in the context of heat transport.We analyze the model in both the transient and steadystate regimes using the numerically exact HEOM [55,56,57,58,59,60].This enables us to capture higher-than-second-order effects at weak system-bath coupling, based on exact expressions for the average heat current and noise obtained via a generating functional approach [33,38].Furthermore, we directly compare its transport properties to the quantum Rabi model of [50,51], which is recovered as a special case of the g-NESB when the inter-mode coupling between oscillators is set to zero.
Our key findings are summerized by the following: (i) in the nonequilibrium steadystate the heat current exhibits a λ 4 -scaling to lowest order in λ, contrasting with the λ 2 -scaling of the quantum Rabi model.This reveals that the addition of a second oscillator considerably suppresses the stationary current in the weak coupling limit.(ii) The steady-state noise in both models exhibits the same λ 2 -scaling to lowest order in λ and-unlike the currents-does not exhibit turnover behavior at strong system-bath coupling.(iii) The steady-state current of the g-NESB obeys Fourier's law at large temperature bias outside the linear response regime.
Additionally, to elucidate the quartic scaling behavior we construct perturbative expressions for the stationary current and noise based on the HEOM.We find that since the spectral densities of the effective baths exhibit a localized gap at the oscillator frequencies, the leading-order contribution to the system current vanishes when the system tunnelling frequency is resonant with the gap.Hence, in this case a secondorder treatment fails to describe the correct physics even in the asymptotically weak coupling limit.However, it should be emphasized that the quartic scaling here stems from a different transport mechanism to that explored in [31,37], which as mentioned above, is due to the noncommutativity of the system-bath coupling operators.
The remainder of the paper is organized as follows.In section 2, we introduce the models under study and derive an analytic expression for the spectral densities of the effective baths.Section 3 outlines technical details relating to the implementation of the HEOM, as well as exact and perturbative expressions for the relevant quantities derived within this framework.In section 4, we present numerical simulations of both the heat current and its noise in the transient and steady state regimes, comparing our results to those obtained from the corresponding quantum Rabi model [50,51].Finally, a summary and conclusions is provided in section 5.

Model
We consider the model of an unbiased two-level system coupled to two structured bosonic baths B j , which are initially in equilibrium at temperatures T j , depicted schematically in figure 1.Each bath consists of a pair of coupled harmonic oscillators M 1 j and M 2 j that linearly interact with both the system and an additional thermal reservoir.The Hamiltonian of the system reads with ∆ a tunnelling frequency and σ x,y,z the Pauli spin-1 2 operators.Throughout we choose natural units such that h = k B = 1.The Hamiltonians describing each pair of harmonic oscillators and system-oscillator coupling are given by (i = 1, 2) Here, a i j (a † i j ) is the annihilation (creation) operator of the ith oscillator in the j-th bath with frequency ω 0 , and X i j = a i j + a † i j .The parameters λ j and g denote the coupling strength of S to the second oscillator of the j-th bath and the coupling between oscillator pairs, respectively.We emphasize that the system is only directly coupled to one of the oscillators in each pair and we assume the inter-mode coupling satisfies g ≪ ω 0 .
In addition to the coupling between oscillator pairs, each oscillator interacts with its own thermal reservoir R i j .As such, the Hamiltonian of each bath B j can be written as where and bi,k j ( b † i,k j ) are the corresponding annihilation (creation) operators of the reservoir R i j .The spectral density of each reservoir is defined as Ji j (ω) = π k g2 i,k j δ(ω − ωk j ).By convention, we have also added a renormalization term H C j to each bath Hamiltonian (4), which accounts for the shift in the minimum of each oscillator potential caused by the oscillator-reservoir coupling [43,61].These terms are explicitly given by so that the full Hamiltonian of the model reads For later discussion we refer to the g = 0 case as the ungapped nonequilibrium spin-boson model (ug-NESB).

Equivalent nonequilibrium spin-boson model
Our goal is to understand how the coupling between each pair of intermediary oscillators influences the properties of the heat flux through the system.To achieve this, we focus on a path integral treatment in which dissipative effects from each of the baths on the system are prescribed through the Feynman-Vernon influence functional F j [q, q ′ ] [39, 43,62].The influence functional appears in the propagator J (q f , t|q i , 0) describing the time evolution of the reduced system under the effect of the baths, where S S [q] is the classical system action, and q, q ′ label the forward and backward system paths.Specifically, these paths represent the possible spin trajectories connecting the boundary states q(t) = q f and q(0) = q i , with q ∈ [−1, 1] labelling the possible spin values in the eigenbasis of the system Hamiltonian (1).
The procedure to derive a closed form expression for F j [q, q ′ ] relies on the functional integrals (9) being Gaussian.For non-interacting harmonic environments in thermal states, this is indeed the case and the influence functional reduces to the well-known exponential form with Φ j [q, q ′ ] represented as a double time integral over the forward and backward system paths [see (14)].Since for non-interacting environments the influence functionals can be determined exactly, it will prove convenient for us to then employ a normal mode representation of the baths B in order to facilitate HEOM calculations of the heat currents ‡.The Hamiltonians of the system and baths in such a representation are ‡ We emphasize that such a representation can be freely chosen for the baths B provided that (i) the system-bath interaction remains linear in annihilation and creation operators, (ii) the bath Hamiltonian remains quadratic in these operators, and (iii) the initial bath states are thermal states at temperatures T j .
written as Here, b k j (b † k j ) is the bosonic annihilation (creation) operator for the k-th normal mode of the j-th bath with frequency ω k j .g k j denote the couplings between each normal mode oscillator of the baths and the system.Given that the total initial state is factorized as with T j = 1/β j the bath temperature and Z B j = Tr[e −β j H B j ], the corresponding path integral expressions for the influence functionals F j [q, q ′ ] are [62] where are the bath correlation functions, and are the spectral densities characterizing the interaction between the effective baths B and system S.
Our task now is to derive an analytic expression for the spectral densities J j (ω) through a suitable choice of Ji j (ω).To this end, we rely on a similar ansatz employed by Garg et al [63], in which a single harmonic oscillator interacting with an Ohmic reservoir was shown to map to an effective harmonic bath with Lorentzian spectral density [64,65,66].We therefore choose each of the spectral densities Ji j (ω) characterizing the M − R interactions to be Ohmic, where γ i j is a dimensionless coupling constant, and ω c a cutoff frequency.Notice that our choice of cutoff function is different from the usual exponential-type used in references [63,64]: this is motivated by fact that it not only gives an effective spectral density which is practical to implement with the HEOM for ω c → ∞, but also removes the dependence of J j (ω) on the renormalization terms in (4) [67].The exact details underpinning the transformation between representations and the derivation of the corresponding spectral density ( 16) are included in Appendix A.
Assuming the first oscillator of each bath is decoupled from its respective reservoir, γ 1 j = 0, the effective spectral density is given by the following analytic function of ω, This expression fully characterizes the interaction between the system and baths in terms of the oscillator parameters ω 0 , g, λ j and Γ j = 4γ 2 j ω 0 ; in particular, if g is set to zero so that interactions with the reservoirs are mediated through a single oscillator, we recover a Lorentzian spectral density in accordance with [63].Our result therefore extends the ansatz of Garg et al [63] to a more general, non-Lorentzian environment, which admits a representation in terms of pairs of damped oscillators.
A key feature of ( 18) compared to the single oscillator case is the localized gap introduced into the bath spectrum at the oscillator frequency ω 0 [68].This is depicted in figure 1(b), where J j (ω) is plotted alongside the corresponding Lorentzian spectral density obtained for g = 0. We investigate below how this feature plays a direct role in modulating the heat flux when the system tunnelling energy is tuned to the gap.
The bath correlation functions associated with J j (ω) may be readily obtained by way of contour integration.The resulting expressions are given by where z jk are complex exponents derived from the pole locations of the bath spectra in the lower half complex plane, and c jk are determined from their corresponding residues.
To implement the HEOM in practice it is necessary to truncate the expansion (19) to a finite number of terms.To this end, we proceed to decompose the bath correlations as where C 0 j (t) comprises a finite sum over the pole contributions from J j (ω), and M j (t) an infinite sum over the Matsubara frequencies ν j n = 2πn β j (n ∈ Z + ): Since terms with exponents ν j n t ≫ 1 make almost no contribution to C j (t) over the relevant time scales, we may neglect them so that C j (t) ≈ K j k=0 c jk e −iz jk t .Given also that the Matsubara frequencies ν j n are proportional to T j , one can expect the number of terms needed to accurately approximate (19) to grow with decreasing bath temperature (in fact, in the limit β j → ∞ the number explodes since the sum in (21) converges to an integral).In all numerical simulations we fix the bath temperatures to be large enough so that truncating M j (t) to a single term produces converged results.

Hierarchical equations of motion
Based on the finite expansion of the correlation functions (19), the system-bath interactions may be treated nonperturbatively using the HEOM technique [55,56,57,58,59,60].In this approach, the non-Markovian effects of the baths B j on S are encoded in a collection of auxiliary density operators (ADOs).These together with the reduced system density matrix ρ S (t) = Tr B [ρ tot (t)] obey a set of coupled equations of motion.The ADOs in the path integral representation are defined as where Here, each of the ADOs is labelled by a pair of multi-index vectors {⃗ n, ⃗ m} = {(n 11 , ...) T , (m 11 , ...) T }, with (n jk , m jk ) ∈ {0, ..., N max j }, and N max j being a maximum cut-off that truncates the depth of the hierarchy associated with each bath.The reduced system density matrix is obtained as the first member of each hierarchy with the index values |⃗ n| = | ⃗ m| = 0, i.e. ρ S (t) = ρ ⃗ 0, ⃗ 0 , and |⃗ n| = jk n jk .The ADOs corresponding to the L-level of the hierarchy are those with indices satisfying |⃗ n| + | ⃗ m| = L. Furthermore, the length of each of the vectors ⃗ n and ⃗ m is determined as the sum K 1 +K 2 of the number of terms kept in the expansion of each bath correlation function.
If we further introduce a set of scaled ADOs ρ′ in accordance with [69], then the HEOM can be derived in the form (the hat indicates a Liouville superoperator representation) Here, L S ρ⃗ n, ⃗ m = [H S , ρ⃗ n, ⃗ m ] is the system Liouvillian, and are superoperators that relate ADOs with index (⃗ n, ⃗ m) to those higher (+) or lower (−) in the hierarchy with indices n jk and m jk incremented by one, i.e. ⃗ n ± jk = (n 11 , ..., n jk ± 1, ...) T .

Heat currents
We now discuss the calculation of the relevant thermodynamic quantities within the HEOM framework.Introducing the current operators I B j = i[H B j , H I j ] and I S j = −i[H S , H I j ], we identify heat with changes in the average energy of each bath [70,33], where defines the mean heat current for the j-th bath, and is the corresponding system heat current.Note that the time dependence of operators inside the averages ⟨...⟩ is generated in the Heisenberg picture.The minus sign in (27) indicates that a positive heat current is associated with a net flow of energy from the bath into the system.Using a generating functional approach [33,38,71], the heat currents may be evaluated directly in terms of the ADOs extracted from the first level of the HEOM [35,72], as detailed in Appendix B: Note that in the nonequilibrium steady state ρ tot → ρ ss tot we have ∂ t ⟨H I j (t)⟩ = 0, and so ⟨I S j ⟩ ss = ⟨I B j ⟩ ss .From the first law of thermodynamics, ⟨I B 1 ⟩ ss = −⟨I B 2 ⟩ ss , it follows that the symmetrized heat current obtained in the asymptotic limit t → ∞, is positive by definition.According to the second law of thermodynamics, the bath currents must also satisfy such that the rate of entropy production in the steady state is non-negative.Higher order moments of the operators I S j and I B j , which represent fluctuations in the heat currents, may also be calculated via the same generating functional approach.In particular, based on the general expressions for the n-th order moments derived in Appendix B, the noises relating to the currents ⟨I S j ⟩ and ⟨I B j ⟩ are determined as [35] where ⃗ 0 + jk,jk ′ denotes ADOs with indices 0 jk and 0 jk ′ raised by one, and Hence, the summations in ( 33) and ( 34) run over all ADOs extracted from the second level of the HEOM.This trend also extends to higher moments; namely, to evaluate ⟨I n S j ⟩ and ⟨I n B j ⟩, ADOs up to the n-th level of the HEOM are needed in general.

Perturbative treatments
The hierarchical structure of the HEOM enables a perturbative solution to the heat currents to be systematically constructed up to a specified order in the bath coupling [73].To demonstrate this, the ADOs for |⃗ n| + | ⃗ m| ̸ = 0 may first be expanded as a formal series in the dimensionless parameter λ, where the tilde indicates a representation in the interaction picture (within Liouville space).Accordingly, this gives a perturbative expansion for the steady state current [see (B.16) in Appendix B]; Note that the leading-order contribution to the steady state current can be derived by restricting the HEOM to first-tier ADOs only.This is because ADOs scale as O(λ 2(|⃗ n|+| ⃗ m|) ) to lowest order in λ [this can be checked by replacing c jk → λ 2 c jk in ( 22)], and hence for |⃗ n| + | ⃗ m| ≥ 2 only contribute to (37) beyond second order.As such, we may determine I (1)  ss by solving the truncated first-tier equations, d dt ρ(1) Notice we have also replaced ρS (τ ) → ρS (t) inside each integral since this falls at the same level of approximation made in equation (38a), i.e. ρS (τ ) = ρS (t) + O(λ 2 ).By then inserting these expressions into (29), we obtain the leading term in the expansion for I ss , where is the second-order dissipator of the j-th bath.Equation ( 40) agrees with a Redfield treatment of the currents [1] and describes local heat flows between the system and baths in the nonequilibrium steady state when the baths act additively.The Redfield equation for ρ S (t) can be written as [74] ∂ρ Higher order contributions capturing non-additive effects may also be evaluated [75], with the corresponding fourth-order expressions given in Appendix C. Defining the second-order noise quantity S ss ≡ ⟨I 2 S 1 ⟩, we may construct another perturbative expansion akin to (37), where to order λ 2 , This result is consistent with the prediction for the current noise from the fluctuationdissipation theorem [43]; that is, S (1)  ss represents the contribution to the noise from equilibrium fluctuations.
In the following we look to compare the predictions of the HEOM with the perturbative approach to assess the role of higher order effects (i.e.those beyond second order) on the behavior of the heat flux.

Results
In this section we present numerical results for both the heat currents ( 29)-( 30) and second-order moments (33), (43) obtained for the g-NESB and ug-NESB.As stated above, these models correspond to when the system is coupled to the reservoirs either through a pair of interacting oscillators, or a single oscillator, respectively.For simplicity we will assume the couplings to the two baths to be equal in both cases, λ 1,2 = λ, and for the system interactions to be resonant with each of the baths, ω 0 = ∆ §.All simulations are performed using the HEOM solver of the Python package QuTiP [76,77].

Transient currents
In figures 2 and 3, we show plots of the transient heat flux for different values of λ, with the spin initalized in the positive z-direction, ⟨σ z (0)⟩ = 1. Figure 4 shows the time evolution of the system coherences for the same set of bath parameters.At strong coupling √ πλ = 0.5∆, we see that the energy currents for the g-NESB and ug-NESB show qualitatively similar behavior; in particular, the finite contribution of § Far from resonance, the currents in the two models are expected to exhibit the same behavior since both the gapped and ungapped spectral densities fall off as O(1/ω 3 ).the interaction energy at these coupling strengths plays a similar role in the heat current dynamics for the two cases.The rate at which the heat currents reach their stationary values, however, is generally slower for the baths with gapped spectra.This behaviour is also reflected in the system coherences which undergo relaxation on a similar time scale as the heat currents.
Conversely, as the coupling strength is further decreased to √ πλ = 0.01∆, transient characteristics of the heat flux begin to deviate for the different bath configurations; see figure 3.Here we only plot the system heat currents since the interaction energies are negligible for all times.Interestingly, we now observe a decoupling between time scales determining the relaxation of the heat currents and the spin coherences for the g-NESB-i.e. the coherences relax much more slowly towards their steady-state than the corresponding current.This is exemplified in figure 4(b) where the real part of the off-diagonal elements of ρ S (t) is displayed alongside it's asymptotic value (red dashdotted line).Additionally, at weak coupling, the gapped baths induce a much smaller rate of energy transfer in the long time limit, contrasting with the ungapped case in which the HEOM and Redfield results are in good agreement (not shown).
Figure 5 displays the scaled dynamical quantity ⟨I 2 S j ⟩/C j (0) characterizing the strength of out-of-equilibrium fluctuations relative to those in equilibrium.Note that our analysis here is restricted to second moments of I S j only, since bath moments beyond  ⟨I B j ⟩ diverge for the spectral densities under consideration.At moderate (strong) system-bath coupling, the scaled noise for the two models exhibits similar short time behavior, but reaches its steady-state value on different time scales.The difference in relaxation times becomes more substantial in the weak coupling regime, although in contrast with the heat currents in figure 3, the noise relaxes faster for the g-NESB than for the ug-NESB.The transient noise for the g-NESB additionally shows quantum beats which are not present for the ug-NESB.

Steady state current
To gain further insight into the differences in the behavior of the heat currents for the two setups, in figure 6(a (37).The other bath parameters are the same as in figure 2.
the coupling strength.The solid lines represent the 'exact' results obtained from the converged HEOM, while open points and dashed lines represent the perturbative results of the truncated HEOM approaches.The acronyms HEOM-1 and HEOM-2 reference the HEOM result evaluated up to second and fourth order in the coupling-i.e.I (1)  ss and I (2)  ss -respectively.
In the weak and ultra-weak coupling regimes, the scaling of the stationary current with λ is noticeably different for the g-NESB and ug-NESB.In the latter case, the current scales quadratically in the coupling strength, I ss ∼ λ 2 , and thus agrees with the HEOM-1 results which treat the heat current up to second-order in the coupling.On the other hand, for baths with gapped spectra, the current scales as I ss ∼ λ 4 .The HEOM-1 therefore cannot capture the steady-state behavior of the current for this case, which must be treated to at least fourth-order in λ.In particular, a second-order approach predicts no heat transfer here at all, since the dissipator (41) vanishes at the system Bohr frequencies, i.e J j (±∆) = 0 and J j (0) = 0.
At strong system-bath coupling, the energy currents exhibit a turnover behavior around √ πλ = 0.1∆, at which point the perturbative approaches no longer agree with exact results.Deeper into the turnover region the two solid lines also converge, indicating that spin transitions mediating the heat flow become largely insensitive to the gap in the bath spectra when λ becomes large.Next, we display the steady state noise S ss in figure 6(b) plotted against the coupling strength λ for the two bath configurations.It can be seen that the g-NESB and ug-NESB display similar noise characteristics in contrast with the different scaling profiles observed for the heat currents.In particular, at weak system-bath coupling the noise The solid blue and green lines correspond to the exact HEOM results for g = 0.1∆ and g = 0, while the black dashed line represents the Redfield results.The inset shows the result of the simulated g-NESB against Fourier's law I F ss = κ∆T .All other parameters are the same as in figure 6.
scales quadratically with λ in both cases as predicted by the fluctuation-dissipation theorem.Furthermore, the noise unlike the current I ss does not exhibit turnover behavior in either of the two models at strong coupling.As such, the perturbative and exact expressions for S ss also adopt a similar functional form in the turnover region of the current: discrepancies only occur between approaches for coupling strengths above λ ≈ 0.1∆.Lastly, figure 7 illustrates how the stationary currents for each model depends on the temperature bias ∆T = T 1 − T 2 , assuming the temperature of the cold bath to be fixed to T 2 = ∆.We compare the exact HEOM simulations with those obtained from the second-order Redfield equation (41) for moderate (strong) and weak system-bath coupling.When the baths are close to equilibrium ∆T ≪ 1, we find the current scales linearly with ∆T for moderate and weak coupling as would be expected from Fourier's law, I F ss = κ∆T , where is the thermal conductance.As the temperature difference grows in magnitude, however, the two models display differing behavior.Within the ug-NESB, the heat current saturates as ∆T is tuned away from the linear response regime.This saturation also occurs in the g-NESB but the current approaches it's asymtoptic value at a slower rate, thereby causing Fourier's law to hold over a wider temperature range than for the ug-NESB, particularly at weaker coupling.We note further that while the Redfield equation ( 40) predicts a similar turnover as the HEOM, it tends to substantially overestimate the heat flow in the ug-NESB.

Summary
We have analyzed heat transport in a structured-environment version of the NESB, where energy exchange between the system and Ohmic reservoirs is mediated by pairs of coupled harmonic oscillators.Compared to the case of a single mediating oscillator considered in [50,51], the transport properties are found to be significantly altered in both the transient and steady-state regimes (with the exception of the steady-state current noise), specifically at weak system-bath coupling.In particular, the steadystate heat current displays a non-trivial scaling that cannot be captured using standard perturbative techniques including the second-order master equation (42).To rigorously treat higher-than-second-order effects we have used the numerically exact HEOM, which facilitates the calculation of arbitrary moments of system and bath current operators in terms of ADOs.For example, averages depend on ADOs extracted from the firsttier of the HEOM [72].The calculation of higher moments has also been achieved through extending the generating functional technique used by Kato and Tanimura [33] to moments beyond first-order on par with the path integral treatment of [35].It should be further noted that our implementation of the HEOM is distinct from that developed by Cerrillo et al [34], which is based on the full counting statistics formalism [47].Overall, our results highlight the possibility of using reservoir engineering to nontrivially influence the properties of thermal transport through a quantum system.In this aspect, practical applications may be found in the design of nanoscale heat devices.In particular, analogous models investigating heat transport through superconducting circuits, where a similar resonator-qubit-resonator architecture is considered, indicate the potential for our model to be implemented on such a platform.
Future work could potentially examine the effect of asymmetry on the transport properties of the model, including rectification of the currents that arises with noncommutative system-bath coupling operators [31,37] and-or asymmetric coupling strengths [1,3], and in particular, whether the diode effect is enhanced or diminished for environments possessing gapped spectra.Alternatively, generalisations of the model to include couplings between components of the different baths in the spirit of [52,53] could be analyzed.Finally, applications to continuous thermal machines, such as autonomous refrigerators and heat engines, may also be of interest.between the two bath representations is exact.Therefore, under analytic continuation of K(ν) to the lower half complex plane, ν = ω − iε, the spectral density J(ω) is determined as

Appendix B. Quantum heat currents in terms of ADOs
In this appendix we outline the derivation of the heat currents following a generating functional approach.The approach we implement has previously been employed in [33] to derive analogous expressions for averages of system and bath heat currents within the HEOM framework.Here, we extend their formulation to also include the derivation of higher-order moments in accordance with [35,72].

Appendix B.1. Generating functional
To fix the notation, we first note that the total density matrix in the interaction picture evolves according to ρtot (t) = T {e −i j t 0 ds LI j (s) }ρ tot (0), (B.1 where T is the chronological time ordering operator [39], LI j (t)ρ = [ HI j (t), ρ] is the Liouville superoperator associated with the interaction Hamiltonian HI j (t) = Ṽj (t) Bj (t) [see (12)], and V j = V † j an arbitrary system operator.We recall that the tilde indicates a time dependence within the interaction picture, Õ(t) = e i(H S +H B )t Oe −i(H S +H B )t .Moreover, by defining the left and right acting superoperators, The evaluation of the functional derivative above requires us to initially obtain a closed form expression for ρB S (ϕ, t).To this end, we note that since each of the baths obeys Gaussian statistics, the righthand side of (B.5) can be evaluated using Wick's theorem, Here, we've used the idempotency of the time ordering superoperator, T 2 = T , as well as that ⟨ j (•)⟩ B = j ⟨•⟩ B j , and ⟨ LI j (t)⟩ B j = 0.By now using Having obtained an exact expression for ρB S (ϕ, t), one can next compute it's functional derivative according to [80] Here, η j (t) is a suitable test function, and F[ϕ] is the Feynman-Vernon influence functional incorporating the source terms ϕ, Let us denote by R = {(s, u) : 0 ≤ s ≤ t, 0 ≤ u ≤ s} the triangular region of integration in (B.13).Over this region the integrals containing the terms η j (u) may be rewritten as so that by substituting these expressions back into (B.13),we can directly read off the functional derivative of ρB S (ϕ, t).In doing so, we get δ δϕ j (t) ρB S (ϕ, t) = −T which in turn gives us the following expression for the system heat current, We can now insert the exponential decomposition of the bath correlation functions (19) into the above to obtain The derivation of all higher order currents can be achieved through evaluating derivatives of the moment generating functional (B.5) for n > 1.In particular, the n-th order moment of the system heat current is determined from We now proceed to evaluate ⟨I n B j ⟩ systematically starting with the n = 2 case.In this regard, we note that the second derivative of ρS (ϕ, t) may be formally written as [see (B.15)] It is straightforward to verify that By then inserting this expression into (B.17) with n = 2, one obtains with Λ j [ϕ j = 0] = Λ i (t).The quantity Λ 2 j (t) may be expanded using the binomial theorem and substituted into (B.23) to give where the 0 jk,jk ′ notation references second-tier ADOs with indices 0 jk and 0 jk ′ raised by one.
Considering moments beyond second order, let us employ the following shorthand notation The (n + 1)-th derivative of ρB S (ϕ, t) can be shown to satisfy The terms Λ ℓ j [ϕ j ] may then be expanded in a similar way to Λ 2 j (t) to obtain ⟨I n S j ⟩ as a sum over ADOs.By doing so, one finds such that n-th order moments of the system current depend on ADOs up to the n-th level of the HEOM.

Appendix B.3. Bath heat currents
To calculate the bath heat current (30) under the same approach, we insert Õj (t) = Ãj (t) + into (B.6) with A j = i j ω k g k j (b † k j − b k j ).This allows us to write I B j in terms of the generating functional ρA S (ϕ, t), The evaluation of ρA S (ϕ, t) follows the same procedure used above to derive a closed form expression for ρB S (ϕ, t) [35].After some algebra, we end up with The functional derivative of ρA S (ϕ, t) may also be evaluated in a similar way to ρB S (ϕ, t), resulting in δ δϕ j (t) ρA S (ϕ, t) = −T To now write ⟨I B j ⟩ in terms of the ADOs, we can expand D j (t) and D * j (t) into sums of exponentials using (19)   where evidently the same relationship holds between the n-th order bath current ⟨I n B j ⟩ and ADOs as with ⟨I n S j ⟩.

Figure 1 .
Figure 1.(a) Schematic of nonequilibrium heat transport model with structured baths.The system S couples to hot and cold thermal reservoirs R with Ohmic spectral densities (17) through an intermediate pair of coupled oscillators M .Each oscillatorpair and their reservoirs are assumed to be intialized in a composite Gibbs state ρ B (0) = j=1,2 e −βj H B j /Z Bj .(b) Plots of the bath spectral densities corresponding to the ungapped NESB (ug-NESB) and gapped NESB.

1 I S 2 Figure 3 .
Figure 3. Transient system currents for (a) the g-NESB and (b) the ug-NESB at weak system-bath coupling √ πλ = 0.01∆, where g = 0.1∆ and g = 0, respectively.The black dotted lines in the inset of panel (a) indicate the steady-state values of the corresponding currents.All other parameters are the same as in figure 2.

Figure 4 .
Figure 4. System coherence for the g-NESB and ug-NESB at (a) strong system-bath coupling √ πλ = 0.5∆, and (b) weak coupling √ πλ = 0.01∆.Re[ρ 10 ] and Im[ρ 10 ] are the real and imaginary parts of the system coherence in the computational basis.The red dash-dotted line indicates the steady-state value of Re[ρ 10 ] in the g-NESB.All other parameters are the same as in figure 2.

Figure 6 .
Figure 6.Stationary (a) heat current and (b) noise S ss = ⟨I 2 S1 ⟩ ss as a function of the coupling strength λ for the g-NESB and ug-NESB.In (a), the green open circles and black open triangles represent perturbative results for the steady state current obtained from(37).The other bath parameters are the same as in figure2.

Figure 7 .
Figure 7. Stationary heat current as a function of the temperature difference ∆T = T 1 − T 2 for (a) weak coupling √ πλ = 0.01∆, and (b) moderate (strong) coupling √ πλ = 0.1∆, with T 2 = ∆.The solid blue and green lines correspond to the exact HEOM results for g = 0.1∆ and g = 0, while the black dashed line represents the Redfield results.The inset shows the result of the simulated g-NESB against Fourier's law I F ss = κ∆T .All other parameters are the same as in figure 6.