Heat transport via a local two-state system near thermal equilibrium

Heat transport in spin-boson systems near the thermal equilibrium is systematically investigated. An asymptotically exact expression for the thermal conductance in a low-temperature regime wherein transport is described via a co-tunneling mechanism is derived. This formula predicts the power-law temperature dependence of thermal conductance $\propto T^{2s+1}$ for a thermal environment of spectral density with the exponent $s$. An accurate numerical simulation is performed using the quantum Monte Carlo method, and these predictions are confirmed for arbitrary thermal baths. Our numerical calculation classifies the transport mechanism, and shows that the noninteracting-blip approximation quantitatively describes thermal conductance in the incoherent transport regime.


Introduction
Heat transport via small systems has recently attracted considerable attention because a lot of intriguing phenomena can emerge reflected from the properties of a system and the surrounding environment. For instance, quantized thermal conductances have been observed in heat transport by phonons [1,2] and photons [3] in a manner similar to electric transport [4]. Thermal rectification [5,6] and thermal transistors [7] have also been theoretically proposed in analogy to electronic devices. Heat transport via a quasi one-dimensional material, e.g., carbon nanotubes, shows neither diffusive nor ballistic transport, which is currently categorized as anomalous transport [8]. Heat transport due to magnetic excitation is now a key ingredient in the field of spintronics [9]. Studying the general properties of thermal transport using typical systems is clearly an important subject not only for theoretical development but also for future experiments.
The spin-boson system is one of most common and important systems for describing a local discrete-level system embedded in a bosonic thermal environment [10,11]. This system has numerous applications, e.g., it is used to describe molecular junctions [12], superconducting circuits [13], and photonic waveguides with local two-level systems [14]. Hence, it is regarded as a minimal model for describing a zero-dimensional object with discrete quantum levels surrounded by a bosonic environment. One of the important problems here is to clarify the dissipative dynamics of the system near the equilibrium situation [10]. Depending on the properties of the thermal environment, the behavior of the autocorrelation function of the system changes from coherent oscillation to incoherent decay as a function of time. Intriguingly, at zero temperature, a quantum phase transition occurs when the coupling strength between the system and the environment is changed [15,16]. The sub-ohmic environment induces a second-order phase transition [17,18,19,20,21,22], while the ohmic case shows a Kostelitz-Thoulesstype phase transition [11,23,24]. The super-ohmic case does not have a distinct phase transition but exhibits a crossover. In addition, the ohmic environment induces the Kondo effect [25] at sufficiently low temperatures [10,11,26,27]. From this background in an equilibrium situation, it is quite natural to ask what happens if one considers heat transport in this system. Herein, we present systematic studies of heat transport via the spin-boson system and derive some exact results for this case.
A number of studies have investigated heat transport via spin-boson systems [28,29,30,31,32,33,34,35,36]. Segal et al. introduced an iterative path-integral technique for numerical calculations to investigate the far-from-equilibrium regime [28]. Ruokola and Ojanen studied low-temperature properties using a perturbation method and discussed co-tunneling mechanisms [29]. However, their methods do not seem to succeed in reproducing low-temperature properties, e.g., the Kondo effect. Two of the present authors (TK and KS) have focused on the transport properties in an ohmic environment and found several Kondo signatures [30], including the T 3 -temperature dependence of the thermal conductance. Herein, we advance in this direction and cover arbitrary types of environments. We consider a general picture for understanding the transport properties at extremely low temperatures for the whole regime of spectral densities and quantitatively characterize the transport mechanism for all temperature regimes.
We present our findings in this paper to distinguish them from existing literature. First, we derived an asymptotically exact expression for the thermal conductance in the extremely low-temperature regime, reproducing the aforementioned T 3 -temperature dependence of thermal conductance in the ohmic case. Our formula is asymptotically exact in the co-tunneling transport regime and predicts power-law temperature dependences ∝ T 2s+1 for the thermal environment of spectral density with the exponent s. Second, we performed accurate numerical calculations to investigate thermal conductance over the entire temperature regime. We confirmed the temperature dependencies predicted by our expressions for the co-tunneling and the sequential tunneling transport regimes. Furthermore, we found that the noninteracting-blip approximation (NIBA) [10] describes thermal conductance in the incoherent tunneling regime accurately. In table 1 the transport mechanisms for each regime are summarized and relevant analytical descriptions are presented. In the table, sequential tunneling, co- Table 1. Summary of the relevant transport process. Here, ∆ eff is an effective tunneling amplitude [see equations (23) and (24)] and T * is the crossover temperature [see equation (31)]. The last column shows the temperature dependences of the thermal conductance, where "Schottky" indicates a Schottky-type temperature dependence proportional to e − ∆ eff /kBT /T 2 . The temperature dependence of NIBA is complex in general, and the symbol ( * ) indicates the high-temperature limit.

Exponent
Condition Transport process Dependence

