The Hamiltonian of the longitudinal two-mode qubit-resonator system, composed of one two-level qubit longitudinally coupled with two optical resonators in Figure 1A, reads H ̂ s = ∑ i = 1,2 ω i a ̂ i † a ̂ i + t a ̂ 1 † + a ̂ 1 a ̂ 2 † + a ̂ 2 + ε 2 σ ̂ z + σ ̂ z ∑ i = 1,2 λ i a ̂ i † + a ̂ i . (1)

Where a ̂ i † ( a ̂ i ) creates (annihilates) one photon in the i-th resonator with the frequency ω _i , t describes inter-resonator hopping strength, σ ̂ α ( α = x , y , z ) denotes the Pauli operator of the qubit, composed of two states |↑⟩ and |↓⟩, e.g., σ ̂ x = | ↑ 〉 〈 ↓ | + | ↓ 〉 〈 ↑ | and σ ̂ z = | ↑ 〉 〈 ↑ | − | ↓ 〉 〈 ↓ | , ɛ is the splitting energy of the qubit, and λ _i is the longitudinally coupling strength between the i-th resonator and the qubit. In this paper, we set ℏ = 1 and select the identical frequency case of resonators ω _i = ω _a by default. The generalization to the distinct case is straightforward. The resonator-resonator linear interaction could be approximately formed in several ways, e.g., collective spin-photon interaction in quantum Dicke model at normal phase under the Holstein-Primakoff transformation [48, 49] and the coupling between two coplanar waveguides in cQED [40, 43]. While for the mechanical oscillators, the resonator-resonator interaction, i.e., V mech = t ′ ( x ̂ 1 − x ̂ 2 ) 2 , should generally include the reorganization terms ( x ̂ 1 2 and x ̂ 2 2 ), which may become negligible compared to ω _a in the weak interaction limit.

Next, we try to obtain the eigenvalues and eigenstates of the Hamiltonian (Section 1). Specifically, if we include the coordinate and momentum operators of the resonators as x ̂ i = 1 2 ω i ( a ̂ i † + a ̂ i ) and p ̂ i = i ω i 2 ( a ̂ i † − a ̂ i ) , the Hamiltonian H ̂ s can be reexpressed as H ̂ s = ∑ i = 1,2 ( ω i 2 x ̂ i 2 + p ̂ i 2 ) / 2 + 2 t ω 1 ω 2 x ̂ 1 x ̂ 2 − ( ω 1 + ω 2 ) / 2 + ε 2 σ ̂ z + σ ̂ z ∑ i = 1,2 λ i 2 ω i x ̂ i . To diagonalize H ̂ s , we further introduce the canonical coordinate operators q ̂ 1 = cos θ 2 x ̂ 1 − sin θ 2 x ̂ 2 and q ̂ 2 = sin θ 2 x ̂ 1 + cos θ 2 x ̂ 2 , with the angle tan ⁡ θ = 4 t ω 1 ω 2 / ω 2 2 − ω 1 2 . (2) Hence, the Hamiltonian (Section 1) is reexpressed as H ̂ s = ∑ i = 1,2 E i d ̂ i † d ̂ i + ε 2 σ ̂ z + σ ̂ z ∑ i = 1,2 Λ i d ̂ i † + d ̂ i , (3) where d ̂ i = ( E i / 2 q ̂ i + i 1 / 2 E i p ̂ q i ) is the new bosonic operator, with the momentum operator p ̂ q i = i E i / 2 ( d ̂ i † − d ̂ i ) , which creates (annihilates) one photon with the eigen-mode energy E i = 1 2 ω 1 2 + ω 2 2 + − 1 i ω 2 2 − ω 1 2 2 + 16 t 2 ω 1 ω 2 , (4) under the bound t ≤ ω 1 ω 2 / 2 as shown in Figure 1B, and the modified qubit-resonator coupling strengths are given by. Λ 1 = ∑ i λ i ⁡ cos θ i 2 ω i / E 1 , (5a) Λ 2 = ∑ i λ i ⁡ sin θ i 2 ω i / E 2 , (5b) With θ ₁ = θ, θ ₂ = θ + π. The constant coefficient in Eq. 3 is ignored for safety for dissipative dynamics and steady state. Hence, the eigenstates are expressed as | ψ n σ 〉 = | σ 〉 ⊗ Π i exp − Λ i σ / E i d ̂ i † − d ̂ i d ̂ i † n i n i ! | 0 〉 d i , (6) with σ = {↑, ↓}, Λ i ↑ = Λ i , Λ i ↓ = − Λ i , and the vacuum state d ̂ i | 0 〉 d i = 0 . And the corresponding eigenvalues are shown as E n , σ = ∑ i = 1,2 n i E i − Λ i 2 / E i + ε σ / 2 . (7) with ɛ _↑ = ɛ, ɛ _↓ = −ɛ, and n = [n ₁, n ₂] (n _i = 0, 1, 2, ⋯ ). In absence of t, the eigen-energy E _i (4) is reduced to ω _i , and the modified coupling strengths become Λ i σ = λ i σ . Then, the eigenstates are characterized as | ψ n , σ 〉 = | σ 〉 ⊗ Π i = 1,2 exp [ − ( λ i σ / ω i ) ( a ̂ i † − a ̂ i ) ] ( a ̂ i † ) n i n i ! | 0 〉 a i , and the corresponding eigenvalues are shown as E n , σ = ∑ i ( n i ω i − λ i 2 / ω i ) + ε σ / 2 .

