Strong coupling effects in quantum thermal transport with the reaction coordinate method

We present a semi-analytical approach for studying quantum thermal energy transport beyond the weak system-bath coupling regime. Our treatment, which results in a renormalized, effective Hamiltonian model is based on the reaction coordinate method. In our technique, applied to the nonequilibrium spin-boson model, a collective coordinate is extracted from each environment and added into the system to construct an enlarged system. After performing additional Hamiltonian's truncation and transformation, we attain an effective two-level system with renormalized parameters, which is weakly coupled to its environments, thus can be simulated using a perturbative Markovian quantum master equation approach. We compare the heat current characteristics in our method to other techniques, and demonstrate that we properly capture strong system-bath signatures such as the turnover behavior of the heat current as a function of system-bath coupling strength. We further investigate the thermal diode effect and demonstrate that strong couplings moderately improve the rectification ratio relative to the weak coupling limit. The effective Hamiltonian method that we developed here offers fundamental insight into the strong coupling behavior, and is computationally economic. Applications of the method towards studying quantum thermal machines are anticipated.

We present a semi-analytical approach for studying quantum thermal energy transport beyond the weak system-bath coupling regime. Our treatment, which results in a renormalized, effective Hamiltonian model is based on the reaction coordinate method. In our technique, applied to the nonequilibrium spin-boson model, a collective coordinate is extracted from each environment and added into the system to construct an enlarged system. After performing additional Hamiltonian's truncation and transformation, we attain an effective twolevel system with renormalized parameters, which is weakly coupled to its environments, thus can be simulated using a perturbative Markovian quantum master equation approach. We compare the heat current characteristics in our method to other techniques, and demonstrate that we properly capture strong system-bath signatures such as the turnover behavior of the heat current as a function of system-bath coupling strength. We further investigate the thermal diode effect and demonstrate that strong couplings moderately improve the rectification ratio relative to the weak coupling limit. The effective Hamiltonian method that we developed here offers fundamental insight into the strong coupling behavior, and is computationally economic. Applications of the method towards studying quantum thermal machines are anticipated.

I. INTRODUCTION
Nonequilibrium dissipative quantum systems play a central role in quantum technologies. However, simulating the dynamics and transport properties of open quantum systems is challenging given the important role of the surrounding environment. Perturbative treatments in the language of Green's functions 1 or quantum master equations (QMEs) 2-4 have proved fruitful, bringing much insights. QME methods are perturbative in the coupling energy of the system to its surroundings, further typically assuming a Markovian-memoryless system's dynamics. Particularly for second-order QMEs, the rate constants responsible for transitions between quantum states are linear in the system-bath interaction energy (quantified e.g. by the reorganization energy).
Nevertheless, many open-system processes fall outside the regime of validity of perturbation techniques, ranging from biological charge and exciton transport processes 5,6 to quantum devices 7 .
In contrast to the linear-monotonous nature of the weak coupling regime, strong system-bath couplings can enact rich functional behavior relying on involved phenomena such as time nonlocal effects, cooperative interactions between thermal baths, and the buildup of correlations between a system and its surroundings. Significant efforts are currently directed to understandand harness-the nontrivial outcome of strong system-bath couplings on the performance of quantum heat machines; a small sample of such studies include [8][9][10][11][12][13][14] .
The nonequilibrium spin-boson (SB) model, with a two-level system (spin) coupled to two independent harmonic-oscillator (HO) baths was suggested for exploring the fundamentals of thermal energy transport in anharmonic nanojunctions 15,16 . Remarkably, the model was recently realized in superconducting quantum circuits demonstrating a thermal diode effect 17 . The heat current in the SB model displays a turnover behavior at strong coupling: While in the weak coupling limit the steady state heat current grows monotonically with the system-bath interaction energy, at strong coupling the current is suppressed [18][19][20][21][22][23] .
Numerous techniques were developed to investigate the thermal transport characteristics of the nonequilibrium SB model, specifically the turnover behavior mentioned above, by going beyond the weak-coupling Born-Markov Redfield (BMR) equation. For example, the non-interacting blip approximation (NIBA), which is perturbative in the tunneling splitting (nonadiabatic parameter) can properly describe strong system-bath coupling effects for Ohmic baths 15,16,19,21,24 . To interpolate between the weak-coupling Born-Markov Redfield equation and the strong-coupling NIBA for general spectral functions, a nonequilibrium polaron-transformed Redfield equation has been developed in Refs. 25,26 . Keldysh nonequilibrium Green's function (NEGF) methods were furthermore developed for the SB model in Refs. [27][28][29][30] .
The reaction coordinate (RC) method offers an alternative, balanced approach between low-order perturbation theory treatments and brute-force numerical methods. In this semi-analytic method, a collective coordinate of the bath, termed reaction coordinate, is identified and included into the quantum system. This step, or mapping, is exact. Since the RC incorporates system-bath interaction energies, the dynamics of the enlarged, "super-system" captures effects beyond weak coupling and beyond the Markovian approximation even when evolved approximately using a second-order Markovian QME 48,49 . Obviously, the principle of the RC method can be generalized with several discrete modes extracted from the baths and included either as part of the system or as a primary non Markovian environment, further coupled to a secondary Markovian bath [44][45][46][47][48][49][50][51][52] .
Apart from its computational simplicity, the reaction-coordinate quantum master equation (RC-QME) procedure can bring insights on the system's dynamics and the evolution of system-environment correlations 53,54 . The RC-QME method has been applied onto both bosonic and fermionic systems [55][56][57][58] . Furthermore, recent studies with this method examined the role of strong coupling effects in the performance of autonomous 8 and four-stroke 57 thermal machines.
In this paper we develop, test, and benchmark a RC quantum master equation method for nonequilibrium quantum heat transport problems. Beyond numerical simulations, our objective is to gain insight on the role of strong system-bath couplings in the operation of thermal devices.
The turnover behavior of the heat current with system-bath coupling energy has been demonstrated with several methods, interpreted using a picture analogous to polaron transformation 15,16,19,21,25,26 . Here, we explain the turnover behavior of the heat current with an alternative picture: Employing the RC framework, we show that with a series of transformations and truncation we can convert the original-standard spin-boson model with strong system-bath couplings into an effective spin-boson model (EFF-SB). In the new model, the spin splitting is renormalized (suppressed) with coupling energy. Other parameters that are renormalized or transformed are the coupling parameters of the spin to the heat baths and the reservoirs' spectral density functions. The new, effective SB model captures nontrivial-nonlinear effects-such as the suppression of the heat current at strong coupling-even at the level of a Markovian second-order QME. Overall, our approach handles the nontrivial strong coupling regime with minimal cost.
Our work addresses the following questions regarding heat transport in nanojunctions: (i) What is the origin of the heat current suppression at strong coupling? (ii) How does the current scale with different system's parameters? (iii) What is the extent of the thermal rectification (diode) behavior? Does strong coupling assist the rectification effect?
The paper is organized as follows. In Section II, we describe the standard spin-boson model. We further describe the RC procedure, which takes us to an effective spin-boson model with renormalized parameters. We describe the RC-QME method in Sec. III. Simulations are presented in Sec. IV where we focus on gaining insights on the role of strong coupling in transport and on benchmarking the RC-QME technique against other methods. In Sec. V, we present additional simulations of the asymmetric spin-boson junction, and analyze its operation as a thermal diode in the strong coupling regime. We conclude and discuss future directions in Section VI. Figure 1. A two-level system interacting with two bosonic reservoirs held at different temperatures is mapped to an enlarged supersystem comprising the original two-level system coupled to two collective modes (one for each reservoir), which are in turn coupled to residual heat baths.