Sequential tunneling Schottky
tunneling, and NIBA imply the analytical descriptions based on the approximate form [equation (29)], the asymptotically exact expression [equation (33)], and the analytical descriptions based on the NIBA expression[equation (22)] with equations (41) and (42). The paper is organized as follows. In section 2, we introduce the model and explain the Meir-Wingreen-Landauer-type formula. In section 3, we classify the transport mechanism and derive an asymptotically exact expression that is valid in the cotunneling transport regime. We perform numerical calculation using the quantum Monte Carlo method, and compare the results with analytic approximations in section 4. In section 5, we summarize our work.

Model
We consider heat transport via a local quantum system coupled to two reservoirs denoted by L and R. The model Hamiltonian is given by Figure 1. Symmetric double-well potential of the local system. An energy spacing of quantum levels in each well is ω 0 (indicated by the blue sold lines), and an energy splitting due to quantum tunneling (indicated by the red dashed lines) is ∆ = E e −E g , where E g and E e are the ground-state energy and the first-excited-state energy, respectively.
where H S , H ν , and H I,ν describe the local system, the reservoir ν (= L, R), and the interaction between them, respectively. The operators p and x are the momentum and position for the local system, respectively, and V (x) is the potential energy. The reservoirs comprise multiple phonon (or photon) modes, which are described in general by harmonic oscillators with frequency ω νk and mass m νk , where the subscript denotes the phonon (photon) wavenumber k in the reservoir ν. The momentum and position of an individual oscillator are denoted by p νk and x νk , respectively. For simplicity, the system-reservoir coupling H I,ν is consider as a bilinear form of x and x νk , and the interaction strength is denoted by C νk . The second term of H I,ν is a counter term to cancel the potential renormalization due to the reservoirs. In this study, the potential energy V (x) of the local system is considered as a doublewell potential as shown in figure 1. We assume that the barrier height of the double-well potential is sufficiently large in comparison with ω 0 , where ω 0 is the frequency of a small oscillation at the potential minima x = ±x 0 /2. Then, quantum tunneling between the two wells induces small energy splitting ∆ ( ω 0 ) between the ground-state energy E g and the first excited energy E e .
After truncating the local system into two states by considering the two lowest energy eigenstates, we obtain the spin-boson model [11,10]: Here, σ i (i = x, y, z) is the Pauli matrix, b νk is an annihilation operator defined by and λ νk = x 0 C νk / √ 2 m νk ω νk . In the present model, we assign the localized states at the left (right) well as |↓ (|↑ ). Throughout this study, we examine the symmetric doublewell potential (ε = 0) and only use the bias term εσ z to define the static susceptibility where · · · implies an equilibrium average. For the symmetric case (ε = 0), the system Hamiltonian H S describes the tunneling splitting ∆ between the ground state (σ x = +1) and the first excited state (σ x = −1). The properties of the reservoirs are characterized by the spectral function which is considered to be continuous assuming that the number of phonon (photon) modes is large. For simplicity, we assume the following simple for the spectral function [11,10]: where α ν is the dimensionless coupling strength between the two-state system and the reservoir ν. To cut off high-frequency excitation, we introduced the exponential cutoff function e −ω/ωc , where ω c is the cutoff frequency, which is considerably larger than other characteristic frequencies, e.g., ∆, ε/ , and k B T / . The exponent s in equation (13) is crucial for determining the properties of the reservoirs. The case s = 1 is called "ohmic," whereas the cases s > 1 and s < 1 are called "super-ohmic" and "sub-ohmic," respectively.

Thermal conductance
The heat current flowing from reservoir ν into the local two-state system is defined as follows: Using the standard technique of the Keldysh formalism [37,38,39], one can derive the Meir-Wingreen-Landauer-type formula [40] for the nonequilibrium steady-state heat current J L = − J R ≡ J as follows [7,6,30,41]: where α = α L + α R , γ = 4α L α R /α 2 is an asymmetric factor, n ν (ω) is the Bose distribution function in reservoir ν, and χ(ω) is the dynamical susceptibility of the two-state system defined by Equation (15) is derived in Appendix A. The linear thermal conductance is defined as Using the exact formula [equation (15)], the linear thermal conductance is given as where χ(ω) is evaluated for the thermal equilibrium and β = 1/(k B T ). Thus, we need to calculate the dynamical susceptibility χ(ω) for evaluating the linear thermal conductance.
For convenience of discussion, we also introduce a symmetrized correlation function and its Fourier transformation: From the fluctuation-dissipation theorem [10], the imaginary part of the dynamical susceptibility is related to S(ω) as The thermal conductance is then rewritten using the correlation function S(ω) as