In practice, the quantum system inevitably interacts with surrounding environments. Hence, the qubit and optical resonators may individually interact with bosonic thermal baths. Accordingly, the total Hamiltonian, including the quantum system, thermal reservoirs, and their interactions, is expressed as H ̂ tot = H ̂ s + H ̂ r + ∑ i = 1 2 V ̂ r i + H ̂ q + V ̂ q . (8) Specifically, the r-th thermal reservoir connecting is described as H ̂ r = ∑ k ω k , r a ̂ k , r † a ̂ k , r , where a ̂ k , r † ( a ̂ k , r ) creates (annihilates) one boson in the r-th thermal reservoir with the frequency ω _k,r . The interaction between the i-th resonator and the reservoir is given by V r i = a ̂ i † + a ̂ i ∑ k g k , r i a ̂ k , r † + H . c . , (9) where g k , r i denotes the coupling strength. Here, we consider each optical resonator is separately coupled with the r-th thermal reservoir. Hence, the interaction between the i-th resonator and the reservoir can be characterized as the spectral function γ r i ( ω ) = 2 π ∑ k | g k , r i | 2 δ ( ω − ω k ) . And the r-th reservoir induced interference is ignored [50], i.e., 2 π ∑ k g k , r i g k , r j ∗ δ ( ω − ω k ) = 0 for i ≠ j. While the q-th thermal reservoir coupled with the qubit is described as H ̂ q = ∑ k ω k , q a ̂ k , q † a ̂ k , q , where a ̂ k , q † ( a ̂ k , q ) means the bosonic creation (annihilation) operator with the frequency ω _k,q . The qubit-reservoir interaction is given by V ̂ q = σ ̂ x ∑ k g k , q a ̂ k , q † + H . c . , (10) with g _k,q being the coupling strength, which is characterized as the spectral functions γ _q (ω) = 2π|g _k,q |² δ(ω − ω _k ). We select the Ohmic case of the spectral functions in this work, i.e., γ _u (ω) = α _u ω exp( − |ω|/ω _c ), with α _u the dissipation strength and ω _c the cutoff frequency.

Then, we assume the interactions between the quantum system and thermal reservoirs are weak. Under the Born-Makov approximation, we may properly perturb them to obtain the quantum master equation within the eigenstate basis of H ̂ s . Since we focus on the steady-state behavior of the qubit-resonator hybrid system, we properly obtain the DME as [51–53]. d d t ρ ̂ s t = − i H ̂ s , ρ ̂ s t + ∑ u = r 1 , r 2 , q ; n , n ′ , σ , σ ′ × Γ u + E n , σ n ′ , σ ′ L | ψ n ′ σ ′ 〉 〈 ψ n σ | ρ ̂ s t + Γ u − E n , σ n ′ , σ ′ L | ψ n σ 〉 〈 ψ n ′ σ ′ | ρ ̂ s t , (11)