II. MODEL: DERIVATION OF THE EFFECTIVE SPIN-BOSON HAMILTONIAN
We begin with the nonequilibrium spin-boson model allowing for strong coupling between the quantum system (spin) and the baths, which are characterized by peaked, Brownian oscillator spectral functions. After a series of transformations and truncations, the model is reduced into an effective SB model, which now can be assumed to be weakly coupled to the bathsdescribed by modified spectral functions. The energy spacing of the effective spin and its coupling parameters to the baths are renormalized, thus providing understanding on the role of strong couplings in thermal transport.

A. The standard SB model
We start with a two-level system coupled to two heat baths (i = h, c) comprising collections of harmonic oscillators (h = 1, k B = 1), Here, ε is the energy splitting of the spin, ∆ is the tunneling frequency, σ x,y,z are the Pauli matrices, c † k,i (c k,i ) are the bosonic creation (annihilation) operators for the hot and cold baths (i = h, c) for a mode of frequency ν k,i . The displacements of the harmonic oscillators are coupled to the spin polarization with strength f k,i , and these interactions are captured by the spectral density functions, In what follows, we refer to the model Hamiltonian Eq. (1) as the "standard spin-boson model" (SSB).
Technically, it is convenient to first specify the spectral function of the residual baths, after the RC mapping (as we do in Sec. II B), then work our way back and derive the spectral function of the SSB model. For completeness, we already present here the spectral function of the SSB model: It can be shown that an Ohmic function in the RC representation corresponds to a Brownian-oscillator model for the SSB model, (2) We elaborate on this mapping in Appendix A. The physical picture emerging from this spectral function is that, in each bath, there is a collection of modes with frequencies centered around Ω i that are strongly coupled to the system. The energy parameter λ i quantifies the coupling strength of the spin to the ith bath. The combination γ i Ω i , with γ i a dimensionless parameter, provides a measure for the width of the spectral function. Since bath modes around Ω i dominate the Brownian spectral function, by extracting an oscillator around this frequency from the ith bath, and including it into the system, we can capture strong systembath coupling effects.