Classification of Transport Processes
The dynamics of dissipative two-state systems have long been studied using a number of approximations [11,10]. In this section, we re-examine such analytic approximations from the viewpoint of heat transport. In section 3.1, we first consider the effective tunneling amplitude and discuss a quantum phase transition driven by strong systemreservoir coupling. Next, we consider the three mechanisms, which we call "sequential tunneling" (section 3.2), "co-tunneling" (section 3.3), and "incoherent tunneling" (section 3.4) following in the previous literatures [11,10,29,42]. We derive analytic expressions for the thermal conductance in each transport process. We also introduce NIBA in section 3.5.
In this section, we show two novel results of our study. The first concerns the co-tunneling process. We derive an asymptotically exact formula for the co-tunneling (a) (b) Figure 2. Schematics of the ground-state wavefunction (a) below the transition (0 ≤ α < α c ) and (b) above the transition (α c < α). The former state is delocalized, whereas the latter is localized at one of the two wells. For the localized state, quantum tunneling between the two wells is forbidden since the overlap integral between the states in the two wells vanishes.
process by utilizing the generalized Shiba relation. This formula always holds at low temperatures for an arbitrary exponent (s) as long as the ground state of the system is delocalized. The second result is related to the incoherent tunneling. In particular, we find that the Markov approximation is inadequate to describe the thermal conductance in the incoherent tunneling regime. Instead, the thermal conductance in this regime is well described by NIBA, which considers the non-Markovian properties of stochastic dynamics. We show that NIBA quantitatively explains numerical calculations in section 4.

Effective tunneling amplitude and quantum phase transition
One important effect of the system-reservoir coupling is renormalization of the tunneling amplitude ∆. In this subsection, we briefly show the effective tunneling amplitude results obtained via adiabatic renormalization [11,10]. A detailed derivation is given in Appendix B.
For the ohmic case (s = 1), the effective tunneling amplitude is given by This result indicates a phase transition at zero temperature, for which the critical value of the system-reservoir coupling is α = 1 [16,15]. For system-reservoir couplings below the transition (0 ≤ α < 1), the ground state is non-degenerate, as shown in figure 2 (a), indicating the coherent superposition of the two localized states |↑ and |↓ . We call this ground state "delocalized." For strong system-reservoir couplings above the transition (α > 1), the coherent superposition of the two localized states is completely broken, leading to the doubly-degenerate ground states shown in figure 2 (b). We call this ground state "localized." In this localized regime, quantum tunneling between the wells is forbidden at zero temperature since there is no mixing (∆ eff = 0) between the two localized states. Thus, the present quantum phase transition can be recognized as a "localization" transition that separates the delocalized and localized regimes at zero Figure 3. Schematic of the sequential tunneling process. Heat transport occurs by a combination of (a) phonon (photon) absorption and (b) phonon (photon) emission.

temperature.
For the sub-ohmic case (s < 1), the adiabatic renormalization always leads to an effective tunneling amplitude of zero (∆ eff = 0). This is correct in the limit ∆/ω c → 0, as discussed in a previous study [11]. However, for a finite value of ∆/ω c , the naive adiabatic renormalization procedure yields incorrect results and should be improved. In subsequent theoretical studies [43,44], it was found that the localization transition actually occurred at a critical system-reservoir coupling (α = α c ), where the critical value α c depended on both s and ∆/ω c . The existence of the localization transition was also confirmed via numerical calculations [18,19]. In summary, for the sub-ohmic case, the ground state is delocalized for 0 ≤ α < α c , as shown in figure 2 (a), and localized for α c < α, as shown in figure 2 For the super-ohmic case (s > 1), the effective tunneling amplitude is always finite: where Γ(z) is the Gamma function. Therefore, there is no localization transition and the ground state is always delocalized, as shown in figure 2 (a).

Sequential tunneling
For weak system-reservoir couplings (α 1), the system and the reservoirs are almost decoupled and the interaction Hamiltonian H I,ν can be regarded as a perturbation. For the second-order perturbation, the system dynamics are described by a stochastic transition between the ground state (σ x = +1) and the excited state (σ x = −1), as shown in figure 3. The transition from the ground state to the excited state involves phonon (photon) absorption, and the inverse transition involves phonon (photon) emission. A combination of these two processes induces heat transport. We refer to this type of transport process as "sequential tunneling" by analogy with the electronic transport process through quantum dots. The transition rates for the process of phonon (photon) absorption and emission are calculated based on Fermi's golden rule as follows [11]: where I(ω) = I L (ω) + I R (ω) and n B (ω) = (e βω − 1) −1 is a Bose distribution function. Using these transition rates, the stochastic dynamics of the system are described using the Lindblad equation where ρ(t) is a density matrix of the system, L e = σ + x ≡ (σ z − iσ y )/2, and L a = σ − x ≡ (σ z + iσ y )/2. By solving this equation, we obtain the symmetrized correlation function as where Γ = (Γ e + Γ a )/2. The correlation function S(ω) has two peaks at ω = ±∆, reflecting the coherent system dynamics. Because Γ ∆ always holds in the weakcoupling regime, the correlation function is approximated as where δ(x) is the delta function. The thermal conductance for the weak coupling regime is obtained by substituting equation (28) into equation (22) as follows: This result is identical to the formula derived in previous research [5] and [30] using the master equation approach and is consistent with the perturbation theory [29]. For actual comparison with the numerical simulation in section 4, we improve the approximation by replacing ∆ with ∆ eff using adiabatic renormalization (see section 3.1).
The formula for sequential tunneling [equation (29)] is valid when For the sub-ohmic case (s < 1), this condition is never satisfied, indicating the absence of a sequential tunneling regime. For the ohmic case (s = 1), the condition is equivalent to α 1, whereas for the super-ohmic case (s > 1), the condition is always satisfied for a moderate temperature (k B T ∼ ∆ eff ). At high temperatures (k B T ∆ eff ), the condition is always satisfied for s ≥ 2, whereas for 1 < s < 2, it becomes where T * is the crossover temperature.
The formula for sequential tunneling [equation (29)] predicts the exponential decrease in the thermal conductance as the temperature is lowered. At low temperatures, the thermal conductance behaves as κ ∝ e − ∆ eff /k B T /T 2 ; this is because the transition from the ground state to the excited state is strongly suppressed if the thermal fluctuation is smaller than the effective energy splitting, i.e., when k B T ∆ eff . When the sequential tunneling process is strongly suppressed at low temperatures, equation (29) becomes invalid since another process becomes dominant, as discussed in the next subsection.