Where ρ ̂ s ( t ) is the reduced density operator of the multi-mode qubit-resonator system and the dissipator is given by L D ̂ ρ ̂ s t = D ̂ ρ ̂ s t D ̂ † − 1 2 D ̂ † D ̂ ρ ̂ s t + ρ ̂ s t D ̂ † D ̂ . (12) The transition rate contributed by the i-th resonator and r-th reservoir is described as Γ r i ± E n , σ n ′ , σ = δ n 1 , n 1 ′ − 1 δ n 2 , n 2 ′ γ r i ± E 1 n r ± E 1 × n 1 ′ ω i E 1 cos 2 θ i 2 + δ n 1 , n 1 ′ δ n 2 , n 2 ′ − 1 γ r i ± E 2 n r ± E 2 × n 2 ′ ω i E 2 sin 2 θ i 2 , (13)

By considering relations a ̂ j † + a ̂ j = cos θ j 2 ω j E 1 ( d ̂ 1 † + d ̂ 1 ) + sin θ j 2 ω j E 2 ( d ̂ 2 † + d ̂ 2 ) , with γ _u ( − ω)n _u ( − ω) = γ _u (ω)[1 + n _u (ω)], E n , σ n ′ , σ ′ = E n ′ , σ ′ − E n , σ , and δ _n,n′ = 1 for n = n′ and δ _n,n′ = 0 for n ≠ n′. The nonzero rate assisted by the q-th reservoir is given by Γ q ± E n , σ n ′ , σ ̄ = Θ E n , σ n ′ , σ ′ γ q ± E n , σ n ′ , σ ̄ n q ± E n , σ n ′ , σ ̄ × D n 1 , n 1 ′ 2 2 Λ 1 / E 1 D n 2 , n 2 ′ 2 2 Λ 2 / E 2 , (14)

Where Θ(x) = 1 for x > 0 and Θ(x) = 0 for x ≤ 0, the coherent state overlap coefficient is specified as [54, 55]. D n m x = e − x 2 / 2 ∑ l = 0 min n , m − 1 l n ! m ! x n + m − 2 l n − l ! m − l ! l ! , (15) and the Bose-Einstein distribution function is given by n u ( E n , σ n ′ , σ ′ ) = 1 / [ exp ( E n , σ n ′ , σ ′ / k B T u ) − 1 ] , with k _B the Boltzmann constant and T _u the temperature of the u-th thermal reservoir. Both Γ u + ( E n , σ n ′ , σ ̄ ) and Γ u − ( E n , σ n ′ , σ ̄ ) become zero as E n , σ n ′ , σ ̄ = 0 , for the phonon with ω = 0 can not support the heat transport.

The transition rate 13) involved with the r-th thermal reservoir and the i-th resonator is individual contributed by two eigen-mode channels, (i.e., d ̂ 1 and d ̂ 2 ). For each channel, the local detailed balance relation is still valid, i.e., γ r i ( E j ) n r ( E j ) / γ r i ( E j ) [ 1 + n r ( E j ) ] = exp ( − E j / k B T r ) . While the rate 14) into the q-th reservoir via the qubit is cooperatively determined by qubit-photon scattering processes, which are characterized as the coherent state overlap coefficient (15).

Based on DME (11), the steady-state heat current mediated by the subsystem component (resonators and qubit) into the corresponding thermal reservoir is expressed as J u = ∑ n , n ′ , σ , σ ′ Θ E n , σ n ′ , σ ′ E n , σ n ′ , σ ′ × Γ u − E n , σ n ′ , σ ′ ρ n ′ , σ ′ − Γ u + E n , σ n ′ , σ ′ ρ n , σ , (16)

Where the dynamical equation of the density matrix element is given by ρ n , σ = 〈 ψ n σ | ρ ̂ s | ψ n σ 〉 and the steady state is obtained by dρ _n,σ /dt = 0. During the numerical calculation, we need to truncate the photon number, i.e., n _i ≤ 40, to practically obtain the steady state and the heat flow.