B. Reaction Coordinate mapping
We follow a procedure analogous to that described in Ref. 53 , but extended here to the case of a system strongly coupled to two heat baths. In brief, we apply a normal-mode transformation to the harmonic part of Eq. (1) and define a collective coordinate for each environment, termed reaction coordinate. The two RCs (i = h, c) are each coupled directly to the spin, as well as to their residual harmonic bath. This mapping is sketched in Fig. 1. The resulting reaction coordinate (RC) Hamiltonian is given by The two collective coordinates are defined such that, Here, λ i is the newly-defined coupling of each RC to the spin system, Ω i are the corresponding frequencies of the reaction coordinates. The residual baths are denoted by creation (annihilation) operators b † k,i (b k,i ), which now couple only to the ith RC. The spectral density functions of the residual heat baths are defined as J RC,i (ω) = ∑ k g 2 k,i δ (ω − ω k,i ). It is important to note that this mapping is exact. To complete the mapping, the spectral density functions of the SSB and the RC models need to be coordinated such that the spin system follows a corresponding dynamics. For simplicity, we assume that the spectral functions of the residual baths are Ohmic, Here, γ i is a dimensionless coupling parameter of the RC to the residual bath. This parameter is assumed to be small in our simulations allowing us to use the second-order QME for the model after the RC mapping. Λ i is the cutoff frequency of the bath.
The reaction coordinate Hamiltonian, Eq. (3), comprises the following components: The enlarged system (supersystem), now including the original spin, two additional harmonic oscillators as the RCs, and the coupling of the spin to these two oscillators. Other terms in Eq. (3), appearing in the second line are the two residual heat baths (i = h, c) and their couplings to the two reaction coordinates. The last element (third line) is a nontrivial quadratic term, . In Ref. 53 it was treated as part of the interaction Hamiltonian and shown to be partially eliminated from the Markovian quantum master equation. In Ref. 59 , the exactly solvable model of a chain of harmonic oscillators coupled to harmonic baths was treated with the Langevin Equation, then compared to the approximate Lindblad method. It was pointed out that standard QME methods (such as in Lindblad form or Redfield) have to be executed without the quadratic term.
Essentially, an exact quantum equation of motion for the Hamiltonian (3) is reduced to a Markovian QME in which the bare parameters of the RC and its coupling to the system and bath (Ω, λ , g k,i ) are to be used, rather than renormalized parameters defined once the quadratic term is absorbed 60 . In Appendix B we perform a rotation transformation and absorb the quadratic term, by renormalizing the parameters of the reaction coordinate Hamiltonian. However, as pointed out in Refs. 53,59,60 , Markovian QME simulations in this rotated basis do not correspond to the original Hamiltonian.

C. Truncation of the RC Hamiltonian
Simulations of the RC model, Eq. (3) are performed by truncating it to include only M equi-distant levels. We count these levels using an index l, For simplicity, we include the same number of states for each truncated RC. After the truncation, the Hamiltonian describes an enlarged system (ES) with a spin coupled to two M-level manifolds, with each M manifold coupled to a residual heat bath, Explicitly, the enlarged system Hamiltonian is As before, the residual baths are given as a collection of harmonic oscillators, where we identify the system operator V M ES,i and the coupled bath operator B i as The Hamiltonian H M ES and system's operators V M ES,i are matrices of dimension 2M 2 × 2M 2 .

