Nonequilibrium thermal transport in the two-mode qubit-resonator system

Nonequilibrium thermal transport in circuit quantum electrodynamics emerges as one interdisciplinary field, due to the tremendous advance of quantum technology. Here, we study steady-state heat flow in a two-mode qubit-resonator model under the influence of both the qubit-resonator and resonator-resonator interactions. The heat current is suppressed and enhanced by tuning up resonator-resonator interaction strength with given weak and strong qubit-resonator couplings respectively, which is cooperative contributed by the eigen-mode of coupled resonators and qubit-photon scattering. Negative differential thermal conductance and significant thermal rectification are exhibited at weak qubit-resonator coupling, which are dominated by cycle transition processes. Moreover, the heat flow through the resonator decoupled from the qubit can be dramatically enhanced via the resonator-resonator interaction, which is attributed by the generation of eigen-mode channels of resonators.

Introduction

The microscopic interaction between bosonic fields (e.g., light and vibration) and quantum matter plays a cornerstone role in the study of quantum optics [1–4], quantum phononics [5, 6], and quantum technology [7–9]. The significant advance of superconducting and quantum-dot circuit quantum electrodynamics (cQED) fertilizes nonequilibrium heat transport in light-matter interacting systems, by integrating cQED components (e.g., resonator and qubit) with thermal and electric reservoirs [10–16]. Bounded by the second law of thermodynamics, the heat flow can be modulated via the nonequilibrium temperature bias [10], geometric heat pump [17, 18], and quantum correlations [19].

The seminal quantum Rabi model (QRM) is one generic system to describe quantum light-matter interaction, i.e., a single-mode photon field transversely coupled with a two-level qubit. QRM nowdays has been extensively realized in various quantum systems, ranging from quantum electrodynamics [2, 20], trapped ions [21, 22], and semiconductors with spin-orbital coupling [23]. Theoretically, though the eigensolution is rather difficult to obtain, D. Braak [24] and Chen et al. [25] separately proposed mapping ways a decade ago, trying to achieve the exact solution of QRM. While considering the interface between the subsystems in cQED and mesoscopic thermal reservoirs, the powerful platform to detect steady-state heat transport has been established [10, 12] to deepen the understanding of quantum thermodynamics, condensed-matter physics, and functional thermal devices. The analogous setup nowadays can also be found in optomechanics [26]. As a consequence, theoretical studies of quantum thermal transport in two-bath QRM have been preliminarily carried out [27–29].

Meanwhile, the longitudinal qubit-resonator model, as an alternative representative system, was also experimentally realized based on cQED devices [30–32], which is able to reach quite strong qubit-resonator coupling. Such longitudinal interaction leads to tremendous intriguing applications, e.g., generation of nonclassical photon states [33, 34], scalable quantum computation [35], and efficient quantum readout [36, 37]. Moreover, as a basic member in the multi-mode qubit-resonator family [30, 38, 39], the two-mode qubit-resonator model has been designed in cQED with longitudinal coupling [40–42], which is able to exhibit quantum switch effect [43]. While for nonequilibrium thermal transport in the longitudinal qubit-resonator system, an inspiring cooling by heating effect was reported, which causes the photon antibunchinig [44]. Besides, several typical nonlinear thermal effects were also proposed [45–47]. Hence, it should be intriguing to explore the transport pictures and microscopic mechanisms in the dissipative two-mode qubit-resonator model.

In this work, we apply the quantum dressed master equation (DME) combined with coherent photon states to study the steady-state heat current of the dissipative two-mode qubit-resonator model. The main points of this study are listed as follows: 1) The steady-state heat flow is suppressed under the influence of the resonator-resonator interaction with given weak qubit-resonator coupling, whereas it is enhanced by increasing the resonator-resonator interaction strength at strong qubit-resonator coupling, which is unavailable for the single-mode qubit-resonator model [29, 45]. The underlying mechanisms are illustrated based on the eigen-modes of resonators (4) and coherent state overlap coefficient (15). 2) The nonlinear thermal effects, e.g., negative differential thermal conductance and thermal rectification are unraveled in weak qubit-resonator coupling regime. They are attributed by two kinds of cycle transition processes (20a-20 b), which fertilizes microscopic transport pictures in the single-mode qubit-resonator model. 3) The resonator-resonator interaction induced indirect heat transport is also investigated. It is interesting to find that at the condition of identical resonators, even weak resonator-resonator interaction can significantly enhance heat current via the resonator isolated from the qubit, due to the eigen-modes of resonators.

The paper is organized as follows: In Section 2 we introduce the two-mode qubit-resonator model, obtain the corresponding eigenvalues and coherent photon eigenstates, and derive the DME. In Section 3 we study steady-state heat current under the influence of both qubit-resonator interaction and resonator-resonator coupling. Typical nonlinear thermal effects and resonator-resonator induced indirect heat transport is discussed. Finally, we give a summary in Section 4.

Model and method

Two-mode qubit-resonator model

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

$\begin{array}{ccc}\hfill {\widehat{H}}_s& =\hfill & \sum _{i=\mathrm{1,2}}{\omega }_i{\widehat{a}}_i^{†}{\widehat{a}}_i+t\left({\widehat{a}}_1^{†}+{\widehat{a}}_1\right)\left({\widehat{a}}_2^{†}+{\widehat{a}}_2\right)\hfill \\ \hfill & \hfill & +\frac{\epsilon }{2}{\widehat{\sigma }}_z+{\widehat{\sigma }}_z\sum _{i=\mathrm{1,2}}{\lambda }_i\left({\widehat{a}}_i^{†}+{\widehat{a}}_i\right).\hfill \end{array}(1)$

FIGURE 1

FIGURE 1. (Color online) (A) Schematic illustration of the nonequilibrium two-mode qubit-resonator model. The yellow circle marked with ${\widehat{a}}_i$ denotes the i-th optical resonator, the black wave curve describes the interaction between nearest-neighbouring resonators, the blue circle embedded with two horizontal lines and marked with $\widehat{\sigma }$ shows the two-level qubit, the double-arrowed line with λ_i means the interaction between the i-th resonator and the qubit, and the orange (blue) rectangles is the r (q)-th bosonic thermal reservoir, characterized as the temperature T_r (T_q). (B) eigen-energies E₁/ω_a and E₂/ω_a at Eq. 4 as functions of t/ω_a, with ω₁ = ω₂ = ω_a.