Reservoir L
Reservoir R Figure 4. Schematic of the co-tunneling process. At k B T ∆ eff , heat transport via a virtual excitation in the local system is dominant.

Co-tunneling and an asymptotically exact formula
At low temperatures, heat transport via the virtual excitation of the local two-state system becomes dominant (see figure 4); this transport process is known as "cotunneling" by analogy with the electronic transport process through quantum dots. In a previous study [29], an analytical expression for thermal conductance was derived using the fourth-order perturbation theory with respect to the interaction H I,ν . However, in this calculation the renormalization of the tunneling amplitude at a low temperature has not been considered.
Here, we derive a new asymptotically exact formula for the thermal conductance without any approximations. For this purpose, we focus on an asymptotically exact relation called the generalized Shiba relation [45,46]: where χ 0 is the static susceptibility defined in equation (10). This exact relation holds at low temperatures (k B T ∆ eff ) for arbitrary environments and arbitrary systemreservoir couplings. At low temperatures (k B T ∆ eff ), the dominant contribution to the integral of equation (18) comes from the low-frequency part (0 ≤ ω k B T ∆ eff ) due to the factor of the Bose distribution function. By substituting the low-frequency asymptotic form S(ω) πα( χ 0 /2) 2Ĩ (ω) into equation (18), we obtain This expression is similar to the co-tunneling formula in previous studies [29,47,42] but significantly differs in terms of static susceptibility, χ 0 , which considers higher-order processes. Equation (33) can be rewritten as Reservoir L Reservoir R Figure 5. Schematic of the incoherent tunneling process. The wavefunction is localized in the two wells, and a stochastic transition occurs between them.
where F (s) is a dimensionless function of s. Thus, we find that the thermal conductance κ is proportional to T 2s+1 at low temperatures. The same temperature dependence has been derived by the perturbation theory [29,47,42]. However, the perturbation theory cannot treat renormalization effect due to higher-order processes on the static susceptibility, and fails in predicting a correct prefactor including χ 0 . In contrast, the present result given in equation (33) is asymptotically exact, incorporating the renormalization effect appropriately. The co-tunneling formula [equation (33)], a new formula that is first derived in the present study, holds universally at low temperatures for an arbitrary exponent, s, as long as the ground state of the system is delocalized (∆ eff > 0) In a previous study [30], the thermal conductance in the ohmic case (s = 1) was shown to be proportional to T 3 , which is consistent with equation (33), and this T 3 -dependence was discussed in terms of the emergence of the Kondo effect. However, it is worth nothing that the power-law temperature dependences are derived in an unified way even in non-ohmic cases. These temperature dependences result from nontrivial many-body effects due to strong mixing between the system and the reservoirs.

Incoherent tunneling: the Markov approximation
For a strong reservoir-system coupling, the coherent superposition of the two localized states is completely broken. In such a situation, heat transport is induced by stochastic dynamics between the two localized states |↑ and |↓ , as shown in figure 5. We call this transport process "incoherent tunneling." Within the Markov approximation [48,49,50], the stochastic dynamics of the system are described by the master equation where P L (t) and P R (t) (= 1 − P L (t)) are the probabilities that the wavefunctions of the system are localized at the well on the left-hand side (σ z = −1) and that on the right-hand side (σ z = 1), respectively, at time t. The transition rate Γ is calculated via second-order perturbation with respect to the Hamiltonian H S as follows [11]: Note that this expression for the transition rate of incoherent tunneling is valid when Γ k B T [50]. By solving the master equation [equation (36)], the symmetrized correlation function is calculated as In contrast to sequential tunneling, S(ω) has only one peak at ω = 0 with a width of 2Γ, indicating the destruction of the superposition of the two localized states. The long-term dynamics are well described by the Markov approximation [11]. Therefore, one may expect that the thermal conductance in the incoherent tunneling regime would be well approximated by substituting equations (37)-(40) into equation (22). However, the results of the Markov approximation show clear deviation from the numerical results, as discussed in section 4. The reason for this is summarized as follows. Note that incoherent tunneling occurs when Γ k B T . Under this condition, the integrand of equation (22) is proportional to ω s−2 for Γ ω k B T / since S(ω) ∝ ω −2 [see equation (40)]. Then, the integral in equation (22) diverges if the high-frequency cut-off occurring due to the Bose distribution function is absent. This indicates that the high-frequency part of the integral in equation (22) makes the dominant contribution to the thermal conductance. Although the Markov approximation yields reasonable results for the low-frequency behavior of S(ω), it fails to reproduce the accurate high-frequency behavior of S(ω) in general, leading to incorrect results for the thermal conductance.