We first study the influence of resonator-resonator coupling t on the behavior of steady-state heat current into the q-th thermal reservoir in Figure 2A. Before the analysis, we approximately characterize the weak couplings as t/ω _a ≲0.01 and λ/ω _a ≲0.01, which may fulfill the relations t/ω _a ≪1 and λ/ω _a ≪1. While for strong couplings, we quantify as t/ω _a ≳0.2 and λ/ω _a ≳1, where multi-photon processes may be included. In the weak qubit-resonator coupling regime with a given λ/ω _a , J _q is monotonically suppressed by tuning the resonator-resonator interaction strength. On the contrary, at strong qubit-resonator coupling the heat current is dramatically enhanced by increasing t. Furthermore, J _q generally exhibits nonmonotonic behavior with the increase of the qubit-resonator interaction strength. We also give a comprehensive picture of the effects of both qubit-resonator and resonator-resonator couplings on the heat current in Figure 2B. It is shown that at strong λ and large t the optimal regime of J _q is broadened compared to the one with t = 0. Hence, we conclude that the resonator-resonator interaction indeed significantly contributes to the steady-state heat current.

Next, we try to analyze the heat current from the analytical view. Based on the case of identical resonators we set in Figure 2, the modified qubit-resonator coupling strengths (5a) and (5 b) are specified as Λ₁ = 0 and Λ 2 = 2 λ . Hence, the bosonic eigen-mode d ̂ 1 is effectively decoupled from the qubit, as shown in Figure 2C, which shows no contribution to the steady-state heat current. The heat transport system can be reduced to the single-mode qubit-resonator model [45].

At weak qubit-resonator coupling (λ/ω _a ≪1), the coherent state overlap coefficient 15) with Λ₂/E ₂≪1 is simplified as D n m ( 2 Λ 2 / E 2 ) ≈ ( − 1 ) n [ δ n , m + ( 2 Λ 2 / E 2 ) n + 1 δ n , m − 1 − ( 2 Λ 2 / E 2 ) n δ n , m + 1 ] . Then, the transition rate assisted by the q-th reservoir is approximated as. Γ q ± E n , ↓ n ′ , ↑ ≈ δ n 1 , n 1 ′ δ n 2 , n 2 ′ γ q ± ε n q ± ε + δ n 2 , n 2 ′ − 1 γ q ± E 2 + ε n q ± E 2 + ε × n 2 ′ 2 Λ 2 E 2 2 , (17a) Γ q ± E n , ↑ n ′ , ↓ ≈ δ n 1 , n 1 ′ δ n 2 , n 2 ′ − 1 γ q ± E 2 − ε n q ± E 2 − ε × n 2 ′ 2 Λ 2 E 2 2 . (17b) And the zeroth order of steady state population is obtained as. P n , ↑ 0 ≈ 1 − e − β r E 1 1 − e − β r E 2 e β q ε + 1 e − β r n 1 E 1 + n 2 E 2 , (18a) P n , ↓ 0 ≈ 1 − e − β r E 1 1 − e − β r E 2 e − β q ε + 1 e − β r n 1 E 1 + n 2 E 2 , (18b) With β _u = 1/k _B T _u . Moreover, if we reexpress the populations in the vector form |ρ _s ⟩, the steady-state solution based on Eq. 11 becomes L | ρ s 〉 = 0 . Then, we approximately expand | ρ s 〉 ≈ | ρ s ( 0 ) 〉 + ( Λ 2 / E 2 ) 2 | ρ s ( 1 ) 〉 and L ≈ L ( 0 ) + ( Λ 2 / E 2 ) 2 L ( 1 ) . The steady-state solution is given by L ( 0 ) | ρ s ( 0 ) 〉 = 0 and L ( 0 ) | ρ s ( 1 ) 〉 + L ( 1 ) | ρ s ( 0 ) 〉 = 0 . Accordingly, the heat current based on the expression (16) is can be obtained, which is contributed by two cyclic components J q ≈ 2 Λ 2 E 2 2 × E 2 I 1 E 2 + I 2 E 2 , (19) where. I 1 E 2 = γ q E 2 + ε 2 n q ε + 1 1 + n q ε + E 2 n q ε n r E 2 − n q ε + E 2 1 + n q ε 1 + n r E 2 , (20a) I 2 E 2 = γ q E 2 − ε 2 n q ε + 1 1 + n q E 2 − ε 1 + n q ε n r E 2 − n q E 2 − ε n q ε 1 + n r E 2 , (20b) Which are described in Figure 2D,E. Physically, the loop of I ₁(E ₂) is established as | ψ n 1 , n 2 + 1 ↑ 〉 → | ψ n 1 , n 2 ↓ 〉 → | ψ n 1 , n 2 ↑ 〉 → | ψ n 1 , n 2 + 1 ↑ 〉 , whereas the cycle transitions of I ₂(E ₂) is specified as | ψ n 1 , n 2 + 1 ↑ 〉 → | ψ n 1 , n 2 + 1 ↓ 〉 → | ψ n 1 , n 2 ↑ 〉 → | ψ n 1 , n 2 + 1 ↑ 〉 . Trough these two heat transport processes, the energy quanta E ₂ is delivered from the r-th reservoir to the q-th one. As t is tuned up, both the factor Λ 2 2 / E 2 and the current components (I ₁ + I ₂) are gradually suppressed, which leads to the monontonic suppression of J _q . Moreover, by comparing the current in the single-resonator limit N = 1, i.e., J _q (N = 1) = (4λ ²/ω _a )I ₁(ω _a ), the current without resonator-resonator coupling becomes J q t = 0 = 2 J q ( N = 1 ) . We note that though transport pictures of both the dissipative two-mode qubit-resonator model and the single-mode case are characterized as cycle transition processes, they are microscopically distinct. Specifically, 1) the modified qubit-resonator interaction (Λ _i ) and eigen-mode energy (E _i ) induced by the resonator-resonator coupling dramatically affect the expression of heat current J _q ; (ii) the existence of the transition rate results in the generation of the cycle flux component I ₂.