D. Diagonalization and the effective spin-boson model
As a quick review, we started from the SSB model in Sec. II A. We extracted RCs, one from each bath in Sec. II B and we truncated the RCs in Sec. II C. The last two steps of the procedure are described in this Section. They involve: (i) Diagonalizing the M-level RC extended-system H M ES of Eq. (8). (ii) Truncating the resulting spectrum, leaving only two states-thus strikingly, finishing where we started with a two state system coupled to harmonic oscillator heat baths. Albeit, as we show next, the spin-boson model at the end of the procedure relates to the SSB model through renormalized parameters. In more details: We diagonalize H M ES of Eq. (8) with a unitary transformation U. In general, this procedure is done numerically. This transformation is further applied onto the interaction Hamiltonians in Eq. (10), As an example, in Fig. 2 we present the evolution of the eigenenergies E n of H M ES with the coupling energy, λ . For simplicity, we use RCs with M = 2, that is, the model includes a central spin and two two-level RCs. Overall, there are eight (many-body) occupation-basis states. For simplicity, the model is made to be fully symmetric; λ = λ i and Ω = Ω i .
Since we work with Ω > λ and Ω ∆, the global eigenstates of H M ES can be associated with the local picture of a central spin and two (left and right) RC oscillators. The lowest two eigenstates in Fig. 2(b), are separated by ≈ ∆. They correspond to the two RCs residing in their respective ground states (energy 2 × 1 2 Ω), while the central spin occupies either its ground or excited states. Higher eigenenergies in Fig. 2(a) correspond to excited RCs, and are separated by multiples of Ω. Specifically, there are four levels clustered around 2Ω corresponding to one RC in its ground state and the other in its first excited state, with the spin occupying its ground or excited state. The top two levels in Fig. 2 are associated with both RCs being excited, and the spin occupying its ground or excited state. One can similarly explain level structure for M > 2.
We now discuss our choice of parameters. In previous studies e.g. in Ref. 53 , the parameters of the RC model Eq.
(3) were such that Ω ∼ ∆, and Ω ∼ T . Here, with the objective to uncover the role of strong couplings, we use a high frequency RC mode, Ω ∆ and Ω T such that we can safely truncate the spectrum of each RC and include few levels when performing heat transport calculations. This choice of parameters allows us to focus on the strong-coupling limit while bypassing issues of convergence involved in the truncation of the harmonic oscillator manifold. As a result of this truncation, we are able to construct an effective model, which exposes the essence of strong interactions, and is highly accurate as we show in section IV. In Appendix D we include additional simulations of the heat current at strong coupling using smaller values of Ω, down to Ω = 5∆.
Inspecting Fig. 2, we make a useful observation: Since in simulations we work with T Ω, excited states beyond the first two levels are thermally inaccessible as there is a dramatic energy separation between the first two levels (gapped by ∆), and the rest of the eigenenergy manifold. Based on this observation, we truncate the spectrum of H D ES and write down an effective spin-boson Hamiltonian, with We denote the pair λ i with a vector notation, λ = (λ h , λ c ). E n ( λ ) are the eigenvalues of H D ES with eigenstates |n( λ ) ; n = 1, 2, as depicted e.g. in Fig. 2. We refer to the Hamiltonian (12) as the "effective spin-boson (EFF-SB) model".
We highlight the dependence of the spin splitting on the interaction energies λ . The dimensionless function f i ( λ ) arises from the transformation of the enlarged system Hamiltonian into the diagonal basis and it describes the dressing of the system-bath coupling energies. It depends on both coupling energies, λ h and λ c and this fact immediately suggests that cooperative, two-bath interactions should play a role in the dynamics and steady state transport behavior, as indeed is the case based on the polaronic approach 15,16,19,21 .
To understand the impact of strong system-bath couplings, we depict in Fig. 3 the effective spin splitting and the coupling function f i ( λ ). Since henceforth we use equal couplings (until Sec. V), we relieve the vector notation, and identify by λ the (equal) couplings of the spin to both reservoirs. On the one hand, strong coupling effects should facilitate transport by enhancing the coupling of the spin to the baths. Indeed, f (λ ) is an increasing function of λ . On the other hand, the detrimental impact of strong coupling clearly emerges through the suppression of ∆ e f f (λ ), the tunneling splitting between the levels of the effective spin. Roughly, as we reduce the spin splitting, less energy is transmitted per quanta transferred. In Appendix C we present analytical derivations on a simpler model to support the observed suppression of ∆ e f f with λ . When focusing on quantum dynamics, the suppression of ∆ e f f (λ ) would eventually lead to the localization of the spin 62 . In calculations of steady state heat transport as we perform here, the substantial suppression of ∆ e f f (λ ) shows up in the rapid decline of the heat current at strong coupling.
We fit the functions in Figure 3(a) and find that, approximately, λ 2 m and f (λ ) ∝ λ , with λ m a parameter with the dimension of energy. Importantly, since f (λ ) is of order 1, as long as γ 1 [see Eq. (5)], we can still safely employ a weak coupling scheme to treat the dynamics and transport characteristics of the effective-SB Hamiltonian Eq. (12).
We present the combination Fig. 3(b) and reveal a non-monotonic behavior: The product first grows approximately as λ 2 then it decays rapidly as exp −λ 2 /λ 2 m , with a maximum at λ m / √ 2. As we discuss in more details in Sec. IV, the steady state heat current in the SB model is proportional to this product for Ohmic baths under the weak system-bath coupling approximation 15,16,21,22 . Thus, the trend observed in Fig. 3(b) is representative of the heat current in this regime.
While the condition of weak system-bath coupling does not hold for the SSB model since λ /∆ can be large, it is valid for the effective SB model since its coupling to the bath is dictated by the width parameter γ i , which is assumed small. This is the main physical and computational outcome of the mapping approach: The effective spin-boson model captures strong coupling effects, as inherited from the SSB model, but now embedded in renormalized parameters while allowing a weak coupling treatment to the baths. The effective SB model described in this subsection is highly beneficial as it hands over analytic understanding on transport characteristics. It is important to appreciate that the EFF-SB model should be identified only after constructing a RC that converges with respect to M, the number of levels representing the RC. This aspect is clearly manifested in Fig. 3.
While ∆ e f f (λ ) and f (λ ) readily converge with M at weak-intermediate couplings, at strong coupling increasingly more levels should be included to converge the parameters of the effective SB model.