NIBA
To study the short-term (high-frequency) dynamics in the incoherent tunneling regime, we introduce the NIBA, which is a natural extension of the Markov approximation in the previous subsection [11,51]. In NIBA, the symmetrized correlation function is calculated in a manner same as that followed in a previous study [10]: where Σ(λ = −iω) is the frequency-dependent self-energy defined as Here, Q 1 (τ ) and Q 2 (τ ) are given by equations (38) and (39), respectively. The thermal conductance is then calculated by substituting equations (41) and (42) into equation (22). From the definition, it is easy to check that NIBA reproduces the Markov approximation if we neglect the frequency dependence of the self-energy and replace it with the zero-frequency value Σ(0) = 2Γ. Since NIBA appropriately considers the non-Markovian properties, it is suitable to describe the thermal conductance in the incoherent tunneling regime. The condition for NIBA is well known [11,10]. As expected from the fact that NIBA is an extension of the Markov approximation, it works well for the incoherent tunneling regime. Roughly, the incoherent tunneling mechanism becomes crucial in a regime wherein both the sequential tunneling formula and the co-tunneling formula fail. (a) NIBA holds at moderate-to-high temperatures in the sub-ohmic (s < 1) and ohmic cases (s = 1). (b) It holds for T > T * in the super-ohmic case of 1 < s < 2, where T * is the crossover temperature discussed in section 3.2. Note that NIBA never holds for s ≥ 2 since the crossover temperature T * diverges.
Here, the NIBA has been introduced to improve the Markov approximation in the incoherent regime. This introduction of the NIBA may give impression to the readers that the NIBA is a good approximation only in the incoherent regime. However, the NIBA is known to be applicable for a wider parameter region not restricted to the incoherent regime [10]. The NIBA holds also in the weak coupling regime (α 1) at arbitrary temperature for the unbiased case (ε = 0), where the interblip interaction is shown to be much weaker than the the intrablip interaction (for detailed discussion, see Sec. 21.3 in Ref. [10]). For this reason, NIBA yields almost the same result as the sequential tunneling formula or the co-tunneling formula if the system-reservoir coupling is sufficiently weak.
In section 4, we show that NIBA is an excellent approximation for reproducing the numerical results for a wide region of the parameter space at moderate-to-high temperatures. Thus, the short-term (high-frequency) non-Markovian behavior in the system dynamics is important for calculating the thermal conductance in the incoherent tunneling regime.

Numerical Results and Comparison with Analytical Formulas
While the analytical approaches discussed in the previous section are sufficiently powerful for clarifying the mechanism of heat transport in a two-state system, the detailed conditions justifying each approximation are not trivial. To understand all features of heat transport, unbiased numerical simulation without any approximation would be helpful. In this section, we therefore perform numerical simulations based on the quantum Monte Carlo method and compare the simulation results with the analytical formulas introduced in section 3. After briefly describing the numerical method in section 4.1, we separately consider the ohmic (section 4.2), sub-ohmic (section 4.3), and super-ohmic cases (sections 4.4 and 4.5).
The dynamics of the spin-boson model has been studied by using various numerical methods [52, 53, 54, 55, 56, 57, 58]. However no systematic comparisons between analytical approximations and numerical simulations has been performed in the context of heat transport near thermal equilibrium. This comparison allows us to discuss the validity of various approximations critically.

Numerical method
For numerical simulations, we employ the continuous-time quantum Monte Carlo (CTQMC) algorithm proposed in a previous study [19]. According to this algorithm, the partition function is rewritten in path-integral form with respect to an imaginary time path, σ z (τ ), and the weight of this path is defined. Then, we apply the Monte Carlo method to this representation using the cluster update algorithm [59]. The details of the CTQMC method are given in Appendix C.
Using the CTQMC method, we evaluate the imaginary time spin correlation function C(τ ) and its Fourier transform as follows: where σ z (τ ) = e τ H/ σ z e −τ H/ . The dynamical susceptibility χ(ω) is obtained from C(iω n ) via analytical continuation as follows: Analytical continuation is performed by Padé approximation [60,61] or by fitting the imaginary time spin correlation function's Fourier transform to the Lorentzian function [58]. For details, see Appendix C.