Where ${\widehat{a}}_i^{†}({\widehat{a}}_i)$ creates (annihilates) one photon in the i-th resonator with the frequency ω_i, t describes inter-resonator hopping strength, ${\widehat{\sigma }}_{\alpha }(\alpha =x,y,z)$ denotes the Pauli operator of the qubit, composed of two states |↑⟩ and |↓⟩, e.g., ${\widehat{\sigma }}_x=|\uparrow \rangle \langle \downarrow |+|\downarrow \rangle \langle \uparrow |$ and ${\widehat{\sigma }}_z=|\uparrow \rangle \langle \uparrow |-|\downarrow \rangle \langle \downarrow |$, ɛ 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_{\mathit{mech}}=t^{\prime }{({\widehat{x}}_1-{\widehat{x}}_2)}^2$, should generally include the reorganization terms (${\widehat{x}}_1^2$ and ${\widehat{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 ${\widehat{x}}_i=\frac{1}{\sqrt{2{\omega }_i}}({\widehat{a}}_i^{†}+{\widehat{a}}_i)$ and ${\widehat{p}}_i=i\sqrt{\frac{{\omega }_i}{2}}({\widehat{a}}_i^{†}-{\widehat{a}}_i)$, the Hamiltonian ${\widehat{H}}_s$ can be reexpressed as ${\widehat{H}}_s={\sum }_{i=\mathrm{1,2}}({\omega }_i^2{\widehat{x}}_i^2+{\widehat{p}}_i^2)/2+2t\sqrt{{\omega }_1{\omega }_2}{\widehat{x}}_1{\widehat{x}}_2-({\omega }_1+{\omega }_2)/2+\frac{\epsilon }{2}{\widehat{\sigma }}_z+{\widehat{\sigma }}_z{\sum }_{i=\mathrm{1,2}}{\lambda }_i\sqrt{2{\omega }_i}{\widehat{x}}_i$. To diagonalize ${\widehat{H}}_s$, we further introduce the canonical coordinate operators ${\widehat{q}}_1=\cos \frac{\theta }{2}{\widehat{x}}_1-\sin \frac{\theta }{2}{\widehat{x}}_2$ and ${\widehat{q}}_2=\sin \frac{\theta }{2}{\widehat{x}}_1+\cos \frac{\theta }{2}{\widehat{x}}_2$, with the angle

$\tan \theta =4t\sqrt{{\omega }_1{\omega }_2}/\left({\omega }_2^2-{\omega }_1^2\right).(2)$

Hence, the Hamiltonian (Section 1) is reexpressed as

${\widehat{H}}_s=\sum _{i=\mathrm{1,2}}{E}_i{\widehat{d}}_i^{†}{\widehat{d}}_i+\frac{\epsilon }{2}{\widehat{\sigma }}_z+{\widehat{\sigma }}_z\sum _{i=\mathrm{1,2}}{\mathrm{\Lambda }}_i\left({\widehat{d}}_i^{†}+{\widehat{d}}_i\right),(3)$

where ${\widehat{d}}_i=(\sqrt{E_i/2}{\widehat{q}}_i+i\sqrt{1/2E_i}{\widehat{p}}_{q_i})$ is the new bosonic operator, with the momentum operator ${\widehat{p}}_{q_i}=i\sqrt{E_i/2}({\widehat{d}}_i^{†}-{\widehat{d}}_i)$, which creates (annihilates) one photon with the eigen-mode energy

$E_i=\sqrt{\frac{1}{2}\left[\left({\omega }_1^2+{\omega }_2^2\right)+{\left(-1\right)}^i\sqrt{{\left({\omega }_2^2-{\omega }_1^2\right)}^2+16t^2{\omega }_1{\omega }_2}\right]},(4)$

under the bound $t\le \sqrt{{\omega }_1{\omega }_2}/2$ as shown in Figure 1B, and the modified qubit-resonator coupling strengths are given by.

${\mathrm{\Lambda }}_1=\sum _i{\lambda }_i\cos \frac{{\theta }_i}{2}\sqrt{{\omega }_i/E_1},(\mathrm{5a})$

${\mathrm{\Lambda }}_2=\sum _i{\lambda }_i\sin \frac{{\theta }_i}{2}\sqrt{{\omega }_i/E_2},(\mathrm{5b})$

With θ₁ = θ, θ₂ = θ + π. The constant coefficient in Eq. 3 is ignored for safety for dissipative dynamics and steady state. Hence, the eigenstates are expressed as

$|{\psi }_{\mathbf{n}}^{\sigma }\rangle =|\sigma \rangle \otimes {\mathrm{\Pi }}_i\left\{\exp \left[-\left({\mathrm{\Lambda }}_i^{\sigma }/E_i\right)\left({\widehat{d}}_i^{†}-{\widehat{d}}_i\right)\right]\frac{{\left({\widehat{d}}_i^{†}\right)}^{n_i}}{\sqrt{n_i!}}|0{\rangle }_{d_i}\right\},(6)$

with σ = {↑, ↓}, ${\mathrm{\Lambda }}_i^{\uparrow }={\mathrm{\Lambda }}_i$, ${\mathrm{\Lambda }}_i^{\downarrow }=-{\mathrm{\Lambda }}_i$, and the vacuum state ${\widehat{d}}_i|0{\rangle }_{d_i}=0$. And the corresponding eigenvalues are shown as

$E_{\mathbf{n},\sigma }=\sum _{i=\mathrm{1,2}}\left[n_i{E}_i-{\mathrm{\Lambda }}_i^2/E_i\right]+{\epsilon }_{\sigma }/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 ${\mathrm{\Lambda }}_i^{\sigma }={\lambda }_i^{\sigma }$. Then, the eigenstates are characterized as $|{\psi }_{\mathbf{n},\sigma }\rangle =|\sigma \rangle \otimes {\mathrm{\Pi }}_{i=\mathrm{1,2}}\left\{\exp [-({\lambda }_i^{\sigma }/{\omega }_i)({\widehat{a}}_i^{†}-{\widehat{a}}_i)]\frac{{({\widehat{a}}_i^{†})}^{n_i}}{\sqrt{n_i!}}|0{\rangle }_{a_i}\right\}$, and the corresponding eigenvalues are shown as $E_{\mathbf{n},\sigma }={\sum }_i(n_i{\omega }_i-{\lambda }_i^2/{\omega }_i)+{\epsilon }_{\sigma }/2$.

Quantum dressed master equation

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

${\widehat{H}}_{\mathit{tot}}={\widehat{H}}_s+\left({\widehat{H}}_r+\sum _{i=1}^2{\widehat{V}}_{r_i}\right)+\left({\widehat{H}}_q+{\widehat{V}}_q\right).(8)$

Specifically, the r-th thermal reservoir connecting is described as ${\widehat{H}}_r={\sum }_k{\omega }_{k,r}{\widehat{a}}_{k,r}^{†}{\widehat{a}}_{k,r}$, where ${\widehat{a}}_{k,r}^{†}({\widehat{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}=\left({\widehat{a}}_i^{†}+{\widehat{a}}_i\right)\sum _k\left(g_{k,r_i}{\widehat{a}}_{k,r}^{†}+\mathrm{H}.\mathrm{c}.\right),(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 ${\gamma }_{r_i}(\omega )=2\pi {\sum }_k|g_{k,r_i}{|}^2\delta (\omega -{\omega }_k)$. And the r-th reservoir induced interference is ignored [50], i.e., $2\pi {\sum }_k{g}_{k,r_i}{g}_{k,r_j}^{\ast }\delta (\omega -{\omega }_k)=0$ for i ≠ j. While the q-th thermal reservoir coupled with the qubit is described as ${\widehat{H}}_q={\sum }_k{\omega }_{k,q}{\widehat{a}}_{k,q}^{†}{\widehat{a}}_{k,q}$, where ${\widehat{a}}_{k,q}^{†}({\widehat{a}}_{k,q})$ means the bosonic creation (annihilation) operator with the frequency ω_k,q. The qubit-reservoir interaction is given by

${\widehat{V}}_q={\widehat{\sigma }}_x\sum _k\left(g_{k,q}{\widehat{a}}_{k,q}^{†}+\mathrm{H}.\mathrm{c}.\right),(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 ${\widehat{H}}_s$. Since we focus on the steady-state behavior of the qubit-resonator hybrid system, we properly obtain the DME as [51–53].

$\begin{array}{ccc}\hfill \frac{d}{dt}{\widehat{\rho }}_s\left(t\right)& =\hfill & -i\left[{\widehat{H}}_s,{\widehat{\rho }}_s\left(t\right)\right]+\sum _{u=r_1,r_2,q;\mathrm{n},{\mathrm{n}}^{\prime },\sigma ,{\sigma }^{\prime }}\times \hfill \\ \hfill & \hfill & \left\{{\mathrm{\Gamma }}_u^{+}\left(E_{\text{n},\sigma }^{{\text{n}}^{\prime },{\sigma }^{\prime }}\right)\mathcal{L}\left[|{\psi }_{{\text{n}}^{\prime }}^{{\sigma }^{\prime }}\rangle \langle {\psi }_{\text{n}}^{\sigma }|\right]{\widehat{\rho }}_s\left(t\right)\right.\hfill \\ \hfill & \hfill & \left.+{\mathrm{\Gamma }}_u^{-}\left(E_{\text{n},\sigma }^{{\text{n}}^{\prime },{\sigma }^{\prime }}\right)\mathcal{L}\left[|{\psi }_{\text{n}}^{\sigma }\rangle \langle {\psi }_{{\text{n}}^{\prime }}^{{\sigma }^{\prime }}|\right]{\widehat{\rho }}_s\left(t\right)\right\},\hfill \end{array}(11)$

Where ${\widehat{\rho }}_s(t)$ is the reduced density operator of the multi-mode qubit-resonator system and the dissipator is given by

$\mathcal{L}\left[\widehat{D}\right]{\widehat{\rho }}_s\left(t\right)=\widehat{D}{\widehat{\rho }}_s\left(t\right){\widehat{D}}^{†}-\frac{1}{2}\left[{\widehat{D}}^{†}\widehat{D}{\widehat{\rho }}_s\left(t\right)+{\widehat{\rho }}_s\left(t\right){\widehat{D}}^{†}\widehat{D}\right].(12)$

The transition rate contributed by the i-th resonator and r-th reservoir is described as

$\begin{array}{ccc}\hfill {\mathrm{\Gamma }}_{r_i}^{\pm }\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\sigma }\right)& =\hfill & {\delta }_{n_1,n_1^{\prime }-1}{\delta }_{n_2,n_2^{\prime }}{\gamma }_{r_i}\left(\pm E_1\right)n_r\left(\pm E_1\right)\times \hfill \\ \hfill & \hfill & n_1^{\prime }\frac{{\omega }_i}{E_1}{\cos }^2\frac{{\theta }_i}{2}\hfill \\ \hfill & \hfill & +{\delta }_{n_1,n_1^{\prime }}{\delta }_{n_2,n_2^{\prime }-1}{\gamma }_{r_i}\left(\pm E_2\right)n_r\left(\pm E_2\right)\times \hfill \\ \hfill & \hfill & n_2^{\prime }\frac{{\omega }_i}{E_2}{\sin }^2\frac{{\theta }_i}{2},\hfill \end{array}(13)$

By considering relations ${\widehat{a}}_j^{†}+{\widehat{a}}_j=\cos \frac{{\theta }_j}{2}\sqrt{\frac{{\omega }_j}{E_1}}({\widehat{d}}_1^{†}+{\widehat{d}}_1)+\sin \frac{{\theta }_j}{2}\sqrt{\frac{{\omega }_j}{E_2}}({\widehat{d}}_2^{†}+{\widehat{d}}_2)$, with γ_u( − ω)n_u( − ω) = γ_u(ω)[1 + n_u(ω)], $E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}=E_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}-E_{\mathbf{n},\sigma }$, 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

$\begin{array}{ccc}\hfill {\mathrm{\Gamma }}_q^{\pm }\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\bar{\sigma }}\right)& =\hfill & \mathrm{\Theta }\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right){\gamma }_q\left(\pm E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\bar{\sigma }}\right)n_q\left(\pm E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\bar{\sigma }}\right)\times \hfill \\ \hfill & \hfill & D_{n_1,n_1^{\prime }}^2\left(2{\mathrm{\Lambda }}_1/E_1\right)D_{n_2,n_2^{\prime }}^2\left(2{\mathrm{\Lambda }}_2/E_2\right),\hfill \end{array}(14)$