While at strong qubit-resonator coupling, e.g., λ/ω _a ≳1, qubit-photon scattering effects dominate the transition processes, which are characterized as the coefficient (15). Such analogous effects generally suppress the heat transport, as shown in Figure 2A and previous works, e.g., nonequilibrium spin-boson model [56–59] and quantum Rabi model [29]. In particular, the crucial factor to quantify such scattering is Λ₂/E ₂. Thus, by increasing t it is found that the factor shows monotonic reduction. As a consequence, the scattering processes are weakened accordingly, which finally enhances the heat current. Therefore, we explain the underlying mechanisms of the influence of the resonator-resonator interaction on the heat current in weak and strong qubit-resonator coupling limits.

Inspired by the cycle transition features of the steady state heat current with weak qubit-resonator interaction, we first study the negative differential thermal conductance (NDTC) in the two-mode qubit-resonator model. NDTC is widely considered as one generic nonlinear feature in nonequilibrium heat transport [60–64], which is characterized as the heat current, e.g., J _q , being reduced with the increase of the reservoir temperature bias ΔT = T _r − T _q . In phononics, NDTC was originally introduced by B. Li and his colleges to analyze heat flow in classical phononic lattices [63]. As a consequence, tremendous open quantum systems have been unraveled to exhibit NDTC. In particular, D. Segal [65] and Cao et al. [66] individually applied the noninteracting blip approximation and canonically transformed Redfield equation to explore the effect of NDTC with strong system-bath interactions. Ren et al. [67] included the spin-boson-fermion systems to propose spin-NDTC based on the Landauer formula expression of the heat current. Moreover, Giazotto et al. [68] designed the superconducting devices to measure NDTC by modulating the superconducting gap with the temperature.

Here, we show NDTC in the two-mode qubit-resonator model with identical resonators in Figure 3A. In absence of the resonator-resonator interaction, the heat current is reduced to J q t = 0 = 2 ( λ 2 / ω a ) I 1 ( ω a ) , which is fully determined via single-type cycle transition process in Figure 2D. And the microscopic picture is identical with single qubit-resonator model [45]. Specifically, it exhibits nonmonotonic behavior by tuning up the temperature bias, i.e., initial enhancement and later suppression. Finally the heat current is completely eliminated at large ΔT (e.g., T _r = 2ω _a and T _q = 0), which originates from the empty occupation of bosons in the q-th reservoir, i.e., n _q (ɛ) = 0 and n _q (ɛ + ω _a ) = 0. Then, we tune on the resonator-resonator interaction. It is found that the signature of NDTC persists in the whole temperature bias regime. However, NDTC is determined by two kinds of cycle transition processes as shown in Figure 2D,E, corresponding to I ₁ and I ₂. In particular, I ₂ dramatically contributes to NDTC even at rather weak resonator-resonator coupling, e.g., the blue-solid-triangle line at t = 0.01ω _a in Figure 3B. Moreover, the signature of NDTC is gradually suppressed with the increase of t, which stems from robustness of cycle transition process of I ₂, partially characterized as the resultant part at large bias limit, i.e., I ₂ ≈ γ _q (E ₂ − ɛ)n _r (E ₂) from Eq. 19.

We also study the thermal rectification effect by modulating the temperature bias, which is another representative nonlinear character in nonequilibrium thermal transport. The concept of thermal rectification is introduced as that the heat flow in one direction is larger than the counterpart in the opposite direction [63, 69]. We characterize the rectification by the factor R = | J q Δ T + J q − Δ T | | J q Δ T − J q − Δ T | , (21) where J _q (ΔT) denotes the current under the setting of temperatures T _r = ω _a /k _B + ΔT/2 and T _q = ω _a /k _B − ΔT/2. Such effect becomes perfect as |J _q (ΔT)|≫|J _q ( − ΔT)| (i.e., R = 1) and vanishing as J _q (ΔT) ≈ − J _q ( − ΔT) (i.e., R = 0). Figure 3C clearly exhibits the giant heat amplification factor at large temperature bias, mainly due to the asymmetric behavior of heat current J _q by tuning ΔT from negative to positive regimes, as shown in Figure 3D. Specifically, the disappearance of the current at ΔT = 2ω _a and t = 0 results in the perfect thermal rectifier, whereas the existence of resultant current at finite t comparatively reduces thermal rectification factor. Therefore, we conclude that microscopic cycle transition processes are crucial to exhibit NDTC and significant thermal rectification in two-mode qubit-resonator model.

To show the ability of the resonator-resonator interaction to modulate the heat transport, we tune off the interaction between the 1-st resonator and the qubit, as shown in Figure 4A. Straightforwardly, the heat current J r 1 always keeps vanishing as t = 0. While by tuning on resonator-resonator coupling in case of identical resonators (ω _i = ω _a , i = 1, 2), the modified qubit-photon interaction strengths become Λ 1 = − λ 2 ω 2 / ( 2 E 1 ) and Λ 2 = λ 2 ω 2 / ( 2 E 2 ) , which implies that two eigen-modes of resonators, i.e., d ̂ 1 and d ̂ 2 , will both contribute to the heat transport, as depicted in Figure 4B. It is interesting to find that even with weak resonator-resonator coupling (e.g., t/ω _a = 0.01), the heat current J r 1 in Figure 5A behaves quite similar with J r 2 in Figure 5B. Based on the expression of transition rates at Eqs. 13 with θ = π/2, Γ r 1 ± ( E n , σ n ′ , σ ) = Γ r 2 ± ( E n , σ n ′ , σ ) can be directly obtained. Each eigen-mode channel shows identical contribution to these two rates. This clearly demonstrates that actually J r 1 = J r 2 by combing the expression of the heat current at Eq. 16. It also leads to J q = 2 J r 2 as exhibited in Figure 5C.

Moreover, we analyze the effect of the resonator-resonator interaction on J r 1 with distinct resonators. With weak resonator-resonator coupling [t/(ω ₂ − ω ₁)≪1], θ becomes vanishing. Hence, the heat flow via the 1-th resonator is dramatically blocked. It is expected to see J r 1 is much weaker than J r 2 , i.e., J q ≈ J r 2 , which are also shown in Figures 5D–F. As a result, the two-mode qubit-resonator transport system is reduced to the standard single-mode qubit-resonator case [45]. While by increasing the resonator-resonator coupling, e.g., t = 0.4ω _a , it is interesting to find that J r 1 becomes comparable with J r 2 , due to dramatic hybridization of two resonators (θ = π/2). Therefore, we conclude that though the 1-th resonator is decoupled from the qubit, the inclusion of the resonator-resonator interaction will open heat-exchange channels, which significantly contributes to the heat current J r 1 .