The ohmic case (s = 1)
In figure 6, we show the thermal conductances for α = 0.05, 0.1, 0.5, and 0.7 as functions of temperature. We plot the graph using the normalized temperature k B T / ∆ eff and the normalized thermal conductance κ/(k B γ∆ eff ), where ∆ eff is the effective tunneling amplitude defined in equation (23). As shown in figure 6, the numerical results fall on a universal scaling curve at each value of α regardless of the ratio ∆/ω c ( 1) obtained via this normalization. This universal behavior is characteristic of the Kondolike effect [30]. In figure 6 (c), we also show the exact solution (the Toulouse point) for α = 0.5 (indicated by the brown dot-dashed line) [58,10,30]. The agreement between the numerical results and the exact solution indicates the correctness of the CTQMC simulation.
At low temperatures (k B T ∆ eff ), the numerical results agree well with those of the approximate formula for the co-tunneling process [equation (34); indicated by blue dashed lines in figure 6]. In this regime, the thermal conductance is always proportional to T 3 (= T 2s+1 ), which is consistent with both results of a previous study [30].
At moderate (k B T ∼ ∆ eff ) and high temperatures (k B T ∆ eff ), the numerical results deviate from the co-tunneling formula and agree well with NIBA (indicated by black solid lines in figure 6). Note that the thermal conductance obtained by NIBA is proportional to T 3−2α at low temperatures, as shown in figure 6. NIBA agrees well even with the low-temperature numerical results for the weak system-reservoir coupling (α 1), whereas it deviates from these results as this coupling becomes large. It is remarkable that NIBA agrees well with the numerical results at arbitrary temperatures for α 1, as shown in figure 6 (a). In figures 6 (a) and (b), we also show the approximate formula for sequential tunneling (indicated by green dot-dashed lines). As shown in this figure, the sequential tunneling formula at moderate temperatures (k B T ∼ ∆ eff ) agrees with the numerical results of the weak system-reservoir coupling (α 1). However, note that NIBA agrees with the numerical results for a wider temperature region than the sequential tunneling formula.
The Markov approximation for incoherent tunneling, indicated by orange dotted lines in figure 6, clearly deviates from the numerical results for α = 0.05, 0.1, and 0.7, indicating the importance of the non-Markovian properties of the system. The Toulouse point α = 0.5 is an exception, as shown in figure 6 (c); NIBA coincides with the Markov approximation since at this point the self-energy in NIBA becomes independent of the frequency for the unbiased case [10]. A detailed discussion on the failure of the Markov approximation is given in section 4.3.
As described in section 3.1, quantum phase transition occurs at α c = 1 for the ohmic case. For α c ≥ 1, the effective tunneling amplitude ∆ eff becomes zero, indicating complete destruction of the superposition of the two localized states. Therefore, heat transport is induced by incoherent tunneling at arbitrary temperatures. In figure 7, we show the thermal conductance for α = 1.0, 1.5, and 2.0 as a function of temperature. As indicated by the black solid lines in the figure, the numerical results agree well with NIBA formula for arbitrary temperatures. Note that for α ≥ 1, the condition for the co-tunneling regime k B T ∆ eff is never satisfied. In figure 7, we also show the Markov approximation for incoherent tunneling (indicated by the orange dashed line). For α ≥ 1, the difference between NIBA and the Markov approximation is not considerably large.

The sub-ohmic case (s < 1)
We first discuss the thermal conductance for the sub-ohmic case wherein the systemreservoir coupling is below the critical value for the quantum phase transition. In figure 8 (a), we show the thermal conductance as a function of the temperature for s = 0.9, ∆/ω c = 0.01, and α = 0.1, for which the ground state is delocalized (α < α c (s, ∆)). At moderate and high temperatures, the numerical results agree well with the NIBA, which is shown by the black solid line. We note that the sequentialtunneling formula cannot be applied to the sub-ohmic case. At low temperatures (k B T ∆ eff ), the numerical results agree well with the co-tunneling formula, showing T 2s+1 -dependence.
We also show the results of the Markov approximation for incoherent tunneling by the orange dotted line in figure 8 (a). The Markov approximation clearly deviates from the numerical results. To understand the failure of the Markov approximation, we show  the numerical and analytical result of the symmetrized correlation function S(ω) as a function of ω/ω c for k B T = ω c /64 in figure 8 (b). While the Markov approximation for the incoherent tunneling process agrees with the numerical results at a low frequency, clear deviation is observed at higher frequencies; the numerical result indicates that the high-frequency decay of S(ω) is much faster than that of the Markov approximation, which is proportional to ω −2 (see equation (40)) We note that the numerical result of S(ω) is well reproduced by the NIBA at arbitrary frequencies. These observations indicate that the non-Markovian properties of the system dynamics are important for obtaining correct thermal conductance results for the sub-ohmic case. Next, let us study the effect of the quantum phase transition. Figure 9 shows the phase diagram determined by the CTQMC method. The detailed procedure for the determination of the critical point is given in Appendix D. The obtained critical system-reservoir coupling, α c , for the quantum phase transition is a function of both s and ∆ and is consistent with previous work based on the NRG calculation [17]. The quantum phase transition remarkably affects the temperature dependence of the thermal conductance. In figure 10, we show the thermal conductance as a function of the temperature for s = 0.6 and ∆/ω c = 0.01, for which a quantum phase transition occurs at α = α c = 0.0615. Figure 10 (a) shows the temperature dependence in the delocalized regime (α = 0.02 < α c ), for which ∆ eff remains finite. The numerical results agree well with the co-tunneling formula at low temperatures and with NIBA at moderate-to-high temperatures. This feature is the same as that shown in figure 8. Figure 10 (b) shows the temperature dependence in the localized regime (α = 0.1 > α c ), for which ∆ eff = 0. Reflecting the quantum phase transition, the numerical results agree with NIBA at arbitrary temperatures, as shown in figure 10 (b). Since the condition for the co-tunneling regime, k B T ∆ eff , is never satisfied for ∆ eff = 0, the thermal conductance does not show a universal T 2s+1 -dependence due to the co-tunneling process at low temperatures.