III. METHOD: STEADY STATE ENERGY TRANSPORT WITH THE RC-QME
We succeeded in transforming the SSB model with possibly strong coupling of the spin to harmonic oscillator baths into an effective SB model. The parameter γΩ in the SSB model, which controls the width of the spectral function, dictates the systembath coupling energy in the RC picture. Thus, if we assume a narrow spectral function for the SSB model, in the RC model the extended system becomes weakly coupled to its surroundings. As a result, the dynamics and steady state characteristics of the effective SB can be studied using a QME perturbative in the system-bath coupling.
Here, we employ the Markovian Redfield master equation to study the steady state behavior under the Hamiltonian (7)-(10). This approach relies on three assumptions: (i) The enlarged system is weakly coupled to its residual environments such that a second-order perturbative approach yields proper results. (ii) The residual environments act as Markovian baths. (iii) The initial state of the total system factorizes into a system times bath form. The residual thermal baths are individually prepared in their canonical states characterized by inverse temperatures β i .
Working in the Schrödinger picture and in the energy basis of the subsystem, the Redfield equation for the reduced density matrix of the enlarged system (dimension 2M 2 ) is given by 2,3 , Here, ω mn = E m − E n are the bohr frequencies of the enlarged system in the energy basis. The R super-operator terms in Eq. (14) are given by half Fourier transform of bath autocorrelation functions, The autocorrelation functions of bath operators are evaluated with respect to the thermal state of their respective baths. In what follows, we neglect the imaginary component of the dissipators, ∆ RC , as it can be shown that they cancel out the quadratic term in the RC Hamiltonian, Eq.(3) 53 .
For harmonic environments and a bilinear coupling between the enlarged system and the residual environments, the real part of the R term reduces to Here, n i (ω) is the Bose-Einstein distribution function of the ith bath; J RC,i (ω) is the RC spectral function given by Eq. (5).
To calculate the heat current flowing through the system, we write down the Redfield equation aṡ with the dissipators D h,c (ρ ES ) organized based on Eq. (14). The dissipators are additive given the weak coupling of the RCs to their residual baths (but each dissipator D i depends non-additively on the original coupling parameters, λ i ).
We solve the equation of motion in steady state,ρ ES (t) = 0 and obtain the density matrix of the enlarged system, ρ ss ES . The steady state heat current at the ith bath, calculated from the heat exchange between the enlarged system and this reservoir is given by The heat current is defined positive when flowing from the ith bath towards the system. This procedure, the application of the Redfield-QME onto the RC model is referred to as the reaction-coordinate quantum master equation method, RC-QME. Below, we perform heat transport simulations of the nonequilibrium spin-boson model using four methods: (i) RC-QME framework: We follow the approach described in this Section, while checking convergence by increasing M, the number of levels used to describe the RCs.
(ii) Weak coupling, BMR-QME calculations: We calculate the heat current with the original SSB Hamiltonian Eq. (1) without performing the RC mapping under the (unjustified) weak coupling assumption, so as to highlight signatures of strongcouplings. In this case, we again employ the Redfield Equation Eq. (14) in the diagonal representation but use J RC (ω) − → J SSB (ω). Obviously, this calculation (based on an analytic expression 15,16,21,22 ) is appropriate only when the coupling of the spin to the baths is weak.
(iii) PT-NEGF: We use the polaron-transformed NEGF method developed in Ref. 30 , denoted by "PT-NEGF". This method, tested on the nonequilibrium spin-boson model, is capable of treating nonperturbatively the system-bath coupling term. It also handles the tunneling energy in a perturbative manner. The method combines the polaron transformation with the nonequilibrium Green's function method, where a Majorana-fermion representation is utilized to evaluate the Dyson series.
(iv) EFF-SB model: We obtain the parameters of the effective two-state model based on the diagonalization of H M ES with large enough M. However, since we end up with two states weakly coupled to the surroundings, we can utilize the closed-form expression for the heat current as derived in Refs. 15,16 , assuming a Markovian, Ohmic bath. We write down this expression emphasizing the nontrivial dependence of parameters on the renormalized spin splitting, The Ohmic spectral function, which is evaluated at the renormalized splitting, is dressed by f (λ ). For brevity, we suppress the dependence of ∆ e f f on λ , For a symmetric system in the limit of high Λ the heat current reduces to In Fig. 3(b) we showed that the temperature-independent prefactor of this function displayed a turnover behavior with λ ,

A. System-bath coupling
We begin by studying the behavior of the heat current as a function of the system-bath coupling strength. In the SSB model of Sec. II A, the parameter λ i quantifies the system-bath coupling, see Eq.
(2). In the RC picture, this parameter describes the coupling of the spin to the added reaction coordinates, while the enlarged-system-bath coupling is controlled by γ. Figure 4 depicts the steady state heat current as a function of λ using our RC-QME framework. We further show the BMR-QME, EFF-SB, and the PT-NEGF results. For simplicity, we assume a symmetric setup, λ = λ i , Ω = Ω i , γ = γ i , Λ = Λ i .
We find that the heat current obtained from the BMR-QME grows quadratically with λ ; recall that the spectral function scales as J SSB (ω) ∝ λ 2 . In contrast, the RC-QME depicts a turnover behavior, in a good agreement with the PT-NEGF method (which takes into account strong coupling effects through the Polaron transformation and the Dyson-series summation). Notably, the turnover behavior of the current with the coupling energy is excellently captured by the effective-SB model, Eq. (19) with the renormalized parameters, displayed in Fig. 3. Parameters for the effective model were taken using M = 4. Figure 5 shows the steady state heat current for a smaller value of Ω. Recall, that this parameter controls the peak position and width of the SSB spectral density, Eq. (2). We again observe a very good agreement between strong coupling methods. Notably, we notice that the EFF-SB better agrees with the PT-NEGF technique than the RC-QME method at M = 4. This trend is further exemplified in Appendix D, and discussed in details in Sec. IV D.