Where Θ(x) = 1 for x > 0 and Θ(x) = 0 for x ≤ 0, the coherent state overlap coefficient is specified as [54, 55].

$D_{nm}\left(x\right)=e^{-x^2/2}\sum _{l=0}^{\min \left[n,m\right]}\frac{{\left(-1\right)}^l\sqrt{n!m!}{x}^{n+m-2l}}{\left(n-l\right)!\left(m-l\right)!l!},(15)$

and the Bose-Einstein distribution function is given by $n_u(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }})=1/[\exp (E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}/k_B{T}_u)-1]$, with k_B the Boltzmann constant and T_u the temperature of the u-th thermal reservoir. Both ${\mathrm{\Gamma }}_u^{+}(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\bar{\sigma }})$ and ${\mathrm{\Gamma }}_u^{-}(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\bar{\sigma }})$ become zero as $E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\bar{\sigma }}=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., ${\widehat{d}}_1$ and ${\widehat{d}}_2$). For each channel, the local detailed balance relation is still valid, i.e., ${\gamma }_{r_i}(E_j)n_r(E_j)/\left({\gamma }_{r_i}(E_j)[1+n_r(E_j)]=\exp (-E_j/k_B{T}_r)\right.$. 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

$\begin{array}{ccc}\hfill J_u& =\hfill & \sum _{\mathbf{n},{\mathbf{n}}^{\prime },\sigma ,{\sigma }^{\prime }}\mathrm{\Theta }\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right)E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\times \hfill \\ \hfill & \hfill & \left[{\mathrm{\Gamma }}_u^{-}\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right){\rho }_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}-{\mathrm{\Gamma }}_u^{+}\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right){\rho }_{\mathbf{n},\sigma }\right],\hfill \end{array}(16)$

Where the dynamical equation of the density matrix element is given by ${\rho }_{\mathbf{n},\sigma }=\langle {\psi }_{\mathbf{n}}^{\sigma }|{\widehat{\rho }}_s|{\psi }_{\mathbf{n}}^{\sigma }\rangle$ 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.

Results and discussions

Steady-state heat current

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.

FIGURE 2

FIGURE 2. (Color online) (A) Steady state heat current J_q/ω_a as a function of the qubit-resonator interaction strength λ/ω_a, with typical resonator-resonator couplings t/ω_a. (B) Heat current J_q/ω_a as functions of λ/ω_a and t/ω_a. (C) Schematic illustration of the Hamiltonian at Eq. 3 with Λ₁ = 0. (D) and (E) are cyclic current components in Eq. 20a and Eq. 20b, respectively. Other system parameters are given by ɛ = ω_a, λ₁ = λ₂ = λ, α₁ = α₂ = 0.001ω_a, ω_c = 20ω_a, T_r = 1.5ω_a, and T_q = 0.5ω_a.

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 ${\mathrm{\Lambda }}_2=\sqrt{2}\lambda$. Hence, the bosonic eigen-mode ${\widehat{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_{nm}(2{\mathrm{\Lambda }}_2/E_2)\approx {(-1)}^n[{\delta }_{n,m}+(2{\mathrm{\Lambda }}_2/E_2)\sqrt{n+1}{\delta }_{n,m-1}-(2{\mathrm{\Lambda }}_2/E_2)\sqrt{n}{\delta }_{n,m+1}]$. Then, the transition rate assisted by the q-th reservoir is approximated as.

$\begin{array}{ll}\hfill {\mathrm{\Gamma }}_q^{\pm }\left(E_{\mathbf{n},\downarrow }^{{\mathbf{n}}^{\prime },\uparrow }\right)\approx & {\delta }_{n_1,n_1^{\prime }}\left[{\delta }_{n_2,n_2^{\prime }}{\gamma }_q\left(\pm \epsilon \right)n_q\left(\pm \epsilon \right)\right.\hfill \\ \hfill & +{\delta }_{n_2,n_2^{\prime }-1}{\gamma }_q\left(\pm \left(E_2+\epsilon \right)\right)n_q\left(\pm \left(E_2+\epsilon \right)\right)\times \hfill \\ \hfill & \left.n_2^{\prime }{\left(\frac{2{\mathrm{\Lambda }}_2}{E_2}\right)}^2\right],\hfill \end{array}(\mathrm{17a})$

$\begin{array}{ll}\hfill {\mathrm{\Gamma }}_q^{\pm }\left(E_{\mathbf{n},\uparrow }^{{\mathbf{n}}^{\prime },\downarrow }\right)\approx & {\delta }_{n_1,n_1^{\prime }}{\delta }_{n_2,n_2^{\prime }-1}{\gamma }_q\left(\pm \left(E_2-\epsilon \right)\right)n_q\left(\pm \left(E_2-\epsilon \right)\right)\times \hfill \\ \hfill & n_2^{\prime }{\left(\frac{2{\mathrm{\Lambda }}_2}{E_2}\right)}^2.\hfill \end{array}(\mathrm{17b})$

And the zeroth order of steady state population is obtained as.

$P_{\mathbf{n},\uparrow }^{\left(0\right)}\approx \frac{\left(1-e^{-{\beta }_r{E}_1}\right)\left(1-e^{-{\beta }_r{E}_2}\right)}{e^{{\beta }_q\epsilon }+1}{e}^{-{\beta }_r\left(n_1{E}_1+n_2{E}_2\right)},(\mathrm{18a})$

$P_{\mathbf{n},\downarrow }^{\left(0\right)}\approx \frac{\left(1-e^{-{\beta }_r{E}_1}\right)\left(1-e^{-{\beta }_r{E}_2}\right)}{e^{-{\beta }_q\epsilon }+1}{e}^{-{\beta }_r\left(n_1{E}_1+n_2{E}_2\right)},(\mathrm{18b})$

With β_u = 1/k_BT_u. Moreover, if we reexpress the populations in the vector form |ρ_s⟩, the steady-state solution based on Eq. 11 becomes $\mathcal{L}|{\rho }_s\rangle =0$. Then, we approximately expand $|{\rho }_s\rangle \approx |{\rho }_s^{(0)}\rangle +{({\mathrm{\Lambda }}_2/E_2)}^2|{\rho }_s^{(1)}\rangle$ and $\mathcal{L}\approx {\mathcal{L}}^{(0)}+{({\mathrm{\Lambda }}_2/E_2)}^2{\mathcal{L}}^{(1)}$. The steady-state solution is given by ${\mathcal{L}}^{(0)}|{\rho }_s^{(0)}\rangle =0$ and ${\mathcal{L}}^{(0)}|{\rho }_s^{(1)}\rangle +{\mathcal{L}}^{(1)}|{\rho }_s^{(0)}\rangle =0$. Accordingly, the heat current based on the expression (16) is can be obtained, which is contributed by two cyclic components

$J_q\approx {\left(\frac{2{\mathrm{\Lambda }}_2}{E_2}\right)}^2\times E_2\left[I_1\left(E_2\right)+I_2\left(E_2\right)\right],(19)$

where.

$\begin{array}{ll}\hfill I_1\left(E_2\right)=& \frac{{\gamma }_q\left(E_2+\epsilon \right)}{2n_q\left(\epsilon \right)+1}\left[\left(1+n_q\left(\epsilon +E_2\right)\right)n_q\left(\epsilon \right)n_r\left(E_2\right)\right.\hfill \\ \hfill & \left.-n_q\left(\epsilon +E_2\right)\left(1+n_q\left(\epsilon \right)\right)\left(1+n_r\left(E_2\right)\right)\right],\hfill \end{array}(\mathrm{20a})$

$\begin{array}{ll}\hfill I_2\left(E_2\right)=& \frac{{\gamma }_q\left(E_2-\epsilon \right)}{2n_q\left(\epsilon \right)+1}\left[\left(1+n_q\left(E_2-\epsilon \right)\right)\left(1+n_q\left(\epsilon \right)\right)n_r\left(E_2\right)\right.\hfill \\ \hfill & \left.-n_q\left(E_2-\epsilon \right)n_q\left(\epsilon \right)\left(1+n_r\left(E_2\right)\right)\right],\hfill \end{array}(\mathrm{20b})$

Which are described in Figure 2D,E. Physically, the loop of I₁(E₂) is established as $|{\psi }_{n_1,n_2+1}^{\uparrow }\rangle \to |{\psi }_{n_1,n_2}^{\downarrow }\rangle \to |{\psi }_{n_1,n_2}^{\uparrow }\rangle \to |{\psi }_{n_1,n_2+1}^{\uparrow }\rangle$, whereas the cycle transitions of I₂(E₂) is specified as $|{\psi }_{n_1,n_2+1}^{\uparrow }\rangle \to |{\psi }_{n_1,n_2+1}^{\downarrow }\rangle \to |{\psi }_{n_1,n_2}^{\uparrow }\rangle \to |{\psi }_{n_1,n_2+1}^{\uparrow }\rangle$. 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 ${\mathrm{\Lambda }}_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}=2J_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.

Negative differential thermal conductance and thermal rectification

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({\lambda }^2/{\omega }_a)I_1({\omega }_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.

FIGURE 3

FIGURE 3. (Color online) (A) Steady-state heat current into the q-th reservoir J_q/ω_a; (B) current components I₁ (20a), I₂ (20b), and I₁ + I₂ as functions of k_BΔT/ω_a at t = 0.01ω_a; (C) thermal rectification factor R as a function of the temperature bias and resonator-resonator coupling strength; (D) the asymmetric behavior of heat current J_q/ω_a by tuning k_BΔT/ω_a. Other system parameters are given by ω₁ = ω₂ = ω_a, ɛ = ω_a, λ₁ = λ₂ = 0.01ω_a, α₁ = α₂ = 0.001ω_a, ω_c = 20ω_a, T_r = ω_a/k_B + ΔT/2, and T_q = ω_a/k_B − ΔT/2.

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=\frac{|J_q\left(\mathrm{\Delta }T\right)+J_q\left(-\mathrm{\Delta }T\right)|}{|J_q\left(\mathrm{\Delta }T\right)-J_q\left(-\mathrm{\Delta }T\right)|},(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.

Resonator-resonator coupling enhances heat transport

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 ${\mathrm{\Lambda }}_1=-{\lambda }_2\sqrt{{\omega }_2/(2E_1)}$ and ${\mathrm{\Lambda }}_2={\lambda }_2\sqrt{{\omega }_2/(2E_2)}$, which implies that two eigen-modes of resonators, i.e., ${\widehat{d}}_1$ and ${\widehat{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, ${\mathrm{\Gamma }}_{r_1}^{\pm }(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\sigma })={\mathrm{\Gamma }}_{r_2}^{\pm }(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },\sigma })$ 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=2J_{r_2}$ as exhibited in Figure 5C.

FIGURE 4

FIGURE 4. (Color online) (A) Schematic illustration of the two-mode qubit-resonator model with λ₁ = 0 and (B) the corresponding eigen-mode model.

FIGURE 5

FIGURE 5. (Color online) Steady-state heat currents as functions of qubit-resonator (λ/ω_a) and resonator-resonator (t/ω_a) coupling strengths. We obtain (A–C) at ω₁ = ω₂ = ω_a and (D–F) at ω₁ = 0.9ω_a and ω₂ = 1.1ω_a. Other system parameters are given by ɛ = ω_a, λ₁ = 0, λ₂ = λ, α₁ = α₂ = 0.001ω_a, ω_c = 20ω_a, T_r = 1.5ω_a, and T_q = 0.5ω_a.

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\approx 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}$.

Conclusion

To summarize, we study nonequilibrium thermal transport in the dissipative two-mode qubit-resonator model by applying quantum DME combined with extended coherent photon states, where optical resonators and the qubit are individually coupled with thermal reservoirs. By properly treating the resonator-resonator interaction, we obtain the eigensolution of the two-mode qubit-resonator model. Then, we analyze the influence of the qubit-resonator coupling and resonator-resonator interaction on steady-state behaviors of heat currents. The currents are suppressed with the increase of resonator-resonator interaction strength at weak qubit-resonator coupling. It mainly results from the eigen-mode energies of coupled resonators (4), based on the analytical expression of cycle heat current components (19). In contrast, the currents are monotonically enhanced in strong qubit-resonator coupling regime, which is mainly attributed to the reduction of the effective qubit-resonator coupling in Eq. 15 by tuning up resonator-resonator interaction strength. Hence, the resonator hybridization and the directional cycle transition processes cooperatively contribute to the nontrivial behaviors of steady-state heat currents.

Inspired by cycle transition components of heat current at weak qubit-resonator coupling, we also investigate two representative nonlinear thermal effects, i.e., NDTC and thermal rectification. NDTC is unraveled at large temperature bias and keeps robust even with strong resonator-resonator interaction strength. The microscopic cycle transition processes (20a) and (20b) determine the appearance of NDTC. Meanwhile, the significant thermal rectification effect is observed in a wide regime of qubit-resonator and resonator-resonator couplings, which becomes perfect in absence of the resonator-resonator interaction.

Moreover, we show the effect of the resonator-resonator interaction on indirectly heat transport, by tuning off the interaction between the 1-st resonator and the qubit. It is interesting to find that for identical resonators the heat current $J_{r_1}$ (flowing into r-th reservoir mediated by the 1-th resonator) dramatically becomes identical with $J_{r_2}$ even at weak resonator-resonator coupling. While for the case of distinct resonators, the effective angle θ gradually increases from 0 to π/2 with the increase of resonator-resonator interaction strength, which leads to the monotonic enhancement of $J_{r_1}$. Therefore, the resonator-resonator interaction can be considered as one route to efficiently realize the indirect heat transport.

We hope that the analysis of quantum heat transport and thermal management in the two-mode qubit-phonon system may provide physical insight for smart energy control in photon-based hybrid quantum systems. In the future, it should be intriguing to explore the steady-state heat flows in multi-mode qubit-resonator systems, e.g., three-mode qubit-resonator model [70].

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

F-YW: analytical and numerical work, manuscript writing J-CL: numerical work ZW: analytical work CW: guiding project, analysis, manuscript writing JR: guiding project, analysis, manuscript writing.

Funding

F-YW, L-WD, and CW are supported by the National Natural Science Foundation of China under Grant No. 11704093 and the Opening Project of Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology. J-CL, ZW, and JR acknowledges the support by the National Natural Science Foundation of China (No. 11935010 and No. 11775159), Natural Science Foundation of Shanghai (No. 18ZR1442800 and No. 18JC1410900), and the China Postdoctoral Science Foundation (Grant No. 2020M681376).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Scully MO, Zubairy MS. Quantum optics. London: Cambridge University Press (1997).

Google Scholar

2. Forn-Diaz P, Lamata L, Rico E, Kono J, Solano E. Ultrastrong coupling regimes of light-matter interaction. Rev Mod Phys (2019) 91:025005. doi:10.1103/RevModPhys.91.025005

CrossRef Full Text | Google Scholar

3. Kockum AF, Miranowicz A, De Liberato S, Savasta S, Nori F. Ultrastrong coupling between light and matter. Nat Rev Phys (2019) 1:19. doi:10.1038/s42254-018-0006-2

CrossRef Full Text | Google Scholar

4. Blais A, Girvin SM, Oliver WD. Quantum information processing and quantum optics with circuit quantum electrodynamics. Nat Phys (2020) 16:247–56. doi:10.1038/s41567-020-0806-z

CrossRef Full Text | Google Scholar

5. Clerk AA. To see a SAW. Nat Phys (2012) 8:256–7. doi:10.1038/nphys2243

CrossRef Full Text | Google Scholar

6. Gustafsson MV, Santos PV, Johansson G, Delsing P. Local probing of propagating acoustic waves in a gigahertz echo chamber. Nat Phys (2012) 8:338–43. doi:10.1038/nphys2217

CrossRef Full Text | Google Scholar

7. Wallquist M, Hammerer K, Rabl P, Lukin M, Zoller P. Hybrid quantum devices and quantum engineering. Phys Scr (2009) T137:014001. doi:10.1088/0031-8949/2009/t137/014001

CrossRef Full Text | Google Scholar

8. Kurizki G, Bertet P, Kubo Y, Molmer K, Petrosyan D, Rabl P, et al. Quantum technologies with hybrid systems. Proc Natl Acad Sci U S A (2015) 112:3866–73. doi:10.1073/pnas.1419326112

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Clerk AA, Lehnert KW, Bertet P, Petta JR, Nakamura Y. Hybrid quantum systems with circuit quantum electrodynamics. Nat Phys (2020) 16:257–67. doi:10.1038/s41567-020-0797-9

CrossRef Full Text | Google Scholar

10. Ronzani A, Karimi B, Senior J, Chang YC, Peltonen JT, Chen CD, et al. Tunable photonic heat transport in a quantum heat valve. Nat Phys (2018) 14:991–5. doi:10.1038/s41567-018-0199-4

CrossRef Full Text | Google Scholar

11. Senior J, Gubaydullin A, Karimi B, Peltonen JT, Ankerhold J, Pekola JP, et al. Heat rectification via a superconducting artificial atom. Commun Phys (2020) 3:40. doi:10.1038/s42005-020-0307-5

CrossRef Full Text | Google Scholar

12. Pekola JP, Karimi B. Colloquium: quantum heat transport in condensed matter systems. Rev Mod Phys (2021) 93:041001. doi:10.1103/revmodphys.93.041001

CrossRef Full Text | Google Scholar

13. Gubaydullin A, Thomas G, Golubev DS, Lvov D, Peltonen JT, Pekola JP, et al. Photonic heat transport in three terminal superconducting circuit. Nat Commun (2022) 13:1552. doi:10.1038/s41467-022-29078-x

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Delbecq MR, Schmitt V, Parmentier FD, Roch N, Viennot JJ, Fve G, et al. Coupling a quantum dot, fermionic leads, and a microwave cavity on a chip. Phys Rev Lett (2011) 107:256804. doi:10.1103/physrevlett.107.256804

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Mi X, Cady JV, Zajac DM, Deelman PW, Petta JR. Strong coupling of a single electron in silicon to a microwave photon. Science (2016) 355:156–8. doi:10.1126/science.aal2469

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Lu JC, Wang RQ, Ren J, Kulkarni M, Jiang JH. Quantum-dot circuit-QED thermoelectric diodes and transistors. Phys Rev B (2019) 99:035129. doi:10.1103/physrevb.99.035129

CrossRef Full Text | Google Scholar

17. Ren J, Hanggi P, Li B. Berry-phase-induced heat pumping and its impact on the fluctuation theorem. Phys Rev Lett (2010) 104:170601. doi:10.1103/physrevlett.104.170601

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wang Z, Wang LQ, Chen JZ, Wang C, Ren J. Geometric heat pump: controlling thermal transport with time-dependent modulations. Front Phys (Beijing) (2022) 17:13201. doi:10.1007/s11467-021-1095-4

CrossRef Full Text | Google Scholar

19. Micadei K, Peterson JPS, Souza AM, Sarthour RS, Oliveira IS, Landi GT, et al. Reversing the direction of heat flow using quantum correlations. Nat Commun (2019) 10:2456. doi:10.1038/s41467-019-10333-7

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Haroche S, Brune M, Raimond JM. From cavity to circuit quantum electrodynamics. Nat Phys (2020) 16:243–6. doi:10.1038/s41567-020-0812-1

CrossRef Full Text | Google Scholar

21. Pedernales JS, Lizuain I, Felicetti S, Romero G, Lamata L, Solano E, et al. Quantum Rabi model with trapped ions. Sci Rep (2015) 5:15472. doi:10.1038/srep15472

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Cai ML, Liu ZD, Zhao WD, Wu YK, Mei QX, Jiang Y, et al. Observation of a quantum phase transition in the quantum rabi model with a single trapped ion. Nat Commun (2021) 12:1126. doi:10.1038/s41467-021-21425-8

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Xie QT, Cui S, Cao JP, Amico L, Fan H. Anisotropic rabi model. Phys Rev X (2014) 4:021046. doi:10.1103/physrevx.4.021046

CrossRef Full Text | Google Scholar

24. Braak D. Integrability of the rabi model. Phys Rev Lett (2011) 107:100401. doi:10.1103/physrevlett.107.100401

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Chen QH, Wang C, He S, Liu T, Wang KL. Exact solvability of the quantum rabi model using bogoliubov operators. Phys Rev A (Coll Park) (2012) 86:023822. doi:10.1103/physreva.86.023822

CrossRef Full Text | Google Scholar

26. Yang C, Wei XR, Sheng JT, Wu HB. Phonon heat transport in cavity-mediated optomechanical nanoresonators. Nat Commun (2020) 11:4656. doi:10.1038/s41467-020-18426-4

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yamamoto T, Kato T. Heat transport through a two-level system embedded between two harmonic resonators. J Phys : Condens Matter (2021) 33:395303. doi:10.1088/1361-648x/ac1281

CrossRef Full Text | Google Scholar

28. Xu M, Stockburger JT, Ankerhold J. Heat transport through a superconducting artificial atom. Phys Rev B (2021) 103:104304. doi:10.1103/physrevb.103.104304

CrossRef Full Text | Google Scholar

29. Chen ZH, Che HX, Chen ZK, Wang C, Ren J. Tuning nonequilibrium heat current and two-photon statistics via composite qubit-resonator interaction. Phys Rev Res (2022) 4:013152. doi:10.1103/physrevresearch.4.013152

CrossRef Full Text | Google Scholar

30. Niemczyk T, Deppe F, Huebl H, Menzel EP, Hocke F, Schwarz MJ, et al. Circuit quantum electrodynamics in the ultrastrong-coupling regime. Nat Phys (2010) 6:772. doi:10.1038/nphys1730

CrossRef Full Text | Google Scholar

31. Forn-Diaz P, Lisenfeld J, Marcos D, Garcia-Ripoll JJ, Solano E, Harmans CJPM, et al. Phys Rev Lett (2010) 105:237001.

PubMed Abstract |

32. Yoshihara F, Fuse T, Ashhab S, Kakuyanagi K, Saito S, Semba K, et al. Superconducting qubit–oscillator circuit beyond the ultrastrong-coupling regime. Nat Phys (2017) 13:44–7. doi:10.1038/nphys3906

CrossRef Full Text | Google Scholar

33. Zhao YJ, Liu YL, Liu YX, Nori F. Generating nonclassical photon states via longitudinal couplings between superconducting qubits and microwave fields. Phys Rev A (Coll Park) (2015) 91:053820. doi:10.1103/physreva.91.053820

CrossRef Full Text | Google Scholar

34. Wang X, Miranowicz A, Li HR, Nori F. Multiple-output microwave single-photon source using superconducting circuits with longitudinal and transverse couplings. Phys Rev A (Coll Park) (2016) 94:053858. doi:10.1103/physreva.94.053858

CrossRef Full Text | Google Scholar

35. Richer S, DiVincenzo D. Circuit design implementing longitudinal coupling: a scalable scheme for superconducting qubits. Phys Rev B (2016) 93:134501. doi:10.1103/physrevb.93.134501

CrossRef Full Text | Google Scholar

36. Didier N, Bourassa J, Blais A. Fast quantum nondemolition readout by parametric modulation of longitudinal qubit-oscillator interaction. Phys Rev Lett (2015) 115:203601. doi:10.1103/physrevlett.115.203601

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Grimsmo AL, Smith TB. Majorana qubit readout using longitudinal qubit-resonator interaction. Phys Rev B (2019) 99:235420. doi:10.1103/physrevb.99.235420

CrossRef Full Text | Google Scholar

38. Egger DJ, Wilhelm FK. Multimode circuit quantum electrodynamics with hybrid metamaterial transmission lines. Phys Rev Lett (2013) 111:163601. doi:10.1103/physrevlett.111.163601

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Bosman SJ, Gely MF, Singh V, Bruno A, Bothner D, Steele GA, et al. Multi-mode ultra-strong coupling in circuit quantum electrodynamics. Npj Quan Inf (2017) 3:46. doi:10.1038/s41534-017-0046-y

CrossRef Full Text | Google Scholar

40. Reuther GM, Zueco D, Deppe F, Hoffmann E, Menzel EP, Weißl T, et al. Two-resonator circuit quantum electrodynamics: dissipative theory. Phys Rev B (2010) 81:144510. doi:10.1103/physrevb.81.144510

CrossRef Full Text | Google Scholar

41. Baust A, Hoffmann E, Haeberlein M, Schwarz MJ, Eder P, Goetz J, et al. Tunable and switchable coupling between two superconducting resonators. Phys Rev B (2015) 91:014515. doi:10.1103/physrevb.91.014515

CrossRef Full Text | Google Scholar

42. Baust A, Hoffmann E, Haeberlein M, Schwarz MJ, Eder P, Goetz J, et al. Ultrastrong coupling in two-resonator circuit QED. Phys Rev B (2016) 93:214501. doi:10.1103/physrevb.93.214501

CrossRef Full Text | Google Scholar

43. Mariantoni M, Deppe F, Marx A, Gross R, Wilhelm FK, Solano E, et al. Two-resonator circuit quantum electrodynamics: a superconducting quantum switch. Phys Rev B (2008) 78:104508. doi:10.1103/physrevb.78.104508

CrossRef Full Text | Google Scholar

44. Naseem MT, Mustecaplioglu OE. Antibunching via cooling by heating. Phys Rev A (Coll Park) (2022) 105:012201. doi:10.1103/physreva.105.012201

CrossRef Full Text | Google Scholar

45. Wang C, Wang LQ, Ren J. Managing quantum heat transfer in a nonequilibrium qubit-phonon hybrid system with coherent phonon states. Chin Phys Lett (2021) 38:010501. doi:10.1088/0256-307x/38/1/010501

CrossRef Full Text | Google Scholar

46. Wang C, Wang LQ, Ren J. Impact of counter-rotating-wave term on quantum heat transfer and phonon statistics in nonequilibrium qubit–phonon hybrid system*. Chin Phys B (2021) 30:030506. doi:10.1088/1674-1056/abcfa8

CrossRef Full Text | Google Scholar

47. Wang C, Chen H, Liao JQ. Nonequilibrium thermal transport and photon squeezing in a quadratic qubit-resonator system. Phys Rev A (Coll Park) (2021) 104:033701. doi:10.1103/physreva.104.033701

CrossRef Full Text | Google Scholar

48. Emary C, Brandes T. Chaos and the quantum phase transition in the dicke model. Phys Rev E (2003) 67:066203. doi:10.1103/physreve.67.066203

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Mezzacapo A, Las Heras U, Pedernales JS, DiCarlo L, Solano E, Lamata L, et al. Digital quantum rabi and dicke models in superconducting circuits. Sci Rep (2014) 4:7482. doi:10.1038/srep07482

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Li SW, Cai CY, Sun CP. Steady quantum coherence in non-equilibrium environment. Ann Phys (N Y) (2015) 360:19–32. doi:10.1016/j.aop.2015.05.004

CrossRef Full Text | Google Scholar

51. Beaudoin F, Gambetta JM, Blais A. Dissipation and ultrastrong coupling in circuit QED. Phys Rev A (Coll Park) (2011) 84:043832. doi:10.1103/physreva.84.043832

CrossRef Full Text | Google Scholar

52. Settineri A, Macri V, Ridolfo A, Omar Di Stefano , Kockum AF, Nori F, et al. Dissipation and thermal noise in hybrid quantum systems in the ultrastrong-coupling regime. Phys Rev A (Coll Park) (2018) 98:053834. doi:10.1103/physreva.98.053834

CrossRef Full Text | Google Scholar

53. Le Boite A. Theoretical methods for ultrastrong light–matter interactions. Adv Quan Tech (2020) 3:1900140. doi:10.1002/qute.201900140

CrossRef Full Text | Google Scholar

54. Chen QH, Zhang YY, Liu T, Wang KL. Numerically exact solution to the finite-size Dicke model. Phys Rev A (Coll Park) (2008) 78:051801. doi:10.1103/physreva.78.051801

CrossRef Full Text | Google Scholar

55. Ren J, Zhu JX, Gubernatis JE, Wang C, Li B. Thermoelectric transport with electron-phonon coupling and electron-electron interaction in molecular junctions. Phys Rev B (2012) 85:155443. doi:10.1103/physrevb.85.155443

CrossRef Full Text | Google Scholar

56. Nicolin L, Segal D. Quantum fluctuation theorem for heat exchange in the strong coupling regime. Phys Rev B (2011) 84:161414. doi:10.1103/physrevb.84.161414

CrossRef Full Text | Google Scholar

57. Wang C, Ren J, Cao JS. Nonequilibrium energy transfer at nanoscale: a unified theory from weak to strong coupling. Sci Rep (2015) 5:11787. doi:10.1038/srep11787

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Xu DZ, Cao JS. Non-canonical distribution and non-equilibrium transport beyond weak system-bath coupling regime: a polaron transformation approach. Front Phys (Beijing) (2016) 11:110308. doi:10.1007/s11467-016-0540-2

CrossRef Full Text | Google Scholar

59. Wang C, Ren J, Cao JS. Unifying quantum heat transfer in a nonequilibrium spin-boson model with full counting statistics. Phys Rev A (Coll Park) (2017) 95:023610. doi:10.1103/physreva.95.023610

CrossRef Full Text | Google Scholar

60. Li B, Wang L, Casati G. Negative differential thermal resistance and thermal transistor. Appl Phys Lett (2006) 88:143501. doi:10.1063/1.2191730

CrossRef Full Text | Google Scholar

61. He DH, Buyukdagli S, Hu B. Origin of negative differential thermal resistance in a chain of two weakly coupled nonlinear lattices. Phys Rev B (2009) 80:104302. doi:10.1103/physrevb.80.104302

CrossRef Full Text | Google Scholar

62. He DH, Chan HK, Hu B. Phys Rev E (2010) 81:041131.

63. Li NB, Ren J, Wang L, Zhang G, Hanggi P, Li B, et al. Colloquium: phononics: manipulating heat flow with electronic analogs and beyond. Rev Mod Phys (2012) 84:1045–66. doi:10.1103/revmodphys.84.1045

CrossRef Full Text | Google Scholar

64. Chan HK, He DH, Hu B. Scaling analysis of negative differential thermal resistance. Phys Rev E (2014) 89:052126. doi:10.1103/physreve.89.052126

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Segal D, Nitzan A. Molecular heat pump. Phys Rev E (2006) 73:026109. doi:10.1103/physreve.73.026109

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Cao XF, Wang C, Zheng H, He DH. Quantum thermal transport via a canonically transformed redfield approach. Phys Rev B (2021) 103:075407. doi:10.1103/physrevb.103.075407

CrossRef Full Text | Google Scholar

67. Ren J. Predicted rectification and negative differential spin seebeck effect at magnetic interfaces. Phys Rev B (2013) 88(R):220406. doi:10.1103/physrevb.88.220406

CrossRef Full Text | Google Scholar

68. Fornieri A, Timossi G, Bosisio R, Solinas P, Giazotto F. Negative differential thermal conductance and heat amplification in superconducting hybrid devices. Phys Rev B (2016) 93:134508. doi:10.1103/physrevb.93.134508

CrossRef Full Text | Google Scholar

69. Li B, Wang L, Casati G. Thermal diode: rectification of heat flux. Phys Rev Lett (2004) 93:184301. doi:10.1103/physrevlett.93.184301

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Zhang YY, Hu ZX, Fu LB, Luo HG, Pu H, Zhang XF, et al. Quantum phases in a quantum rabi triangle. Phys Rev Lett (2021) 127:063602. doi:10.1103/physrevlett.127.063602

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix A: Derivation of Eq. 11.

Starting from Eq. 8, we assume that system-bath interactions are weak. Under the Born-Markov approximation, we individually perturb ${\widehat{V}}_{r_i}$ and ${\widehat{V}}_q$ up to the second order to obtain the quantum generalized master equation as [52].

$\begin{array}{ccc}\hfill \frac{\partial {\widehat{\rho }}_s\left(t\right)}{\partial t}& =\hfill & -i\left[{\widehat{H}}_s,{\widehat{\rho }}_s\left(t\right)\right]+\frac{1}{2}\sum _{\omega ,{\omega }^{\prime };u=r_1,r_2,q}\times \hfill \\ \hfill & \hfill & \left\{{\kappa }_u\left({\omega }^{\prime }\right)\left[{\widehat{P}}_u\left({\omega }^{\prime }\right){\widehat{\rho }}_s\left(t\right),{\widehat{P}}_u\left(\omega \right)\right]+\mathrm{H}.\mathrm{c}.\right\},\hfill \end{array}(\mathrm{A1})$

Where the rate is given by κ_u(ω) = γ_u(ω)n_u(ω), with the Bose-Einstein distribution function n_u(ω) = 1/[ exp(ω/k_BT_u) − 1], the projecting operators based on the eigenstates are given by $[{\widehat{a}}_i^{†}(-\tau )+{\widehat{a}}_i(-\tau )]={\sum }_{\omega }{\widehat{P}}_{r_i}(\omega )e^{-i\omega \tau }$ and ${\widehat{\sigma }}_x(-\tau )={\sum }_{\omega }{\widehat{P}}_q(\omega )e^{-i\omega \tau }$, with ${\widehat{P}}_{r_i}(\omega )={\sum }_{n,m}\langle {\psi }_n|({\widehat{a}}_i^{†}+{\widehat{a}}_i)|{\psi }_m\rangle \delta (\omega -E_{nm})|{\psi }_n\rangle \langle {\psi }_m|$ and ${\widehat{P}}_q(\omega )={\sum }_{n,m}\langle {\psi }_n|\left({\widehat{\sigma }}_x|{\psi }_m\rangle \delta (\omega -E_{nm})|{\psi }_n\rangle \langle {\psi }_m|\right.$, the eigensolution is given by ${\widehat{H}}_s|{\psi }_n\rangle =E_n|{\psi }_n\rangle$, and the energy gap becomes E_nm = E_n − E_m.

It should be noted that for the finite-time evolution, the off-diagonal elements of the density operator are generally coupled with the diagonal ones in the eigenspace. However, after a long-time evolution, the off-diagonal elements gradually become ignored in the present model. Hence, the populations are naturally decoupled from the off-diagonal elements, which reduces the pair of projecting operators ${\widehat{P}}_u(\omega )$ and ${\widehat{P}}_u({\omega }^{\prime })$ to ${\widehat{P}}_u(\omega =E_{nm})=\langle {\psi }_n|{\widehat{A}}_u|{\psi }_m\rangle |{\psi }_n\rangle \langle {\psi }_m|$ and ${\widehat{P}}_u(\omega =E_{mn})=\langle {\psi }_m|{\widehat{A}}_u|{\psi }_n\rangle |{\psi }_m\rangle \langle {\psi }_n|$, with ${\widehat{A}}_{r_i}=({\widehat{a}}_i^{†}+{\widehat{a}}_i)$ and ${\widehat{A}}_q={\widehat{\sigma }}_x$. Accordingly, the generalized mater equation is reduced to the dressed master equation [51]

$\begin{array}{ccc}\hfill \frac{\partial {\widehat{\rho }}_s\left(t\right)}{\partial t}& =\hfill & -i\left[{\widehat{H}}_s,{\widehat{\rho }}_s\left(t\right)\right]+\sum _{j,k>j;u=r_1,r_2,q}\times \hfill \\ \hfill & \hfill & \left\{{\mathrm{\Gamma }}_u^{jk}\left(1+n_u\left({\mathrm{\Delta }}_{kj}\right)\right)\mathcal{L}\left[|{\varphi }_j\rangle \langle {\varphi }_k|,{\widehat{\rho }}_s\left(t\right)\right]\right.\hfill \\ \hfill & \hfill & \left.+{\mathrm{\Gamma }}_u^{jk}{n}_u\left({\mathrm{\Delta }}_{kj}\right)\mathcal{L}\left[|{\varphi }_k\rangle \langle {\varphi }_j|,{\widehat{\rho }}_s\left(t\right)\right]\right\},\hfill \end{array}(\mathrm{A2})$

Where the dissipative is given by

$\begin{array}{ccc}\hfill \mathcal{L}\left[|{\psi }_{n^{\prime }}\rangle \langle {\psi }_n|\right]{\widehat{\rho }}_s& =\hfill & |{\psi }_{n^{\prime }}\rangle \langle {\psi }_n|{\widehat{\rho }}_s|{\psi }_n\rangle \langle {\psi }_{n^{\prime }}|\hfill \\ \hfill & \hfill & -\frac{1}{2}\left(|{\psi }_n\rangle \langle {\psi }_n|{\widehat{\rho }}_s+\mathrm{H}.\mathrm{c}.\right).\hfill \end{array}(\mathrm{A3})$

Based on DME (section A2) and specifying |ψ_k⟩ = |ψ_n,σ⟩ and E_k = E_n,σ, the dynamical equation of the density matrix element ${\rho }_{\mathbf{n},\sigma }=\langle {\psi }_{\mathbf{n}}^{\sigma }|{\widehat{\rho }}_s|{\psi }_{\mathbf{n}}^{\sigma }\rangle$ is obtained as

$\begin{array}{ccc}\hfill \frac{d}{dt}{\rho }_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}& =\hfill & \sum _{\mathbf{n},\sigma ,u}\left\{\left[{\mathrm{\Gamma }}_u^{+}\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right){\rho }_{\mathbf{n},\sigma }-{\mathrm{\Gamma }}_u^{-}\left(E_{\mathbf{n},\sigma }^{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right){\rho }_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right]\right.\hfill \\ \hfill & \hfill & +\sum _{\mathbf{n},\sigma ,u}\left\{\left[{\mathrm{\Gamma }}_u^{-}\left(E_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}^{\mathbf{n},\sigma }\right){\rho }_{\mathbf{n},\sigma }-{\mathrm{\Gamma }}_u^{+}\left(E_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}^{\mathbf{n},\sigma }\right){\rho }_{{\mathbf{n}}^{\prime },{\sigma }^{\prime }}\right]\right\},\hfill \end{array}(\mathrm{A4})$

Where the rates are specified in Eqs 13, 14.