The super-ohmic case (1 < s < 2)
In figure 11, we show the numerical thermal conductance results obtained using CTQMC as a function of temperature for s = 1.5. Here, the horizontal and vertical axes are the normalized temperature k B T / ∆ eff and the normalized thermal conductance κ/(k B γ∆ eff (∆ eff /ω c ) 2s−2 ), respectively, where ∆ eff is the effective tunneling amplitude defined in equation (24). Note that there is no quantum transition for the super-ohmic case (s > 1); ∆ eff is finite for arbitrary system-reservoir couplings. At low temperatures (k B T ∆ eff ), the numerical results agree with the co-tunneling formula (indicated by blue dashed lines) and show T 2s+1 -dependence, regardless of the strength of the system-reservoir coupling. As shown in figure 11 (a), the numerical results for α = 0.1 agree with the sequential tunneling formula at moderate temperatures (k B T ∼ ∆ eff ) and with NIBA at high temperatures. However, from figure 11 (b), it is evident that the numerical results for α = 0.5 agree better with NIBA than with the sequential  tunneling formula at moderate-to-high temperatures (k B T ∆ eff ). This change can be explained by the crossover temperature T * , which separates the sequential (T < T * ) and incoherent (T > T * ) tunneling regimes [see equation (31)]. As the system-reservoir coupling α increases, the temperature region for which the numerical results agree with NIBA is widened since the crossover temperature T * is lowered.
The Markov approximation for incoherent tunneling is indicated by orange dotted lines in figure 11. The incoherent tunneling formula clearly deviates from numerical results, indicating the importance of the non-Markovian properties of the system dynamics. The origin of this disagreement is same as that for the sub-ohmic case (see section 4.3).

The super-ohmic case (s ≥ 2)
In figure 12, we show the numerical results of the thermal conductance obtained using the CTQMC method as a function of the temperature for s = 2.0. The normalization of the horizontal and vertical axes as well as the linetypes of the analytical formula are same as those in figure 11. At low temperatures, the numerical results agree well with the co-tunneling formula and show T 2s+1 -dependence, regardless of the strength of the system-reservoir coupling. In contrast to the case of 1 < s < 2, the numerical results agree with the sequential tunneling formula at moderate-to-high temperatures. This is reasonable since the crossover T * becomes of the order of ω c for s = 2.

Summary
We systematically considered heat transport via a local two-state system for all types of reservoirs, i.e., for the ohmic case (s = 1), super-ohmic case (s > 1), and sub-ohmic case (s < 1). We used the exact expression for the thermal conductance obtained from the Keldysh formalism and studied it using both analytic and numerical methods.
First, we considered the approximations of three transport processes: sequential tunneling, co-tunneling, and incoherent tunneling. In particular, we newly derived a universal formula for co-tunneling using the generalized Shiba relation, which predicts the T 2s+1 -dependence of the thermal conductance at low temperatures. We also pointed out that the Markov approximation yielded incorrect results for the thermal conductance in the incoherent tunneling regime since the non-Markovian properties are important. However, for the incoherent tunneling regime, NIBA yielded correct results.
Next, we used a continuous-time Monte Carlo algorithm and systematically compared the numerical results with those of the analytical approximation formulas. We found that all numerical results were well reproduced by one of three formulas, i.e., the sequential tunneling formula, co-tunneling formula, or NIBA. The formulas that yielded correct results are summarized in Table 1. We also showed that for 0 < s ≤ 1, the quantum phase transition between the delocalized and localized phases strongly affected the temperature dependence of the thermal conductance. For the delocalized phase (α < α c ), the thermal conductance is well described by the co-tunneling formula at low temperatures and by NIBA at moderate-to-high temperatures. On the contrary, for the localized phase (α > α c ), NIBA holds at arbitrary temperatures.
Our study is expected to provide a theoretical basis for describing heat transport via nano-scale objects. Herein, we focused on heat transport in a symmetric doublewell-shaped potential near the thermal equilibrium in the limit of ∆ ω c . The effect of asymmetry of system's potential, the cutoff-frequency dependence, and the far-formequilibrium effect constitute an important future problem. The temperature dependence of the thermal conductance in the critical regime near the quantum phase transition is also an intriguing subject for research and will be discussed elsewhere.
contour. By projection onto the real-time axis, the lesser component of equation (A.8) is rewritten as The heat current is then rewritten as where Σ < ν (t, t ) and Σ a ν (t, t ) are the lesser and advanced components, respectively, of the reservoir self-energy which are calculated as respectively. Here, n ν (ω) = (e ω/k B Tν − 1) −1 is the Bose distribution function of phonons (photons) for reservoir ν. The Fourier transformation of equation (A.10) gives where G r σz,σz (ω) and G < σz,σz (ω) are the Fourier transformations of the retarded and lesser components of the nonequilibrium Green function, respectively. Considering the conservation law of energy given by J L = − J R ≡ J , the heat current is rewritten as Here, we used I ν (ω) = α νĨ (ω). Rewriting G r σz,σz (ω) with χ(ω), we finally obtain equation (15).
If, for a moment, we ignore the other low-frequency oscillators, the wavefunctions of the two lowest energy eigenstates for the system-plus-reservoir are described by where |Ψ L and |Ψ R are given by respectively. Here, the prime symbol indicates that the product is in the range pω c < ω νk < ω c . Ψ ± νk is the ground-state wave function of the oscillator k in reservoir ν when the wavefunction of the local system is located at x = ±x 0 /2; it is obtained by translation of the ground-state wavefunction |Ψ 0 νk for the isolated oscillator as Adiabatic renormalization suggests that the tunneling amplitude is renormalized by the overlap between the ground-state wavefunctions of the oscillators for different localized states (σ z = ±1): If the renormalized tunneling amplitude ∆ (p) is less than pω c , the adiabatic renormalization can continue by reducing the factor p. If ∆ (p * ) = p * ω c holds at p = p * , adiabatic renormalization must be stopped there and the finite effective tunneling amplitude ∆ eff = ∆ (p * ) is obtained. On the contrary, if ∆ (p) < pω c holds for an arbitrary value of p, adiabatic renormalization can be completed even at p = 0, yielding an effective tunneling amplitude of zero (∆ eff = 0). For the ohmic case (s = 1), the effective tunneling amplitude is obtained as follows: , (0 ≤ α < 1), 0, (1 ≤ α). (B.8) In this paper, following Ref. [10], we employ a modified effective tunneling amplitude multiplied by a dimensionless function of α: Using this definition, equation (23) is derived. Based on equation (B.7), it is straightforward to show that the effective tunneling amplitude in the super-ohmic case (s > 1) assumes a finite value given by (24) and that it always vanishes for the sub-ohmic case (s < 1).