B. Tunneling splitting
We study the behavior of the heat current as a function of the tunneling splitting ∆ in Fig. 6, separately looking at weak (a), intermediate (b)-(c) and strong (d) system-bath couplings. The dependence of J q on the tunneling splitting was previously investigated analytically with NIBA 19 and the BMR-QME 20 , and numerically, with exact methods 20,32,33 . The general trend observed is of the heat current first growing with ∆ quadratically. However, for larger spin splitting such that ∆ > T , J q begins to decay with ∆. This suppression of the current can be rationalized in the weak coupling picture when only resonant processes contribute: When ∆ becomes large, bath modes at this frequency are not thermally populated, resulting in the suppression of the heat current.
Inspecting Fig. 6, we find that the different methods agree at weak coupling. However, as we move into the intermediate and strong coupling regimes the current predicted by the BMR-QME offsets from the RC-QME. This is expected: We learned from Fig. 3 that strong coupling manifests itself in the renormalization (suppression) of the tunneling splitting. As such, in the  BMR-QME method the position of the peak occurs at ∆ ∼ T while it appears only at a higher splitting for the RC-QME, once ∆ e f f ∼ T . Interestingly, while the BMR-QME method incorrectly positions the peak, the maximal magnitude of the current is the same as in the RC-QME method. As for the PT-NEGF method, while it agrees with the RC-QME at weak-intermediate couplings, see Fig. 6(a)-(c), more substantial deviations show up at strong coupling (up to 50%) and for large splitting, see Fig. 6(d). However, while the magnitude of the current differs between the RC-QME and the PT-NEGF methods, the peak position agrees. This manifests that the essential aspect of gap renormalization is captured in both methods, which is significant given the distinct ways they handle strong coupling effects. Which method, RC-QME or PT-NEGF is more accurate at strong coupling and large spinsplitting? We think that the RC-QME provides a more accurate result in this difficult regime: The PT-NEGF method does not handle exactly the spin splitting at large ∆; while it is more accurate than the Markovian NIBA method of Ref. 19 , which can only capture the J q ∝ ∆ 2 scaling, the self-energy in the PT-NEGF includes ∆ only perturbatively. In contrast, the RC-QME handles both ∆ and λ non-perturbatively. Nevertheless, additional benchmarking is required to determine the exact behavior of the current at strong coupling and at large tunneling splitting.

C. Temperature-difference dependence
We study the behavior of the heat current as a function of the average temperature, T a ≡ (T h +T c )/2 in Figs. 7-8. Inspecting the low-intermediate temperature regime of T a < 2∆ in Fig. 7, we observe that the different methods display a crossover behavior: The current rises with T a as it approaches a resonance condition at ∆ (or ∆ e f f in the RC-QME). However, at higher temperatures the current drops as 1/T a due to the temperature dependence of the scattering rates as captured by Eq. (19). The BMR-QME overestimates the current at strong couplings since it does not capture the suppression of ∆, thus it allows for the transfer of more heat than in the RC-QME method.
As we raise the temperature, T a > 2∆, we observe in Fig. 7, and in more details in Fig. 8, a dramatic enhancement of the  Fig. 7 (c), but continuing here into the high temperature regime. Parameters are the same as in Fig. 4. current, particularly at strong coupling. As a clue to deciphering the origin of this enhancement, we note that it is not captured by the EFF-SB model thus it emerges from the contributions of higher-energy transitions in the manifold of the extended system, see Fig. 2. Naturally, we question whether the enhancement of the current at high temperature is physical, or is it an artifact of the extraction of a single oscillator mode as a RC. However, in Fig. 8 we show that the PT-NEGF method also displays this enhancement, though there are differences in the magnitude of the current.
Altogether, as a function of temperature, we gather from Figs. 7 and 8 that the heat current is peaked first around T a ≈ ∆ then at T a ≈ Ω (there are additional peaks at higher temperature). In the RC picture, we can explain this structure as follows: The current is high when T a roughly matches transition frequency in the extended system. The first peak corresponds to the excitation of the spin, with the RCs in their ground state. The second peak corresponds to energy transfer via the excitation of the RCs.
Interestingly, using the PT-NEGF method we can explain these peaks with a different picture, thinking about resonant and off-resonant (tunneling) transport of phonons. The first peak at T a ≈ ∆ corresponds to resonant transmission of phonons from the bath through the spin. Beyond that at higher temperatures, the transport process occurs increasingly via tunneling: high frequency phonons from the baths must utilize the low frequency spin impurity to transfer heat. This off-resonant process becomes more potent as the temperature rises eventually matching the peak of the baths spectral density function (recall that we use a Brownian function in the SSB model), where the bath density of state is maximized. If one continues and increases the temperature, additional peaks will show up at multiples of Ω. They would represent multi-phonon transfer processes in the PT-NEGF method, and the excitation of high energy levels in the extended-system RC model.  7) as long as T Ω, ∆ Ω. What is the explanation for the remarkable success of the EFF-SB model? We can explain this superior behavior by noting that weak-coupling methods typically overestimate the heat current compared to exact results 20 . In the RC-QME framework, as we include a large number of states for the RC (M > 2), many transitions contribute to the heat current, thereby increasing the error associated with the weak coupling scheme. Thus, while a large M more accurately represents the RC oscillators, at the same time errors associated with the calculation of the current using a weak-coupling method become more significant.
In contrast, the effective SB model offers an excellent, fortunate tradeoff of these two aspects: On the one hand, it is constructed by diagonalizing a large M-RC Hamiltonian, Eq. (7); the renormalized parameters therefore capture the RC manifold. On the other hand, the subsequent truncation of the 2M 2 spectrum, to include only few (two) states limit the over-estimation error associated with weak coupling schemes.

V. THERMAL-DIODE EFFECT
The nonequilibrium spin-boson model was originally introduced to describe a quantum thermal diode (rectifier) 15,16 . In a boson-bath nanojunction, this effect relies on two conditions: anharmonicity, captured here by the impurity spin, and spatial asymmetry, introduced here by coupling the spin asymmetrically to the different baths. Recently, the nonequilibrium spin-boson model was realized in a superconducting quantum circuit demonstrating a diode effect with rectification ratio R = | 17 . An important question that has been touched on, without conclusion is whether strong coupling effects are beneficial for the diode effect, and to what extent. The EFF-SB model offers computation at low cost, and we utilize it here to address this question.
The quantity of interest here is the thermal rectification ratio, R, which quantifies transport asymmetry when reversing the temperature bias. In this definition, R = 1 corresponds to lack of rectification, while R > 1, or R < 1 corresponds to a diode effect. To materialize this effect, spatial symmetry needs to be introduced into the model. This is done here through the use of an asymmetry parameter −1 < χ < 1, with λ h = λ (1 − χ), λ c = λ (1 + χ). In this definition, χ = 0 corresponds to the symmetric model where no rectification is expected.
Using Eq. 19, we obtain an analytic expression for the rectification ratio This equation reveals that the ratio of the effective coupling parameters breaks spatial symmetry in the system, allowing a diode effect to emerge. Indeed, if f L ( λ ) = f R ( λ ), then R = 1 and the system will not exhibit a diode behaviour. Figure 9 examines the effect of asymmetry in the coupling strength on the diode effect. We test the effect at different coupling strengths λ and for two temperature biases. First, in Fig. 9 (a) and (c) we display the heat current itself in the forward direction using the BMR-QME and the EFF-SB methods. Notably, the current is not maximal when the coupling is symmetric, rather, the current reaches an extremum when the qubit couples more strongly to the cold bath. Fig. 9(b) and (d) display the corresponding rectification ratios using the BMR-QME and the effective SB models. At stronger coupling the rectification ratio increases. However, this increase is rather minor. Identifying strategies for significantly enhancing the thermal diode effect in impurity models remain a challenge.

VI. SUMMARY
We presented the reaction coordinate framework, RC-QME for studies of quantum heat transport in nanoscale systems beyond the weak system-bath coupling limit. In this method, degrees of freedom of the baths, which are strongly coupled to the system are extracted and explicitly included into the system. After performing additional, physically-motivated transformations and truncation, transport characteristics of the enlarged system are computed using a perturbative, Markovian QME method. The main results of our work are: (i) On the fundamental side, by continuing with physically-motivated transformations of the RC Hamiltonian we ended up with a spin-boson model, albeit with renormalized parameters. The effective SB model elegantly reveals the impact of strong system-bath coupling on spin dynamics and heat transport behavior, by manifesting competing effects: When increasing the system-bath coupling energy the effective spin splitting is strongly reduced. On the other hand, the coupling of the effective spin to the baths increased approximately linearly. One could generalize the principles behind the effective SB model, construct other effective few-level models with energy parameters renormalized to absorb strong coupling effects, and describe the operation of more complex thermal machines. (ii) On the application side, by comparing our method to another scheme that utilizes the polaron transformation, we can safely conclude that our method captures system-bath coupling effects beyond low order perturbation theory treatments. Specifically, the RC-QME method reproduced the turnover behavior of the heat current with λ , the system-bath interaction parameter, in a quantitative agreement with the polaron-NEGF method. Furthermore, we demonstrated that the RC-QME method took into account off-resonant phonon transmission processes. This was manifested in the behavior of the heat current with temperature: The current showed multiple peaks representing resonant, then enhanced off-resonant transport once the temperature matched the peak of the bath density of states (spectral function). We further looked at the diode behavior of the asymmetric SB model and showed that strong coupling effects could enhance the rectification behavior, albeit to a small extent.
The construction of an effective SB model is one of the central results of this work, allowing us to rationalize strong coupling effects, as well as perform calculations at minimal computational cost. Once reaching the effective SB model, quantum dynamics and transport behavior can be evaluated with a perturbative QME method, as we did here. Moreover, one can further employ the effective SB model in conjunction with more accurate techniques developed for the SB model, such as the NIBA 19 and Majorana-fermion NEGF methods 29,30 . Even more so, once reaching the effective-SB model one could continue iteratively with the RC transformation, to extract additional modes into the system (transforming the spectral function numerically). As such, this technique is not limited to spectral functions that are peaked around a central mode. Our method differs from chain mapping techniques such as TEPODA [44][45][46] , since we are able to study the strong coupling regime semi-analytically, gaining physical intuition and minimizing computational costs.
The RC-QME method described in this paper should be extended in two ways: First, one should test whether it is still accurate in models where the spectral density function of the SSB is not peaked at a single frequency as in the Brownian oscillator model, but rather resembles e.g. an Ohmic model. In such cases, one should find how many oscillators should be extracted from the baths to capture quantitatively strong coupling effects, and how does the reaction coordinate method correspond to the polaron model. Second, the RC method and the effective Hamiltonian principle can be feasibly generalized to describe the coupling of multi-level quantum systems to multiple reservoirs. As such, it is suitable to study the operation of strongly-coupled quantum heat machines, our objective in future studies.
generalize the discussion, we consider a coordinate q moving in a confining potential U(q). We will make use of classical equations of motions but an analogous treatment can be done quantum mechanically. This can be done because the spectral density function only contains information about the harmonic baths and their interaction with the system, which is linear in the oscillators' displacements. Our discussion here follows Ref. 53 . It is included here for completeness, demonstrating its generalization for the case of a system coupled to two heat baths.
Starting from the SSB Hamiltonian Eq. (1), the classical analogue of the SSB model takes the form Here, we introduce the compact form φ k,i ≡ 2ν k,i f k,i . Also, we identify the position and momenta of the environment as: and Using this Hamiltonian, we obtain the classical equations of motion for the coordinate q and each bath mode x k,i : Here, U (q) corresponds to a derivative with respect to the coordinate q. Working in Fourier space, where observables take the formq(z) = ∞ −∞ dte −izt q(t), it is possible to eliminate the environment variables from the equations of motion. Doing so, we obtain an algebraic equation K(z)q(z) = −Ǔ (z). Here, the Fourier space kernel K(z) is defined as: The integral can be solved by the residue theorem making use of the simple pole at ν = z. thus, the kernel takes the form: At this stage, we repeat the same procedure after the RC transformation, and identify the kernel K(z) in terms of the RC spectral densities. Since the mapping is exact, the two kernel functions, before and after the normal mode transformation, should be identical. Following a similar procedure as before, we replace the Hamiltonian Eq. (3) by a classical coordinate moving in a potential U(q). This yields where we define scaled coefficients κ i = √ 2Ω i λ i and ξ k,i = 2g k,i Ω i ω k,i , while x i , p i and X k,i , P k,i are positions and momenta of the RC and the environments, respectively. The Hamiltonian can be compactly written as H s (q, x i ) includes the terms from the first line of Eq. (A8). This Hamiltonian leads to equations of motion of the form Working in Fourier space, it is possible to eliminate both the RC and environment variables from the equation of motion such that the kernel takes the form: where, using the definition of the RC spectral density: If we assume that the RC spectral density is Ohmic, J RC,i = γ i ωe −|ω|/Λ i , where Λ i is a high frequency cut-off, then L i (z) = −z 2 + 2iπzΩ i γ i . By plugging this expression into Eq. (A13), writing z = ω − iε, we equate the two expressions for the kernel K(z). This leads to the expression, which is the associated spectral density function for the "standard" spin-boson model.

Appendix B: Rotated Reaction Coordinate Model
In this Appendix, we perform a rotation transformation on the RC Hamiltonian, Eq. (3), and absorb the quadratic term. This transformation is exact. There are two terms, one for each RC, which are quadratic in the displacement of the reaction coordinate . These terms can be absorbed by renormalizing the corresponding reaction coordinates. To do this, we define a dimensionless parameter η i such that The last equality holds assuming Ohmic spectral functions for the residual baths. We now perform a rotation transformation of the RC oscillators,ã i = a i cosh r i + a † i sinh r i , a † i = a † i cosh r i + a i sinh r i , whereλ i = λ i (1 + 4η i ) 1/4 , with η i = γ i Λ i Ω i . The Hamiltonian in Eq. (B4) includes a spin coupled to two "primary" harmonic oscillators, each coupled to secondary harmonic baths. We select the spectral functions of the residual baths in the rotated-RC basis to be Ohmic, and find after a long algebra that it reduces to The relationship between the SSB model and the rotated-RC model is a useful result of this work as it is an exact mapping, however, it comes with a warning. The behavior of the rotated-RC Hamiltonian cannot be studied with a Markovian QME method, as we discuss in Sec. II B. Other techniques, which take into account environmental memory and strong correlation effects should be used for an accurate treatment of this model. quantum master equation approach becomes less accurate when the width parameter grows. This explains the large overshoot of the current. In order to accurately simulate the current in this regime, one would probably need to extract several modes ("chain mapping") from the reservoirs by iteratively performing the RC transformation.
Remarkably, the EFF-SB treatment successfully captures the heat transport characteristics, compared to the PT-NEGF method, at a much cheaper computational cost, performing better than the M level RC method. We discuss the superior performance of the EFF-SB method in Sec. IV D.