Appendix C. Continuous-time Quantum Monte Carlo Method
In early numerical studies [62,58,63], the Monte Carlo method has been applied directly to the long-range Ising model, which is mapped from the spin-boson model [10,23,64,65,11]. Subsequently, the continuous-time quantum Monte Carlo (CTQMC) algorithm [66,59] has been applied directly to the spin-boson model without mapping [19]. In this section, we describe the CTQMC algorithm employed in the present numerical simulation.
The partition function of the spin-boson model (5) is written in the path-integral form as [10,19] where σ z (τ )(= ±1) is a spin variable defined on the imaginary-time axis, Dσ z (τ ) indicates the integral for all possible paths σ z (τ ), and K(τ ) is a kernel defined as As shown in figure C1 (a), the path σ z (τ ) is assigned by an alternative configuration of kinks (jumps from σ z = −1 to σ z = +1) and anti-kinks (jumps from σ z = +1 to σ z = −1) and described by the positions τ i (i = 1, 2, · · · , 2n) of the kinks (q i = +1) and anti-kinks (q i = −1) as where n is the number of the pairs of kinks and anti-kinks. Note that the kinks and anti-kinks are alternatively located (q i+1 = −q i ). By substituting equation Here, we apply the CTQMC method to this partition function. The present CTQMC algorithm [66] employs a cluster-flip update similar to that in the Swendsen-Wang cluster algorithm [67]. The cluster-flip update is constructed as follows [19] (see figure C1). We consider the initial path σ z (τ ) of figure  as shown in (b-iii), and construct segment clusters. Here, σ z (s i ) is the value of σ z in the segment s i and the positions of the vertices (including the inserted ones) at the two edges of the segment s i are denoted by τ i−1 and τ i , respectively. Finally, we flip each segment cluster with probability 1/2, as shown in (b-iv), and remove the redundant vertices within segments, as shown in (b-v). The final path is then given by figure C1 (c). The Monte Carlo data presented in this paper typically represent averages over 10 3 -10 4 updates at low temperatures and 10 7 -10 8 updates at high temperatures.
To perform this continuation numerically, we usually employ Padé approximation [60,61]. For the weak coupling regime, Padé approximation yields poor results since the imaginary part of the pole nearest to the real frequency axis is small. In this case, we employ another approximation based on the fitting [58]. We assume that the spin correlation function as C(iω n ) ≈ aω 3 0 (ω n + λ) 2 + ω 2 0 + const, (C.10) where a, ω 0 , and λ are the fitting parameters determined using the least-squares method.
It is easy to obtain the dynamical susceptibility Im[χ(ω)] using the fitting function (C.10) with optimized parameters. Note that this fitting method works well for weak couplings since it is compatible with the dynamic susceptibility for the sequential tunneling process. For using the co-tunneling formula (33), we need to calculate the static susceptibility χ 0 . Typically, a simple estimate χ 0 2/( ∆ eff ) yields quantitatively correct results. However, for the sub-ohmic case, χ 0 has nontrivial temperature dependence, even at low temperatures. For this case, we numerically calculate χ 0 using the CTQMC method as follows: