# Tuning nonequilibrium heat current and two-photon statistics via composite
qubit-resonator interaction

Zhe-Huan Chen1 Han-Xin Che1 Zhe-Kai Chen1 Chen Wang1,
<EMAIL_ADDRESS><EMAIL_ADDRESS>Jie Ren2<EMAIL_ADDRESS>1Department of Physics, Zhejiang Normal University, Jinhua 321004, Zhejiang,
P. R. China
2Center for Phononics and Thermal Energy Science, China-EU Joint Center for
Nanophononics,
Shanghai Key Laboratory of Special Artificial Microstructure Materials and
Technology,
School of Physics Sciences and Engineering, Tongji University, Shanghai
200092, China

###### Abstract

Quantum thermal transport and two-photon statistics serve as two
representative nonequilibrium features in circuit quantum electrodynamics
systems. Here, we investigate quantum heat flow and two-photon correlation
function at steady-state in a composite qubit-resonator model, where one qubit
shows both transverse and longitudinal couplings to a single-mode optical
resonator. With weak qubit-resonator interaction, we unravel two microscopic
transport pictures, i.e., cotunneling and cyclic heat exchange processes,
corresponding to transverse and longitudinal couplings respectively. At strong
qubit-resonator coupling, the heat current exhibits nonmonotonic behavior by
increasing qubit-resonator coupling strength, which tightly relies on the
scattering processes between the qubit and corresponding thermal bath.
Furthermore, the longitudinal coupling is preferred to enhance heat current in
strong qubit-resonator coupling regime. For two-photon correlation function,
it exhibits an antibunching-to-bunching transition, which is mainly dominated
by the modulation of energy gap between the first and second excited
eigenstates. Our results are expected to deepen the understanding of
nonequilibrium thermal transport and nonclassical photon radiation based on
the circuit quantum electrodynamics platform.

##  I Introduction

Deep understanding and efficient characterization of nonequilibrium energy
transport via quantum light-matter interaction constitutes an active frontier
for quantum optics and quantum transport gchen2005book ; hxu2016nature ;
plodahl2017nature ; aronzani2018 ; dwwang2019np . The heat flow is considered
as one generic feature of quantum thermal transport, which is bounded by the
second law of thermodynamics. Under thermodynamic bias (e.g., voltage and
temperature bias), the heat current is driven directionally from the hot
source to the cold drain. However, the direction of the current can be
reversed against the thermodynamic bias, e.g., by quantum correlations
kmicadei2019nc  and geometric-phase-induced pump jren2010prl ; zwang2021fop .

Due to the dramatic advance of quantum circuit technology, the circuit quantum
electrodynamics (cQED) systems emerge as one representative platform to
realize quantum light-matter interaction gkurizki2015pnas ; ab2020np ;
aaclerk2020np ; ab2021rmp . The cQED systems are traditionally described by
the seminal quantum Rabi model (QRM) rabi1 ; rabi2 ; braak2011prl ;
qhchen2012pra , i.e., one two-level qubit transversely interacting with a
single-model photon resonator, which is able to describe ultrastrong qubit-
resonator coupling pfdiaz2019rmp ; afkockum2019nrp ; lb2020aqt . QRM has been
extensively investigated in finite-component quantum phase transition
mjhwang2015prl ; mxliu2017prl ; yyzhang2021prl ; mlcai2021nc , quantum
nonlinear optics lgarziano2015pra ; lgarziano2016prl ; xwang2017pra ;
vmacri2021arxiv , and quantum thermodynamics jhwang2012pre ; faltintas2015pra
; sseah2018pre . In particular for two-photon statistics, Ridolfo _et al._
aridolfo2012prl ; aridolfo2013prl ; rstassi2013prl ; lgarziano2017acs
proposed a modified definition of two-photon correlation function to properly
characterize photon nonclassicality at strong qubit-photon coupling.
Meanwhile, an alternative scheme, i.e., longitudinal coupling between the
qubit and the resonator, can also be realized based on the superconducting
circuit engineering tniemczyk2010np ; fyoshihara2017np , which has pronounced
consequences for nonclassical-photon-state generation yjzhao2015pra ;
xinwang2016pra , scalable circuit design pmb2015prb ; sricher2016prb , and
fast nondemolition qubit readout ndidier2015prl ; algrimsmo2019prb .

Recently, quantum thermal transport in cQED systems has attracted increasing
attention, which leads to a flurry of valuable works aronzani2018np ;
jsenior2020cp ; mmajland2020prb ; jppekola2021arxiv ; ylu2021arxiv .
Particularly, Pekola _et al._ aronzani2018np  experimentally detected heat
flow in a hybrid quantum system comprising one transmon-qubit and two
microwave resonators, of which the resonators are individually coupled to two
metallic resistors, respectively. Consequently, the typical thermal
functionalities were realized, e.g., heat valve aronzani2018np , thermal diode
jsenior2020cp , and thermal transistor mmajland2020prb . Moreover, Maguire _et
al._ jsmith2014pra ; jsmith2016jcp ; hmaguire2019prl  theoretically
investigated Franco-Condon physics in noncommutative QRM via the reaction
coordinate mapping approach. Yamamoto _et al._ tyamamoto2021jpcm  unravelled
nontrivial two-peak feature of thermal conductance in QRM at linear response
limit. Wang _et al._ cwang2021cpl ; cwang2021cpb  analyzed nonmonotonic
behavior of the heat current in a longitudinally-coupled qubit-resonator
model. However, nonequilibrium heat flow and the microscopic picture in
dissipative QRM at finite temperature bias currently is lack of exploration,
which is crucial to deepen the understanding of nonequilibrium heat transport
based on the cQED platform. Furthermore, though both the transverse and
longitudinal couplings, between the qubit and the optical resonator can be
flexibly modulated tniemczyk2010np ; fyoshihara2017np ; lgarziano2016prl ;
xinwang2016pra ; nlambert2018prb , the influence of composite qubit-resonator
interaction in the heat current and two-photon statistics has not yet been
reported.

In this paper, we applied dressed master equation (DME) to investigate the
effect of composite qubit-resonator interaction on quantum thermal transport.
At weak qubit-resonator coupling, it is found the heat current exhibits
monotonic enhancement by increasing thermal bath temperature bias in the
dissipative QRM, which stems from cotunneling transport processes. In
contrast, the current is changed to show nonmonotonic feature with
longitudinal qubit-resonator coupling, which is dominated by cyclic heat
exchange transitions. These two distinct microscopic processes are crucial to
unravel the physics pictures of quantum thermal transport in dissipative cQED
systems. At strong qubit-resonator coupling with different composite angles,
the current generally shows nomonotonic behavior, which tightly relies on the
scattering processes between the qubit and bosons in the corresponding thermal
bath. Furthermore, the optimal composite angle to generate the maximal heat
current, gradually changes from the transverse coupling type to the
longitudinal counterpart with increase of qubit-resonator coupling strength.
We also study steady-state two-photon correlation function at strong qubit-
resonator coupling. It is found that by tuning up the composite angle, an
antibunching-to-bunching transition is significantly exhibited, which mainly
originates from increment and reduction of the energy gap between the first
and second excited states. In particular, a giant bunching signature of
photons is unravelled at moderate composite angle.

The rest of this paper is organized as follows: In section II, we present the
composite qubit-resonator model, derive the dressed master equation, and
obtain the expression of the steady-state heat current. In section III, we
investigate the effect of the composite angle on the steady-state heat current
and two-photon correlation function. The microscopic pictures of these
behaviors are also discussed. Finally, we give a summary in section IV.

Figure 1: (Color online) (a) The schematic description of a composite qubit-
resonator model. The red half-circle (top-left) is the optical resonator,
where $\hat{a}$ annihilates one photon with the frequency $\omega_{0}$. The
blue circle (top-right) represents the qubit, which is characterized by the
Pauli operator $\hat{\sigma}_{z}$ and the splitting energy $\varepsilon$. The
double-arrowed rectangle shows the composite qubit-resonator interaction with
the interaction strength $\lambda$ and composite angle $\theta$. The red
(bottom-left) and blue (bottom-right) rectangles denote two thermal baths,
characterized by temperatures $T_{\textrm{R}}$ and $T_{\textrm{Q}}$, and
bosonic annihilators $\hat{b}_{k,\textrm{R}}$ and $\hat{b}_{k,\textrm{Q}}$.
(b) Steady-state heat current $J_{ss}/(\alpha\omega_{0})$ modulated by
temperature bias ${\Delta}T/\omega_{0}$ and the composite angle $\theta$, with
weak qubit-resonator interaction strength $\lambda=0.01\omega_{0}$. Other
parameters are given by $\varepsilon=1.5\omega_{0}$, $\alpha=0.001$,
$\omega_{c}=10\omega_{0}$, $T_{\textrm{R}}=\omega_{0}+{\Delta}T/2$, and
$T_{\textrm{Q}}=\omega_{0}-{\Delta}T/2$. (c) and (d) show cotunneling
processes for first terms of the components $I_{x,1}$ and $I_{x,2}$ in Eqs.
(14a-14b), respectively. (e) and (f) describe cyclic energy exchange processes
for first terms of current components $I_{z,1}$ and $I_{z,2}$ in Eqs.
(16a-16b), respectively. The red solid (blue dashed) arrowed lines denote
transitions between two eigenstates assisted by the R-th (Q-th) thermal bath.

##  II Model and method

###  II.1 Composite qubit-resonator model

We study the dissipative qubit-resonator hybrid model in Fig. 1(a), where one
two-level qubit shows both longitudinal and transverse couplings to a single-
mode resonator, each individually interacting with the corresponding thermal
baths. The Hamiltonian is described as ($\hbar=1$ and $k_{B}=1$)

$\displaystyle~{}\hat{H}=\hat{H}^{\theta}_{\textrm{S}}+\sum_{\mu=\textrm{Q},\textrm{R}}(\hat{H}^{\mu}_{\textrm{B}}+\hat{V}_{\mu}).$
(1)

Specifically, the composite qubit-resonator system is expressed as

$\displaystyle~{}\hat{H^{\theta}_{\textrm{S}}}=\frac{\varepsilon}{2}\hat{\sigma}_{z}+\omega_{0}\hat{a}^{\dagger}\hat{a}+\lambda(\cos\theta\hat{\sigma}_{x}+\sin\theta\hat{\sigma}_{z})(\hat{a}^{\dagger}+\hat{a}),$
(2)

where $\hat{a}^{\dagger}(\hat{a})$ is the creation (annihilation) operator of
one photon in the resonator with the frequency $\omega_{0}$, $\varepsilon$ is
the splitting energy of the two-level qubit,
$\hat{\sigma}_{\alpha}~{}(\alpha=x,y,z)$ is the Pauli operators of the qubit
under the qubit basis $\\{\uparrow,\downarrow\\}$, with
$\hat{\sigma}_{z}|\uparrow{\rangle}=|\uparrow{\rangle}$ and
$\hat{\sigma}_{z}|\downarrow{\rangle}=-|\downarrow{\rangle}$, $\lambda$ is the
qubit-resonator coupling strength, and $\theta$ is the angle to tune the
composite qubit-resonator interaction.
$\hat{H}^{\mu}_{\textrm{B}}=\sum_{k}\omega_{k,\mu}\hat{b}^{\dagger}_{k,\mu}\hat{b}_{k,\mu}$
describes the $\mu$-th bosonic thermal bath, where
$\hat{b}^{\dagger}_{k,\mu}~{}(\hat{b}_{k,\mu})$ creates (annihilates) one
boson with the frequency $\omega_{k,\mu}$ and momentum $k$. The system-bath
interactions are described as

$\displaystyle\hat{V}_{\textrm{R}}=$
$\displaystyle(\hat{a}^{\dagger}+\hat{a})\sum_{k}g_{k,\textrm{R}}(\hat{b}^{\dagger}_{k,\textrm{R}}+\hat{b}_{k,\textrm{R}}),~{}$
(3a) $\displaystyle\hat{V}_{\textrm{Q}}=$
$\displaystyle\hat{\sigma}_{x}\sum_{k}g_{k,\textrm{Q}}(\hat{b}_{k,\textrm{Q}}^{\dagger}+\hat{b}_{k,\textrm{Q}}),~{}$
(3b)

with $g_{k,\mu}$ the coupling strength. The $\mu$-th thermal bath is
characterized by the spectral function
$\gamma_{\mu}(\omega)=2\pi\sum_{k}|g_{k,\mu}|^{2}\delta(\omega-\omega_{k,\mu})$.
In this work, we specify the spectral function as the Ohmic case, i.e.,
$\gamma_{\mu}(\omega)=\pi\alpha\omega\exp(-|\omega|/\omega_{c})$, with
$\alpha$ the dissipation strength and $\omega_{c}$ the cutoff frequency.

Generally, it is difficult to analytically find the eigensoultion of the
qubit-resonator hybrid system at Eq. (2). However, in the limiting case we may
obtain the analytical expression. Specifically, at $\theta=0$ the composite
system with only transverse coupling is simplified as

$\displaystyle\hat{H}^{0}_{\textrm{S}}=\frac{\varepsilon}{2}\hat{\sigma}_{z}+\omega_{0}\hat{a}^{\dagger}\hat{a}+\lambda\hat{\sigma}_{x}(\hat{a}^{\dagger}+\hat{a}),$
(4)

which is the seminal quantum Rabi model rabi1 ; rabi2 . Hence, the eigenvalues
can be mapped to the roots of transcendental G-function by applying the
Bargmann algebra and extended coherent states approaches braak2011prl ;
qhchen2012pra , respectively. In weak qubit-resonator coupling regime, the
quantum Rabi model is reduced to the Jaynes-Cummings model
$\hat{H}_{\textrm{JC}}={\varepsilon}\hat{\sigma}_{z}/2+\omega_{0}\hat{a}^{\dagger}\hat{a}+\lambda(\hat{a}^{\dagger}\hat{\sigma}_{-}+\hat{a}\hat{\sigma}_{+})$,
which is dominated by the rotating-wave-terms. Consequently, the eigenvalues
are given by

$\displaystyle
E_{n,\pm}=(n+1/2)\omega_{0}{\pm}\sqrt{(\varepsilon-\omega_{0})^{2}/4+\lambda^{2}(n+1)},$
(5)

and the corresponding eigenvectors are

$\displaystyle|\phi^{+}_{n}{\rangle}=$
$\displaystyle\cos\frac{\theta_{n}}{2}|n,\uparrow{\rangle}+\sin\frac{\theta_{n}}{2}|n+1,\downarrow{\rangle},~{}$
(6a) $\displaystyle|\phi^{-}_{n}{\rangle}=$
$\displaystyle-\sin\frac{\theta_{n}}{2}|n,\uparrow{\rangle}+\cos\frac{\theta_{n}}{2}|n+1,\downarrow{\rangle},~{}$
(6b)

with $\tan\theta_{n}=2\lambda\sqrt{n+1}/(\varepsilon-\omega_{0})$.

While at $\theta=\pi/2$ the qubit is longitudinally coupled to the resonator
yjzhao2015pra ; sricher2016prb ; xinwang2016pra ; pmb2015prb ; nlambert2018prb
; cwang2021cpl , with the Hamiltonian

$\displaystyle\hat{H}^{\pi/2}_{\textrm{S}}=\frac{\varepsilon}{2}\hat{\sigma}_{z}+\omega_{0}\hat{a}^{\dagger}\hat{a}+\lambda\hat{\sigma}_{z}(\hat{a}^{\dagger}+\hat{a}).$
(7)

Accordingly, the eigenvalues are expressed as

$\displaystyle E_{n,\uparrow}=$
$\displaystyle\omega_{0}n+\varepsilon/2-\lambda^{2}/\omega_{0},$ (8a)
$\displaystyle E_{n,\downarrow}=$
$\displaystyle\omega_{0}n-\varepsilon/2-\lambda^{2}/\omega_{0},$ (8b)

and the eigenvectors are described by the extended coherent boson states

$\displaystyle|\phi^{\uparrow}_{n}{\rangle}=$
$\displaystyle\exp\Big{[}\frac{\lambda}{\omega_{0}}(\hat{a}-\hat{a}^{\dagger})\Big{]}\frac{(\hat{a}^{\dagger})^{n}}{\sqrt{n!}}|0{\rangle}_{a}{\otimes}|\uparrow{\rangle},~{}$
(9a) $\displaystyle|\phi^{\downarrow}_{n}{\rangle}=$
$\displaystyle\exp\Big{[}-\frac{\lambda}{\omega_{0}}(\hat{a}-\hat{a}^{\dagger})\Big{]}\frac{(\hat{a}^{\dagger})^{n}}{\sqrt{n!}}|0{\rangle}_{a}{\otimes}|\downarrow{\rangle},~{}$
(9b)

with the vacuum state of the resonator $\hat{a}|0{\rangle}_{a}=0$.

We should note that the dissipative composite qubit-resonator model is more
than a toy model. It can be experimentally realized in the superconducting
quantum circuit platforms aronzani2018np ; jsenior2020cp ; jppekola2021arxiv ;
mmajland2020prb , where the transmon qubit is able to show both the
longitudinal and transverse interactions with the microwave resonator, and the
bosonic thermal bath is simulated by the LC circuit coupled to a metallic
resistor. Hence, the heat energy naturally flows from the hot source to the
cold drain under the temperature bias.

###  II.2 Quantum master equation

We consider weak interactions between the composite hybrid system and bosonic
thermal baths. We focus on the steady-state properties of the composite qubit-
resonator system, where the off-diagonal elements of the reduced system
density operator in the eigenbasis of $\hat{H}_{\textrm{S}}$ become
negligible. Thus, $\hat{V}_{\textrm{R}}$ and $\hat{V}_{\textrm{Q}}$ at Eqs.
(3a-3b) can be properly perturbed. Under the Born approximation, the total
density operator can be separated as
$\hat{\rho}_{\textrm{tot}}(t){\approx}\hat{\rho}_{\textrm{S}}(t){\otimes}\hat{\rho}_{\textrm{B},\textrm{R}}{\otimes}\hat{\rho}_{\textrm{B},\textrm{Q}}$,
where $\hat{\rho}_{\textrm{S}}(t)$ is the reduced hybrid system density
operator, and
$\hat{\rho}_{\textrm{B},\mu}=\exp(-\hat{H}^{\mu}_{\textrm{B}}/k_{\textrm{B}}T_{u})/\textrm{Tr}_{\textrm{B}}\\{\exp(-\hat{H}^{\mu}_{\textrm{B}}/k_{\textrm{B}}T_{u})\\}~{}(\mu=\textrm{R},\textrm{Q})$
is the density operator of the $\mu$-th thermal bath, with $k_{\textrm{B}}$
the Planck constant and $T_{\mu}$ is the temperature of the $\mu$-th bath.
Then, by further including the Markov approximation, we obtain the quantum
dressed master equation as fbeaudoin2011pra ; asettineri2018pra ; lb2020aqt

$\displaystyle~{}\frac{d}{dt}\hat{\rho}_{\textrm{S}}(t)$ $\displaystyle=$
$\displaystyle
i[\hat{\rho}_{\textrm{S}}(t),\hat{H}^{\theta}_{\textrm{S}}]+\sum_{n,m,\mu}\\{\Gamma^{+}_{\mu}(E_{nm})\hat{\mathcal{L}}_{nm}[\hat{\rho}_{\textrm{S}}(t)]$
(10)
$\displaystyle+\Gamma^{-}_{\mu}(E_{nm})\hat{\mathcal{L}}_{mn}[\hat{\rho}_{\textrm{S}}(t)]\\},$

where the dissipator is given by
$\hat{\mathcal{L}}_{nm}[\hat{\rho}_{\textrm{S}}(t)]=|\phi_{n}{\rangle}{\langle}\phi_{m}|\hat{\rho}_{\textrm{S}}(t)|\phi_{m}{\rangle}{\langle}\phi_{n}|-(|\phi_{m}{\rangle}{\langle}\phi_{m}|\hat{\rho}_{\textrm{S}}(t)+\hat{\rho}_{\textrm{S}}(t)|\phi_{m}{\rangle}{\langle}\phi_{m}|)/2$,
with
$\hat{H}^{\theta}_{\textrm{S}}|\phi_{n}{\rangle}=E_{n}|\phi_{n}{\rangle}$, and
the corresponding transitions rates are given by

$\displaystyle\Gamma^{+}_{\mu}(E_{n,m})=$
$\displaystyle\gamma_{\mu}(E_{n,m})n_{\mu}(E_{n,m})|{\langle}\phi_{n}|\hat{A}_{\mu}|\phi_{m}{\rangle}|^{2},~{}$
(11a) $\displaystyle\Gamma^{-}_{\mu}(E_{n,m})=$
$\displaystyle\gamma_{\mu}(E_{nm})[1+n_{\mu}(E_{n,m})]|{\langle}\phi_{n}|\hat{A}_{\mu}|\phi_{m}{\rangle}|^{2},~{}$
(11b)

with $n_{\mu}(E_{n,m})=1/[\exp(E_{n,m}/T_{u})-1]$ the Bose-Einstein
distribution function, $\hat{A}_{\textrm{R}}=\hat{a}^{\dagger}+\hat{a}$,
$\hat{A}_{\textrm{Q}}=\hat{\sigma}_{x}$, $E_{n,m}=E_{n}-E_{m}$ the energy gap
between two eigenstates $|\phi_{n}{\rangle}$ and $|\phi_{m}{\rangle}$ of
$\hat{H}_{\textrm{S}}$. The rate $\Gamma^{+(-)}_{\mu}(E_{n,m})$ describes the
energy exchange process that the composite qubit-resonator system is excited
(relaxed) from the eigenstate $|\phi_{m}{\rangle}$ to $|\phi_{n}{\rangle}$ by
absorbing (releasing) one photon with the energy $E_{n}-E_{m}$ from (into) the
$\mu$-th thermal bath.

Therefore, after long time evolution, i.e. $d\hat{\rho}_{\textrm{S}}(t)/dt=0$,
we can obtain the steady-steady population distribution $P_{n}$ with
$P_{n}={\langle}\phi_{n}|\hat{\rho}_{\textrm{S}}(t\rightarrow\infty)|\phi_{n}{\rangle}$.
Moreover, from the dressed master equation Eq. (10) we obtain the steady-state
heat current into the Q-th bath (see part A of the Appendix)

$\displaystyle
J_{ss}=\sum_{E_{n}>E_{n^{\prime}}}E_{n,n^{\prime}}[\Gamma^{-}_{\textrm{Q}}(E_{n,n^{\prime}})P_{n}-\Gamma^{+}_{\textrm{Q}}(E_{n,n^{\prime}})P_{n^{\prime}}],$
(12)

where the energy gap is $E_{n,n^{\prime}}=E_{n}-E_{n^{\prime}}$.

##  III Results and discussions

###  III.1 Steady state heat current

####  III.1.1 Weak qubit-resonator interaction

We first investigate the steady-state behavior of heat current at weak qubit-
resonator coupling in Fig. 1(b), which is both modulated by the temperature
bias ${\Delta}T=T_{\textrm{R}}-T_{\textrm{Q}}$ and composite angle $\theta$.
It is found that for small $\theta$, the heat current exhibits monotonic
increase by increasing the temperature bias ${\Delta}T$, particularly in the
limit of $\theta=0$, i.e., the dissipative quantum Rabi model. However, in the
large $\theta$ regime the heat current is changed to shown nonmonotonic
behavior, i.e., the current is first enhanced and later suppressed with the
increase of ${\Delta}T$, which identifies the signature of the negative
differential thermal conductance (NDTC) bli2006apl ; bli2012rmp ; dhhe2009prb
; dhhe2010pre ; hkchan2014pre . The appearance of NDTC is consistent with
previous works cwang2021cpl ; cwang2021cpb . It needs to note that though not
shown here, the similar result can also be found at resonance
($\varepsilon=\omega_{0}$). Hence, we conclude that the composite qubit-
resonator interaction strongly affects the steady-state heat current.

Then, we try to explore microscopic processes of heat transport with weak
qubit-resonator coupling. We admit that to analytically find the microscopic
mechanism with arbitrary composite angle $\theta$ is quite difficult. Here, we
focus on two limits, i.e., $\theta=0$ and $\theta=\pi/2$, to unravel the
representative physical pictures of the heat current. For $\theta=0$, under
the eigenbasis $\\{|\phi^{\pm}_{n}{\rangle}\\}$ of QRM in Eqs. (6a-6b), the
leading order of steady-state heat current at finite energy bias regime [e.g.,
$(\varepsilon-\omega_{0}){\gg}\lambda$] can be analytically expressed as (see
the detail in part B of the Appendix)

$\displaystyle~{}J^{x}_{ss}$ $\displaystyle{\approx}$
$\displaystyle\frac{\lambda^{2}}{(\varepsilon-\omega_{0})^{2}}(\omega_{0}I_{x,1}+{\varepsilon}I_{x,2}),$
(13)

where two components are specified as

$\displaystyle I_{x,1}=$
$\displaystyle\gamma_{\textrm{Q}}(\omega_{0})\\{n_{\textrm{R}}(\omega_{0})[1+n_{\textrm{Q}}(\omega_{0})]-[1+n_{\textrm{R}}(\omega_{0})]n_{\textrm{Q}}(\omega_{0})\\},~{}$
(14a) $\displaystyle I_{x,2}=$
$\displaystyle\frac{\gamma_{\textrm{R}}(\varepsilon)}{2n_{\textrm{Q}}(\varepsilon)+1}\\{n_{\textrm{R}}(\varepsilon)[1+n_{\textrm{Q}}(\varepsilon)]-{[1+n_{\textrm{R}}(\varepsilon)]}n_{\textrm{Q}}(\varepsilon)\\}.~{}$
(14b)

Both $I_{x,1}$ and $I_{x,2}$ are dominated by cotunneling processes.
Specifically, $I_{x,1}$ describes the process such that as the state
$|\phi^{\eta}_{n}{\rangle}$ $(|\phi^{\eta}_{n+1}{\rangle})$ is excited
(relaxed) to $|\phi^{\eta}_{n+1}{\rangle}$ $(|\phi^{\eta}_{n}{\rangle})$ by
absorbing (emitting) energy $\omega_{0}$ from (into) the R-th reservoir, the
transition $|\phi^{\eta}_{n+1}{\rangle}{\rightarrow}|\phi^{\eta}_{n}{\rangle}$
$(|\phi^{\eta}_{n}{\rangle}{\rightarrow}|\phi^{\eta}_{n+1}{\rangle})$
simultaneously occurs by emitting (absorbing) $\omega_{0}$ into (from) the
Q-th reservoir with $\eta=\pm$, which is also shown in Fig. 1(c). While
$I_{x,2}$ shows other typical cotunneling processes, exemplified in Fig. 1(d),
that the excitation (relaxation) transition
$|\phi^{-(+)}_{n}{\rangle}{\rightarrow}|\phi^{+(-)}_{n}{\rangle}$ by absorbing
(emitting) energy $\varepsilon$ from (into) the R-th reservoir is accompanied
by the dual transition
$|\phi^{+(-)}_{n}{\rangle}{\rightarrow}|\phi^{-(+)}_{n}{\rangle}$. It is
interesting to find that the directional cotunneling transport from the R-th
reservoir to Q-th one, i.e. described by first terms of $I_{x,1}$ and
$I_{x,2}$, are monotonically enhanced with increase of the temperature bias,
which mainly contribute to $J^{x}_{ss}$ at finite temperature bias. In
contrast, the opposite transitions (from the Q-th reservoir to R-th one) is
dramatically suppressed. Finally, the current in Eq. (13) becomes significant
at large temperature bias
$(T_{\textrm{R}}{\approx}2\omega_{0},~{}T_{\textrm{Q}}{\approx}0)$, which is
specified as
$J^{x}_{ss}{\approx}[{\lambda}/{(\varepsilon-\omega_{0})}]^{2}[\omega_{0}\gamma_{\textrm{Q}}(\omega_{0})n_{\textrm{R}}(\omega_{0})+{\varepsilon}\gamma_{\textrm{R}}(\varepsilon)n_{\textrm{R}}(\varepsilon)]$.

While for $\theta=\pi/2$, it is found that the transition coefficient
${\langle}\phi^{\sigma}_{n}|\hat{\sigma}_{x}|\phi^{\overline{\sigma}}_{n^{\prime}}{\rangle}$
in Eqs. (11a-11b) under the coherent-state basis
$\\{|\phi^{\uparrow(\downarrow)}_{n}{\rangle}\\}$ is approximated as
${\langle}\phi^{\uparrow}_{n}|\hat{\sigma}_{x}|\phi^{\downarrow}_{n^{\prime}}{\rangle}{\approx}(-1)^{n}\Big{[}\delta_{n,n^{\prime}}+({2\lambda}/{\omega_{0}})(\sqrt{n+1}\delta_{n,n^{\prime}-1}-\sqrt{n}\delta_{n,n^{\prime}+1})\Big{]}$
(also see Refs. qhchen2008pra ; jren2012prb ). Based on the systematic
perturbation (see the full solution in part C of the Appendix), the steady-
state heat current in the bias regime $[(\varepsilon-\omega_{0}){\gg}\lambda]$
is described as

$\displaystyle~{}J^{z}_{ss}\approx\Big{(}\frac{2\lambda}{\omega_{0}}\Big{)}^{2}\omega_{0}(I_{z,1}+I_{z,2}),$
(15)

where these two components are given by

$\displaystyle I_{z,1}=$
$\displaystyle\frac{\gamma_{\textrm{Q}}(\varepsilon+\omega_{0})}{2n_{\textrm{Q}}(\varepsilon)+1}\\{[1+n_{\textrm{Q}}(\varepsilon+\omega_{0})]n_{\textrm{Q}}(\varepsilon)n_{\textrm{R}}(\omega_{0})$
$\displaystyle-
n_{\textrm{Q}}(\varepsilon+\omega_{0})[1+n_{\textrm{Q}}(\varepsilon)][1+n_{\textrm{R}}(\omega_{0})]\\},~{}$
(16a) $\displaystyle I_{z,2}=$
$\displaystyle\frac{\gamma_{\textrm{Q}}(\varepsilon-\omega_{0})}{2n_{\textrm{Q}}(\varepsilon)+1}\\{n_{\textrm{Q}}(\varepsilon-\omega_{0})[1+n_{\textrm{Q}}(\varepsilon)]n_{\textrm{R}}(\omega_{0})$
$\displaystyle-[1+n_{\textrm{Q}}(\varepsilon-\omega_{0})]n_{\textrm{Q}}(\varepsilon)[1+n_{\textrm{R}}(\omega_{0})]\\}.~{}$
(16b)

In sharp contrast to $\theta=0$ limit, $I_{z,1}$ is contributed by two
competing cyclic fluxes, rather not the cotunneling processes. Specifically,
the first term
$[1+n_{\textrm{Q}}(\varepsilon+\omega_{0})]n_{\textrm{Q}}(\varepsilon)n_{\textrm{R}}(\omega_{0})$
describes the loop transition
$|\phi^{\uparrow}_{n+1}{\rangle}{\rightarrow}|\phi^{\downarrow}_{n}{\rangle}{\rightarrow}|\phi^{\uparrow}_{n}{\rangle}{\rightarrow}|\phi^{\uparrow}_{n+1}{\rangle}$
by directionally transferring the energy $\omega_{0}$ into the Q-th thermal
reservoir, which is also depicted in Fig. 1(c). And the second term
$n_{\textrm{Q}}(\varepsilon+\omega_{0})[1+n_{\textrm{Q}}(\varepsilon)][1+n_{\textrm{R}}(\omega_{0})]$
shows the counter loop transition. Similarly, $I_{z,2}$ is composed by other
two opposite cyclic fluxes, where the first loop path is shown in Fig. 1(d),
and characterizes the joint transport process
$|\phi^{\uparrow}_{n+1}{\rangle}{\rightarrow}|\phi^{\downarrow}_{n+1}{\rangle}{\rightarrow}|\phi^{\uparrow}_{n}{\rangle}{\rightarrow}|\phi^{\uparrow}_{n+1}{\rangle}$.
Intriguingly, at large temperature bias ${\Delta}T{\approx}2\omega_{0}$, i.e.,
$T_{\textrm{R}}{\approx}2\omega_{0}$ and $T_{\textrm{Q}}{\approx}0$, All of
cyclic current components in $I_{z,1}$ and $I_{z,2}$ break down, due to
negligible excitation in the Q-th reservoir
[$n_{\textrm{Q}}(\omega>0){\approx}0$]. This directly results in the
suppression of the steady-state heat current, which identifies the emergence
of the NDTC.

Therefore, we exploit two distinct microscopic pictures in limiting composite
angles with weak qubit-resonator interaction, i.e., cotunneling transitions at
$\theta=0$ and cyclic transitions at $\theta=\pi/2$, which are generic to
unravel microscopic mechanisms of quantum thermal transport, e.g., in
nonequilibrium spin-boson model truokola2011prb ; cwang2015sr ; dzxu2016fop ;
cwang2017pra ; tchen2020pre ; crduan2020jpcl , metal-insulator interfaces
jren2013prb ; jren2013prbr , and inelastic thermoelectrics jhjiang2015prb ;
jhjiang2017prapp ; jclu2020prb ; jclu2021prb . The analytical expressions of
the heat current in Eq. (13) and Eq. (15) are obtained for the first time in
dissipative qubit-resonator hybrid model. Moreover, though the transport
picture with longitudinal qubit-resonator coupling was preliminarily reported
in Ref. cwang2021cpl , it was analyzed based on only transition rates.

####  III.1.2 Strong qubit-resonator interaction

Figure 2: (Color online) Steady-state heat current $J_{ss}/(\alpha\omega_{0})$
modulated by temperature bias ${\Delta}T/\omega_{0}$ and qubit-resonator
interaction strength $\lambda/\omega_{0}$ with composite angle (a) $\theta=0$,
(b) $\theta=\pi/4$, and (c) $\theta=\pi/2$. Other parameters are given by
$\varepsilon=1.5\omega_{0}$, $\alpha=0.001$, $\omega_{c}=10\omega_{0}$,
$T_{\textrm{R}}=\omega_{0}+{{\Delta}T}/2$, and
$T_{\textrm{Q}}=\omega_{0}-{{\Delta}T}/2$.  Figure 3: (Color online) Steady-
state heat current $J_{ss}/(\alpha\omega_{0})$ modulated by qubit-resonator
interaction strength $\lambda/\omega_{0}$ and composite angle $\theta$. The
dashed line with red circles shows the maximal $J_{ss}/(\alpha\omega_{0})$ by
searching $\theta{\in}[0,\pi/2]$ with given $\lambda/\omega_{0}$. Other
parameters are given by $\varepsilon=1.5\omega_{0}$, $\alpha=0.001$,
$\omega_{c}=10\omega_{0}$, $T_{\textrm{R}}=2\omega_{0}$, and
$T_{\textrm{Q}}=0$.

Next, we investigate steady-state heat current $J_{ss}/(\alpha\omega_{0})$
beyond weak qubit-resonator coupling with typical composite angles in Fig. 2.
As the qubit shows transverse interaction with the optical resonator, i.e.,
$\theta=0$ in Fig. 2(a), the heat current exhibits nonmonotonic behavior by
increasing qubit-resonator interaction strength, under finite temperature bias
(e.g., ${\Delta}T/\omega_{0}=1$). Due to the effect of counter-rotating-terms,
the eigenstates of the quantum Rabi model beyond Eqs. (6a-6b) will introduce
additional energy exchange transitions, which may effectively enhance heat
current in the regime $\lambda/\omega_{0}{\lesssim}0.5$. While in strong
qubit-resonator coupling limit (e.g., $\lambda/\omega_{0}>1.5$), the
eigenstates become nearly degenerate, i.e.,
$|\phi^{\pm}_{n}{\rangle}{\approx}\exp[\pm\lambda(\hat{a}-\hat{a}^{\dagger})/\omega_{0}][(\hat{a}^{\dagger})^{n}/\sqrt{n!}]|0{\rangle}{\otimes}|\pm{\rangle}$
with $\hat{\sigma}_{x}|\pm{\rangle}=\pm|\pm{\rangle}$, which significantly
prohibits energy exchange between the hybrid system and the Q-th reservoir
(${\langle}\phi^{+}_{n}|\hat{\sigma}_{x}|\phi^{-}_{n}{\rangle}{\approx}0$).
Then, by tuning on the composite angle, e.g., $\theta=\pi/4$ and $\pi/2$, it
is found that the profiles of heat currents modulated by
${\Delta}T/\omega_{0}$ and $\lambda/\omega_{0}$ in Figs. 2(b-c) are generally
similar with limiting angle case $\theta=0$. In particular for the limiting
case $\theta=\pi/2$, the transition coefficient in Eqs. (11a-11b) is expressed
as qhchen2008pra ; jren2012prb

$\displaystyle~{}{\langle}\phi^{\uparrow}_{n}|\hat{\sigma}_{x}|\phi^{\downarrow}_{n^{\prime}}{\rangle}$
$\displaystyle=$
$\displaystyle{(-1)}^{n}\exp{(-2\lambda^{2}/\omega^{2}_{0})}\sqrt{n!n^{\prime}!}$
(17)
$\displaystyle{\times}\sum^{\min[n,n^{\prime}]}_{l=0}\frac{(-1)^{l}\sqrt{(2\lambda/\omega_{0})^{n+n^{\prime}-2l}}}{(n-l)!(n^{\prime}-l)!l!},$

which induces higher-order transitions between
$|\phi^{\uparrow(\downarrow)}_{n}{\rangle}$ and
$|\phi^{\downarrow(\uparrow)}_{n^{\prime}}{\rangle}$ with
$|n-n^{\prime}|{\geq}2$, besides the lowest-order transport processes, i.e.,
cyclic exchange in Figs. 1(c-d). These additional transitions are robust even
at large temperature bias, which mainly result in comparatively large heat
current and disappearance of the NDTC. While at finite temperature bias (e.g.,
${\Delta}T/\omega_{0}=1$), the initial enhancement of
$J_{ss}/(\alpha\omega_{0})$ by increasing $\lambda/\omega_{0}$ stems from
additional transport processes. It may be quantified by
${\langle}\phi^{\uparrow}_{n}|\hat{\sigma}_{x}|\phi^{\downarrow}_{n^{\prime}}{\rangle}{\approx}(-1)^{n}\Big{\\{}[1-(n+1/2)(\frac{2\lambda}{\omega_{0}})^{2}]\delta_{n,n^{\prime}}+(\frac{2\lambda}{\omega_{0}})(\sqrt{n+1}\delta_{n,n^{\prime}-1}-\sqrt{n}\delta_{n,n^{\prime}+1})+\frac{1}{2}(\frac{2\lambda}{\omega_{0}})^{2}[\sqrt{n(n-1)}\delta_{n,n^{\prime}+2}+\sqrt{(n+1)(n+2)}\delta_{n,n^{\prime}-2}]\Big{\\}}$,
where last two terms will enhance the current by forming efficient transition
paths. However, the final decrease of heat current in strong qubit-resonator
coupling limit is mainly attributed to the fact the transition coefficient
${\langle}\phi^{\uparrow}_{n}|\hat{\sigma}_{x}|\phi^{\downarrow}_{n^{\prime}}{\rangle}$
in Eq. (17) and the corresponding transition rates are dramatically weakened.
Consequently, the energy exchange processes are strongly blocked.

Moreover, we analyze the interplay between the composite angle and qubit-
resonator coupling strength on
$\max_{\\{\theta\\}}\\{J_{ss}/(\alpha\omega_{0})\\}$ in Fig. 3. It is
intriguing to find that with increase of the qubit-resonator interaction, the
composite angle dominating the maximal heat current is gradually modified from
$0$ to $\pi/2$, corresponding to the transverse and longitudinal qubit-
resonator couplings, respectively. Hence, we conclude that at strong
$\lambda/\omega_{0}$, the longitudinal coupling type is preferred to enhance
the steady-state heat current.

###  III.2 Two-photon correlation function

Two-photon correlation function describes the correlation between two
temporally separated photon signals from one light source, which is pioneered
by R. J. Glauber to unveiling the optical coherence of quantum theory
rjglauber1963pr . Later, it has been extended to investigate superradiant
spontaneous emission jheberly1970pra ; rbonifacio1971pra ; rbonifacio1971pra2
; dmeiser2010pra ; AAuffeves2011njp ; EMascarenhas2013pra , strongly
interacting photons aimamoglu1997prl ; prabl2011prl ; droy2017rmp , and photon
(phonon) bundle emission qbin2020prl ; qbin2021prl . Alternatively, Ridolfo
_et al._ aridolfo2012prl ; aridolfo2013prl ; rstassi2013prl ; lgarziano2017acs
proposed a modified definition of two-photon correlation function within the
dressed framework. Hence, it can also be safely included to study photon
statistics in the present model in Eq. (1) with strong qubit-resonator
interaction. Specifically, two-photon correlation function at steady state is
defined as

$\displaystyle
G^{(2)}_{\theta}(\tau)=\lim_{t{\rightarrow}\infty}\frac{{\langle}\hat{X}^{+}_{\theta}(t)\hat{X}^{+}_{\theta}(t+\tau)\hat{X}^{-}_{\theta}(t+\tau)\hat{X}^{-}_{\theta}(t){\rangle}}{{\langle}\hat{X}^{+}_{\theta}(t)\hat{X}^{-}_{\theta}(t){\rangle}^{2}},$
(18)

where the measurement operator is
$\hat{X}^{-}_{\theta}=-i\sum_{k>j}\Delta_{kj}X_{jk}|\phi_{j}{\rangle}{\langle}\phi_{k}|$
and $\hat{X}^{+}_{\theta}=(\hat{X}^{-}_{\theta})^{\dagger}$, with
$X_{jk}={\langle}\phi_{j}|(\hat{a}^{\dagger}+\hat{a})|\phi_{k}{\rangle}$ the
energy gap $\Delta_{jk}=E_{j}-E_{k}$, and the eigensolution
$\hat{H}^{\theta}_{\textrm{S}}|\phi_{k}{\rangle}=E_{k}|\phi_{k}{\rangle}$.
Here, we focus on the effect of the composite qubit-resonant interaction on
zero-time two-photon correlation function $G^{(2)}_{\theta}(0)$.

Figure 4: (Color online) Two-photon correlation function as functions of
qubit-resonator interaction strength $\lambda/\omega_{0}$ and temperature
$T_{\textrm{R}}=T_{\textrm{Q}}=T$, with typical composite angle (a)
$\theta=0$, (b) $\theta=\pi/4$, and (c) $\theta=\pi/2$. Other parameters are
given by $\varepsilon=1.5\omega_{0}$, $\alpha=0.001$, and
$\omega_{c}=10\omega_{0}$.  Figure 5: (Color online) (a) Two-photon
correlation function $G^{(2)}_{\theta}(0)$ modulated by qubit-resonator
interaction strength $\lambda/\omega_{0}$ and composite angle $\theta$ in low
temperature regime. (b) Comparison of $G^{(2)}_{\theta}(0)$ via numerical
calculation with that based on Eq. (*), (c) four lowest eigenenergies, (d)
four lowest steady-state populations, and (e) coefficients $A_{n},~{}B_{n}$ as
a function of composite angle $\theta$ at $\lambda=\omega_{0}$. Other
parameters are given by $\varepsilon=1.5\omega_{0}$, $\alpha=0.001$,
$\omega_{c}=10\omega_{0}$, and $T_{\textrm{R}}=T_{\textrm{Q}}=0.1\omega_{0}$.

We first investigate two-photon correlation function by tuning both the qubit-
resonator coupling strength and bath temperature with different composite
angle in Fig. 4. For $\theta=0$, it is found that an intriguing antibunching
behavior can be found with strong $\lambda$, which demonstrates the seminal
two-photon blockade aridolfo2012prl . Then, by tuning on composite angle,
e.g., $\theta=\pi/4$, a giant photon bunching behavior is exhibited at low
temperature. If we further increase the composite angle to $\theta=\pi/2$, it
is shown that $G^{(2)}_{\theta=0}(0){\approx}2$ regardless of $\lambda$ and
$T$ ($T_{\textrm{R}}=T_{\textrm{Q}}=T$), due to the fully thermalization of
the longitudinally coupled qubit-resonator system

$\displaystyle~{}\rho_{\textrm{S}}(\infty)$ $\displaystyle=$
$\displaystyle\frac{\sinh[\omega_{0}/(2T)]}{\cosh[\varepsilon/(2T)]}\sum_{n}e^{-\frac{(n+1/2)\omega_{0}}{T}}$
(19)
$\displaystyle{\times}[e^{-\frac{\varepsilon}{2T}}|\phi^{\uparrow}_{n}{\rangle}{\langle}\phi^{\uparrow}_{n}|+e^{\frac{\varepsilon}{2T}}]|\phi^{\downarrow}_{n}{\rangle}{\langle}\phi^{\downarrow}_{n}|].$

Moreover, we plot Fig. 5(a) to see the influence of the composite angle in
two-photon correlation function. It is found that at strong $\lambda$ (e.g.,
$\lambda/\omega_{0}{\approx}1$), an antibunching-to-bunching transition is
clearly exhibited by increasing $\theta$. Therefore, we conclude that the
modulation of the composite angle is quite important to exhibit the
nonclassical photon statistics, which may provide physical guidance to measure
photon correlation in circuit QED.

Next, we analyze the mechanism of antibunching-to-bunching transition
modulated by the composite angle. At low temperature (e.g.,
$T=0.1\omega_{0}$), the finite spacing distribution of energy levels [see Fig.
5(b)] results in $P_{0}{\gg}P_{1}{\gg}P_{2}{\gg}P_{3}$, exhibited in Fig.
5(c). Hence, one-photon and two-photon terms are approximated as
${\langle}\hat{X}^{+}_{\theta}\hat{X}^{-}_{\theta}{\rangle}{\approx}P_{1}A_{1}$
and
${\langle}(\hat{X}^{+}_{\theta})^{2}(\hat{X}^{-}_{\theta})^{2}{\rangle}{\approx}P_{2}B_{2}$,
where the coefficients are $A_{n}=\sum_{l<k}(\Delta_{kl}X_{kl})^{2}$ and
$B_{n}=\sum_{p<l<k}(\Delta_{kl}\Delta_{lp}X_{kl}X_{lp})^{2}$. Consequently,
the two-photon correlation function is expressed in a concise way

$\displaystyle~{}G^{(2)}_{\theta}(0){\approx}\frac{P_{2}B_{2}}{(P_{1}A_{1})^{2}},$
(20)

which shows agreement with the numerical result in Fig. 5(d) in a wide regime
of $\theta$ at strong qubit-resonator coupling ($\lambda/\omega_{0}=1$).
Particularly in the anti-bunching regime (e.g., $0<\theta{{\lesssim}}\pi/10$),
the large energy gap $(E_{2}-E_{1})$ suppresses the ratio $P_{2}/P^{2}_{1}$.
Moreover, the coefficients $A_{1}$ and $B_{2}$ are nearly flat, as shown in
Fig. 5(e). Hence, the transition
$|\phi_{1}{\rangle}{\rightarrow}|\phi_{2}{\rangle}$ is strongly blocked,
leading to the antibunching behavior of photons. While in the bunching regime
(e.g., $\pi/10{\lesssim}\theta{\leq}\pi/4$), though $A_{1}$ is strengthened by
increasing $\theta$, the reduction of energy gap between $E_{2}$ and $E_{1}$
dramatically enhances steady-state population $P_{2}$ and successive two-
photon excitation process
$|\phi_{0}{\rangle}{\rightarrow}|\phi_{1}{\rangle}{\rightarrow}|\phi_{2}{\rangle}$,
resulting in the bunching behavior of photons, as shown in Fig. 5(b).

##  IV Conclusion

In summary, we investigate the effect of composite qubit-resonator interaction
on quantum thermal transport and zero-delay-delay two-photon correlation
function at steady state. We apply the quantum dressed master equation to
properly treat strong qubit-resonator interaction with arbitrary composite
angle. For heat transport at weak qubit-resonator coupling, it is found that
the heat current with transverse qubit-resonator coupling shows monotonic
behavior by increasing bath temperature bias. It is dominated by the
cotunneling process, which is quantified by Eqs. (14a-14b). While the current
is gradually changed to exhibit nonmonotonic feature by tuning on the
composite angle, signifying the emergence of the NDTC effect, completely
characterized as the cyclic energy exchange processes, and described by Eqs.
(16a-16b). Hence, we unravel two crucial microscopic processes for quantum
thermal transport. In strong qubit-resonator interaction regime, the heat
current with $\theta=0$ exhibits nonmonotonic behavior by increasing qubit-
resonator interaction strength. The initial enhancement of $J_{ss}$ stems from
additional energy transitions due to counter-rotating-terms, whereas the final
suppression of $J_{ss}$ is attributed to the nearly degeneracy of eigenstates,
which prohibits energy exchange between the qubit and Q-th thermal reservoir.
Moreover, it is intriguing to find that the optimal composite angle, which
corresponds to the maximal heat current, switches directionally from
$\theta=0$ (transverse) to $\theta=\pi/2$ (longitudinal) by modulating qubit-
resonator interaction from weak to strong couplings. Hence, the longitudinal
coupling is preferred to enhance the steady-state heat current in strong
qubit-resonator interaction regime.

We also investigate the steady-state two-photon correlation function by
modulating the composite angle. At $\theta=0$, the pronounced anti-bunching
feature is exhibited in the regime of low temperature and strong qubit-
resonator interaction. Then by tuning on $\theta$, we find the giant bunching
signal instead. By further increasing $\theta$ (e.g., $\theta=\pi/2$), two-
photon correlation function is globally around $2$, due to special
distribution of the density operator of the qubit-resonator hybrid system in
Eq. (19). Moreover, We also present the mechanism with approximate expression
in Eq. (20) to explain this antibunching-to-bunching transition. The
antibunching and bunching behaviors of photons are modulated by the large and
reduced energy gap between the first and second excited eigenstates,
respectively. We hope that our results affected by composite qubit-resonator
interaction may deepen the understanding of quantum thermal transport and two-
photon statistics in dissipative QED systems.

## ACKNOWLEDGEMENTS

Z.C., H.C, Z.C, and C.W. 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.R. 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).

## Appendix

###  IV.1 General expression of heat current

From the quantum dressed master equation Eq. (10), the population dynamics is
described as

$\displaystyle~{}\frac{d}{dt}\rho_{nn}(t)=\sum_{n,m,\mu}[\Gamma^{-}_{\mu}(E_{n,m})\rho_{mm}(t)-\Gamma^{+}_{\mu}(E_{n,m})\rho_{nn}(t)],$
(21)

with the transition rates $\Gamma^{\pm}_{\mu}(E_{n,m})$ specified by Eqs.
(11a-11b), and the population
$\rho_{nn}(t)={\langle}\phi_{n}|\hat{\rho}_{\textrm{S}}(t)|\phi_{n}{\rangle}$.
Hence, the expression of steady-state heat current into the $\mu$-th thermal
bath is expressed as

$\displaystyle
J_{\mu}=\sum_{E_{n}>E_{n^{\prime}}}E_{n,n^{\prime}}[\Gamma^{-}_{\mu}(E_{n,n^{\prime}})P_{nn}-\Gamma^{+}_{\mu}(E_{n,n^{\prime}})P_{n^{\prime}n^{\prime}}].$
(22)

Then, we try to approximately investigate the steady state populations
$P_{nn}={\langle}\phi_{n}|\hat{\rho}_{\textrm{S}}(t\rightarrow\infty)|\phi_{n}{\rangle}$
at weak qubit-resonator coupling limit via the perturbation method. As we
reexpress the populations in the vector form $|\rho_{s}{\rangle}$, the steady-
state solution based on Eq. (21) becomes $\mathcal{L}|\rho_{s}{\rangle}=0$. At
weak qubit-resonator coupling, we expand
$|\rho_{\textrm{S}}{\rangle}{\approx}|\rho^{(0)}_{\textrm{S}}{\rangle}+(\lambda/\omega_{0})^{2}|\rho^{(1)}_{\textrm{S}}{\rangle}$
and
$\mathcal{L}{\approx}\mathcal{L}^{(0)}+(\lambda/\omega_{0})^{2}\mathcal{L}^{(1)}$
up to $(\lambda/\omega_{0})^{2}$. Hence, the general solution is given by

$\displaystyle\mathcal{L}^{(0)}|\rho^{(0)}_{\textrm{S}}{\rangle}=$
$\displaystyle 0,~{}$ (23a)
$\displaystyle\mathcal{L}^{(0)}|\rho^{(1)}_{\textrm{S}}{\rangle}+\mathcal{L}^{(1)}|\rho^{(0)}_{\textrm{S}}{\rangle}=$
$\displaystyle 0~{}.$ (23b)

###  IV.2 $\theta=0$ case

At $\theta=0$, the model in the weak qubit-resonator coupling regime becomes
the dissipative Jaynes-Cummings model. We consider the off-resonant regime
$[(\varepsilon-\omega_{0}){\gg}\lambda]$. Under the eigenbasis
$|\phi^{\pm}_{n}{\rangle}$ [Eqs. (6a-6b)] the transition coefficients involved
with the Q-th bath are specified as

$\displaystyle{\langle}\phi^{+}_{n+1}|\hat{\sigma}_{x}|\phi^{+}_{n}{\rangle}=$
$\displaystyle\sin\frac{\varphi_{n}}{2}\cos\frac{\varphi_{n+1}}{2},$ (24a)
$\displaystyle{\langle}\phi^{+}_{n+1}|\hat{\sigma}_{x}|\phi^{-}_{n}{\rangle}=$
$\displaystyle\cos\frac{\varphi_{n}}{2}\cos\frac{\varphi_{n+1}}{2},$ (24b)
$\displaystyle{\langle}\phi^{-}_{n+1}|\hat{\sigma}_{x}|\phi^{+}_{n}{\rangle}=$
$\displaystyle-\sin\frac{\varphi_{n}}{2}\sin\frac{\varphi_{n+1}}{2},$ (24c)
$\displaystyle{\langle}\phi^{-}_{n+1}|\hat{\sigma}_{x}|\phi^{-}_{n}{\rangle}=$
$\displaystyle-\cos\frac{\varphi_{n}}{2}\sin\frac{\varphi_{n+1}}{2},$ (24d)

with $\tan\varphi_{n}=2\lambda\sqrt{n+1}/\Delta$ and
$\Delta=\varepsilon-\omega_{0}$. Here, we further approximately treat the
eigenvalues as $E_{n,+}{\approx}\omega_{0}(n+1/2)+\varepsilon/2$ and
$E_{n,-}{\approx}\omega_{0}(n+3/2)-\varepsilon/2$. And the corresponding
eigenvectors are simplified to
$|\phi^{+}_{n}{\rangle}{\approx}|n{\rangle}{\otimes}|\uparrow{\rangle}$ and
$|\phi^{-}_{n}{\rangle}{\approx}|n+1{\rangle}{\otimes}|\downarrow{\rangle}$.
Hence, the transition rates upper to the order $(\lambda/\Delta)^{2}$ are
given by

$\displaystyle\Gamma^{\pm}_{\textrm{Q}}(E^{n+1,+}_{n,+})\approx$
$\displaystyle\gamma_{\textrm{Q}}(\pm\omega_{0})n_{\textrm{Q}}(\pm\omega_{0})\Big{(}\frac{\lambda}{\Delta}\Big{)}^{2}(n+1),$
(25a) $\displaystyle\Gamma^{\pm}_{\textrm{Q}}(E^{n+1,+}_{n,-})\approx$
$\displaystyle\gamma_{\textrm{Q}}(\pm\varepsilon)n_{\textrm{Q}}(\pm\varepsilon)\Big{[}1-\Big{(}\frac{\lambda}{\Delta}\Big{)}^{2}(2n+3)\Big{]},$
(25b) $\displaystyle\Gamma^{\pm}_{\textrm{Q}}(E^{n+1,-}_{n,+})\approx$
$\displaystyle 0,$ (25c)
$\displaystyle\Gamma^{\pm}_{\textrm{Q}}(E^{n+1,-}_{n,-})\approx$
$\displaystyle\gamma_{\textrm{Q}}(\pm\omega_{0})n_{\textrm{Q}}(\pm\omega_{0})\Big{(}\frac{\lambda}{\Delta}\Big{)}^{2}(n+2),$
(25d)

with
$E^{n,\sigma}_{n^{\prime},\sigma^{\prime}}=E_{n,\sigma}-E_{n^{\prime},\sigma^{\prime}}$.
Similarly, the rates related with the R-th bath are approximated by

$\displaystyle\Gamma^{\pm}_{\textrm{R}}(E^{n+1,+}_{n,+})\approx$
$\displaystyle\gamma_{\textrm{R}}(\pm\omega_{0})n_{\textrm{R}}(\pm\omega_{0})(n+1)\Big{[}1+\Big{(}\frac{\lambda}{\Delta}\Big{)}^{2}\Big{]},$
(26a) $\displaystyle\Gamma^{\pm}_{\textrm{R}}(E^{n+1,+}_{n,-})\approx$
$\displaystyle\gamma_{\textrm{Q}}(\pm\varepsilon)n_{\textrm{Q}}(\pm\varepsilon)\Big{(}\frac{\lambda}{\Delta}\Big{)}^{2},$
(26b) $\displaystyle\Gamma^{\pm}_{\textrm{R}}(E^{n,-}_{n-1,+})\approx$
$\displaystyle 0,$ (26c)
$\displaystyle\Gamma^{\pm}_{\textrm{R}}(E^{n,-}_{n-1,-})\approx$
$\displaystyle\gamma_{\textrm{Q}}(\pm\omega_{0})n_{\textrm{Q}}(\pm\omega_{0})(n+2)\Big{[}1-\Big{(}\frac{\lambda}{\Delta}\Big{)}^{2}\Big{]}.$
(26d)

From Eq. (23b) we have the relation at steady state

$\displaystyle\sum_{n}[\Gamma^{-}_{\textrm{Q}}(E^{n,+}_{n-1,-})P_{n,+}-\Gamma^{+}_{\textrm{Q}}(E^{n,+}_{n-1,-})P_{n-1,-}]$
$\displaystyle=$
$\displaystyle\sum_{n}[\Gamma^{+}_{\textrm{R}}(E^{n,+}_{n-1,-})P_{n-1,-}-\Gamma^{-}_{\textrm{R}}(E^{n,+}_{n-1,-})P_{n,+}]$
$\displaystyle+\sum_{n}[\Gamma^{+}_{\textrm{Q}}(E^{n,+}_{n-1,+})P_{n-1,+}-\Gamma^{-}_{\textrm{Q}}(E^{n,+}_{n-1,+})P_{n,+}]$
$\displaystyle+\sum_{n}[\Gamma^{-}_{\textrm{Q}}(E^{n+1,+}_{n,+})P_{n+1,+}-\Gamma^{+}_{\textrm{Q}}(E^{n+1,+}_{n,+})P_{n,+}].$

Hence, the zeroth order populations are given by

$\displaystyle P^{(0)}_{n,+}\approx$
$\displaystyle\frac{(1-e^{-\beta_{\textrm{R}}\omega_{0}})}{e^{\beta_{\textrm{Q}}\varepsilon}+1}e^{-{n}\beta_{\textrm{R}}\omega_{0}},$
(27a) $\displaystyle P^{(0)}_{n,-}\approx$
$\displaystyle\frac{e^{\beta_{\textrm{Q}}\varepsilon}(1-e^{-\beta_{\textrm{R}}\omega_{0}})}{e^{\beta_{\textrm{Q}}\varepsilon}+1}e^{-{(n+1)}\beta_{\textrm{R}}\omega_{0}},$
(27b)

with $\beta_{u}=1/(T_{u})~{}(u=\textrm{R},\textrm{Q})$. Consequently, we
obtain the expression of heat current as

$\displaystyle J^{x}_{\textrm{Q}}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{\lambda}{\varepsilon-\omega_{0}}\Big{)}^{2}(\omega_{0}I_{x,1}+{\varepsilon}I_{x,2}),$
(28)

where two components are specified as

$\displaystyle I_{x,1}=$
$\displaystyle\gamma_{\textrm{Q}}(\omega_{0})\\{n_{\textrm{R}}(\omega_{0})[1+n_{\textrm{Q}}(\omega_{0})]-[1+n_{\textrm{R}}(\omega_{0})]n_{\textrm{Q}}(\omega_{0})\\},$
(29a) $\displaystyle I_{x,2}=$
$\displaystyle\frac{\gamma_{\textrm{R}}(\varepsilon)}{2n_{\textrm{Q}}(\varepsilon)+1}\\{n_{\textrm{R}}(\varepsilon)[1+n_{\textrm{Q}}(\varepsilon)]-{[1+n_{\textrm{R}}(\varepsilon)]}n_{\textrm{Q}}(\varepsilon)\\},$
(29b)

###  IV.3 $\theta=\pi/2$ case

At $\theta=\pi/2$, Based on the eigenbasis
$|\phi^{\uparrow(\downarrow)}_{n}{\rangle}$ at Eqs. (9a-9b) the transition
coefficient involved with the Q-th bath is simplified to

$\displaystyle{\langle}\phi^{\sigma}_{n}|\hat{\sigma}_{x}|\phi^{\overline{\sigma}}_{n^{\prime}}{\rangle}$
$\displaystyle{\approx}$
$\displaystyle(-1)^{n}\Big{[}\delta_{n,n^{\prime}}+\frac{2\lambda}{\omega_{0}}\sqrt{n+1}\delta_{n,n^{\prime}-1}$
(30)
$\displaystyle-\frac{2\lambda}{\omega_{0}}\sqrt{n}\delta_{n,n^{\prime}+1}\Big{]}.$

Accordingly, the transition rates defined by Eqs. (11a-11b) are approximated
by

$\displaystyle\Gamma^{\pm}_{\textrm{Q}}(E^{n,\uparrow}_{n^{\prime},\downarrow})\approx$
$\displaystyle\delta_{n,n^{\prime}}\kappa^{\pm}_{\textrm{Q}}(\varepsilon)+\delta_{n,n^{\prime}-1}n^{\prime}\Big{(}\frac{2\lambda}{\omega_{0}}\Big{)}^{2}\kappa^{\pm}_{\textrm{Q}}(\varepsilon-\omega_{0})$
(31a)
$\displaystyle+\delta_{n,n^{\prime}+1}n\Big{(}\frac{2\lambda}{\omega_{0}}\Big{)}^{2}\kappa^{\pm}_{\textrm{Q}}(\varepsilon+\omega_{0}),$
$\displaystyle\Gamma^{\pm}_{\textrm{Q}}(E^{n,\downarrow}_{n^{\prime},\uparrow})\approx$
$\displaystyle\delta_{n,n^{\prime}+1}n^{\prime}\Big{(}\frac{2\lambda}{\omega_{0}}\Big{)}^{2}\kappa^{\pm}_{\textrm{Q}}(\omega_{0}-\varepsilon),$
(31b)

with
$\kappa^{+}_{\textrm{Q}}(\omega)=\gamma_{\textrm{Q}}(\omega)n_{\textrm{Q}}(\omega)$
and
$\kappa^{-}_{\textrm{Q}}(\omega)=\gamma_{\textrm{Q}}(\omega)[1+n_{\textrm{Q}}(\omega)]$.
Similarly, the transition rates assisted by the R-th thermal bath are given by
$\Gamma^{+}_{\textrm{R}}(E^{n,\sigma}_{n-1,\sigma})=\gamma_{\textrm{R}}(\omega_{0})n_{\textrm{R}}(\omega_{0})n$
and
$\Gamma^{-}_{\textrm{R}}(E^{n,\sigma}_{n-1,\sigma})=\gamma_{\textrm{R}}(\omega_{0})[1+n_{\textrm{R}}(\omega_{0})]n$.
Hence, the zeroth order of populations based on Eq. (23a) can be directly
obtained as

$\displaystyle P^{(0)}_{n,\uparrow}\approx$
$\displaystyle\frac{(1-e^{-\beta_{\textrm{R}}\omega_{0}})}{e^{\beta_{\textrm{Q}}\varepsilon}+1}e^{-{n}\beta_{\textrm{R}}\omega_{0}},$
(32a) $\displaystyle P^{(0)}_{n,\downarrow}\approx$
$\displaystyle\frac{e^{\beta_{\textrm{Q}}\varepsilon}(1-e^{-\beta_{\textrm{R}}\omega_{0}})}{e^{\beta_{\textrm{Q}}\varepsilon}+1}e^{-{n}\beta_{\textrm{R}}\omega_{0}}.$
(32b)

Moreover, from Eq. (23b) it is known that

$\displaystyle\Big{(}\frac{2\lambda}{\omega_{0}}\Big{)}^{2}\sum_{n}[\Gamma^{-}_{\textrm{Q}}(E^{n,\uparrow}_{n,\downarrow})P^{(1)}_{n,\uparrow}-\Gamma^{+}_{\textrm{Q}}(E^{n,\uparrow}_{n,\downarrow})P^{(1)}_{n\downarrow}]\approx$
$\displaystyle\sum_{n}[-\Gamma^{-}_{\textrm{Q}}(E^{n,\uparrow}_{n-1,\downarrow})P^{(0)}_{n,\uparrow}+\Gamma^{+}_{\textrm{Q}}(E^{n,\uparrow}_{n-1,\downarrow})P^{(0)}_{n-1,\downarrow}$
$\displaystyle-\Gamma^{-}_{\textrm{Q}}(E^{n,\uparrow}_{n+1,\downarrow})P^{(0)}_{n,\uparrow}+\Gamma^{+}_{\textrm{Q}}(E^{n,\uparrow}_{n+1,\downarrow})P^{(0)}_{n+1,\downarrow}$
$\displaystyle-\Gamma^{+}_{\textrm{Q}}(E^{n+1,\downarrow}_{n,\uparrow})P^{(0)}_{n+1,\downarrow}+\Gamma^{-}_{\textrm{Q}}(E^{n+1,\downarrow}_{n,\uparrow})P^{(0)}_{n,\uparrow}].$

Finally, Then, the current is contributed by three components

$\displaystyle
J^{z}_{\textrm{Q}}\approx\Big{(}\frac{2\lambda}{\omega_{0}}\Big{)}^{2}\omega_{0}(I_{z,1}+I_{z,2}+I_{z,3}),$
(33)

where

$\displaystyle I_{z,1}=$
$\displaystyle\theta(\varepsilon+\omega_{0})\gamma_{\textrm{Q}}(\varepsilon+\omega_{0})\frac{1}{2n_{\textrm{Q}}(\varepsilon)+1}$
$\displaystyle{\times}[(1+n_{\textrm{Q}}(\varepsilon+\omega_{0}))n_{\textrm{Q}}(\varepsilon)n_{\textrm{R}}(\omega_{0})$
$\displaystyle-
n_{\textrm{Q}}(\varepsilon+\omega_{0})(1+n_{\textrm{Q}}(\varepsilon))(1+n_{\textrm{R}}(\omega_{0}))],$
(34a) $\displaystyle I_{z,2}=$
$\displaystyle\theta(\varepsilon-\omega_{0})\gamma_{\textrm{Q}}(\varepsilon-\omega_{0})\frac{1}{2n_{\textrm{Q}}(\varepsilon)+1}$
$\displaystyle{\times}[n_{\textrm{Q}}(\varepsilon-\omega_{0})(1+n_{\textrm{Q}}(\varepsilon))n_{\textrm{R}}(\omega_{0})$
$\displaystyle-(1+n_{\textrm{Q}}(\varepsilon-\omega_{0}))n_{\textrm{Q}}(\varepsilon)(1+n_{\textrm{R}}(\omega_{0}))],$
(34b) $\displaystyle I_{z,3}=$
$\displaystyle\theta(\omega_{0}-\varepsilon)\gamma_{\textrm{Q}}(\omega_{0}-\varepsilon)\frac{1}{2n_{\textrm{Q}}(\varepsilon)+1}$
$\displaystyle{\times}[(1+n_{\textrm{Q}}(\omega_{0}-\varepsilon))(1+n_{\textrm{Q}}(\varepsilon))n_{\textrm{R}}(\omega_{0})$
$\displaystyle-
n_{\textrm{Q}}(\omega_{0}-\varepsilon)n_{\textrm{Q}}(\varepsilon)(1+n_{\textrm{R}}(\omega_{0}))].$
(34c)

with the Heviside function $\theta(x)=1$ for $x{\geq}0$, and $\theta(x)=0$ for
$x{<}0$.

## References

  * (1) G. Chen, _Nanoscale energy transport and conversion: a parallel treatment of electrons, molecules, phonons, and photons_ (Oxford University Press, 2005).
  * (2) H. Xu, D. Mason, L. Y. Jinag, and J. G. E. Harris, Nature 537, 80 (2016).
  * (3) P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeute, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, Nature 541, 473 (2017).
  * (4) A. Ronzani, B. Karimi, J. Senior, Y. C. Chang, J. T. Peltonen, C. D. Chen, and J. P. Pekola, Nat. Phys. 14, 991 (2018).
  * (5) D. W. Wang, C. Song, W. Feng, H. Cai, D. Xu, H. Deng, H. K. Li, D. N. Zheng, X. B. Zhu, H. Wang, S. Y. Zhu, and M. O. Scully, Nat. Phys. 15, 382 (2019).
  * (6) K. Micadei, J. P. S. Peterson, A. M. Souza, R. S. Sarthour, Ivan S. Oliveira, G. T. Landi, T. B. Batalhao, R. M. Serra, and E. Lutz, Nat. Commun. 10, 2456 (2019).
  * (7) J. Ren, P. Hanggi, and B. Li, Phys. Rev. Lett. 104, 170601 (2010).
  * (8) Z. Wang, L. Q. Wang, J. Z. Wang, C. Wang, and J. Ren, Frontiers of Physics, 17, 13201 (2022).
  * (9) G. Kurizki, P. Bertet, Y. Kubo, K. Molmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, PNAS 112, 3866 (2015).
  * (10) A. Blais, S. M. Girvin, and W. D. Oliver, Nat. Phys. 16, 247 (2020).
  * (11) A. A. Clerk, K W. Lehnert, P. Bertet, J. R. Petta, and Y. Nakamura, Nat. Phys. 16, 257 (2020).
  * (12) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * (13) I. I. Rabi, Phys. Rev. 49, 324 (1936).
  * (14) I. I. Rabi, Phys. Rev. 51, 652 (1937).
  * (15) D. Braak, Phys. Rev. Lett. 107, 100401 (2011).
  * (16) Q. H. Chen, C. Wang, S. He, T. Liu, and K. L. Wang, Phys. Rev. A 86, 023822 (2012).
  * (17) P. Forn-Díaz, L. Mamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * (18) A. F. Kockum, A. Miranwicz, S. De Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * (19) A. Le Boité, Adv. Quantum Technol. 3, 1900140 (2020).
  * (20) M. J. Hwang, R. Peubla, and M. B. Pleino, Phys. Rev. Lett. 115, 180404 (2015).
  * (21) M. X. Liu, S. Chesi, Z. J. Ying, X. S. Chen, H. G. Luo, and H. Q. Lin, Phys. Rev. Lett. 119, 220601 (2017).
  * (22) Y. Y. Zhang, Z. X. Hu, L. B. Fu, H. G. Luo, H. Pu, and X. F. Zhang, Phys. Rev. Lett. 127, 063602 (2021).
  * (23) M. L. Cai, Z. D. Liu, W. D. Zhao, Y. K. Wu, Q. X. Mei, Y. Jiang, L. He, X. Zhang, Z. C. Zhou, and L. M. Duan, Nat. Commun. 12, 1126 (2021).
  * (24) L. Garziano, R. Stassi, V. Macrí, A. F. Kockum, S. Savasta, and F. Nori, Phys. Rev. A 92, 063830 (2015).
  * (25) L. Garziano, V. Macrí, R. Stassi, O. Di Stefano, F. Nori, and S. Savasta, Phys. Rev. Lett. 117, 04301 (2016).
  * (26) X. Wang, A. Miranowicz, H. R. Li, and F. Nori, Phys. Rev. A 96, 063820 (2017).
  * (27) V. Macrí, F. Minganti, A. F. Kockum, A. Ridolfo, S. Savasta, F. Nori, arXiv:2107.08759.
  * (28) J. H. Wang, Z. Q. Wu, and J. Z. He, Phys. Rev. E 85, 041148 (2012).
  * (29) F. Altintas, A. Ü. C Hardal, and Ö. E. Müstecaplioglu, Phys. Rev. A 91, 023816 (2015).
  * (30) S. Seah, S. Nimmrichter, and V. Scarani, Phys. Rev. E 98, 012131 (2018).
  * (31) A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann, Phys. Rev. Lett. 109, 193602 (2012).
  * (32) A. Ridolfo, S. Savasta, and M. J. Hartmann, Phys. Rev. Lett. 110, 163601 (2013).
  * (33) R. Stassi, A. Ridolfo, O. Di Stefano, M. J. Hartmann, and S. Savasta, Phys. Rev. Lett. 110, 243601 (2013).
  * (34) L. Garziano, A. Ridolfo, S. De Liberato, and S. Savasta, ACS. Photonics 4, 2345 (2017).
  * (35) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, and R. Gross, Nat. Phys. 6, 772 (2010).
  * (36) F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Nat. Phys. 13, 44 (2017).
  * (37) Y. J. Zhao, Y. L. Liu, Y. X. Liu, and F. Nori, Phys. Rev. A 91, 053820 (2015).
  * (38) X. Wang, A. Miranowicz, H. R. Li, and F. Nori, Phys. Rev. A 94, 053858 (2016).
  * (39) P. M. Billangeon, J. S. Tsai, and Y. Nakamura, Phys. Rev. B 91, 094517 (2015).
  * (40) S. Richer and D. DiVincenzo, Phys. Rev. B 93, 134501 (2016).
  * (41) N. Didier, J. Bourassa, and A Blais, Phys. Rev. Lett. 115, 20360 (2015).
  * (42) A. L. Grimsmo and T. B. Smith, Phys. Rev. B 99, 235420 (2019).
  * (43) A. Ronzani, B. Karimi, J. Senior, Y. C. Chang, J. T. Peltonen, C. D. Chen, and J. P. Pekola, Nat. Phys. 14, 991 (2018).
  * (44) J. Senior, A. Gubaydullin, B. Karimi, J. T. Peltonen, J. Ankerhold, and J. P. Pekola, Comm. Phys. 3, 40 (2020).
  * (45) M. Majland, K. S. Christensen, and N. T. Zinner, Phys. Rev. B 101, 184510 (2020).
  * (46) J. P. Pekola and B. Karimi, arXiv:2107.12936.
  * (47) Y. Lu, N. Lambert, A. F. Kockum, K. Funo, A. Bengtsson, S. Gasparinetti, F. Nori, and P. Delsing, arXiv:2107.12700.
  * (48) J. Lles-Smith, N. Lambert, and A. Nazir, Phys. Rev. Lett. 90, 032114 (2014).
  * (49) J. Lles-Smith, A. G. Dijkstra, N. Lambert, and A. Nazir, J. Chem. Phys. 144, 044110 (2016).
  * (50) H. Maguire, J. Lles-Smith, and A. Nazir, Phys. Rev. Lett. 123, 093601 (2019).
  * (51) T. Yamamoto and T. Kato, J. Phys.: Condens. Matter 33, 395303 (2021).
  * (52) C. Wang, L. Q. Wang, and J. Ren, Chin. Phys. Lett. 38, 010501 (2021).
  * (53) C. Wang, L. Q. Wang, and J. Ren, Chin. Phys. B 30, 030506 (2021).
  * (54) N. Lambert, M. Cirio, M. Delbecq, G. Allison, M. Marx, S. Tarucha, and F. Nori, Phys. Rev. B 97, 125429 (2018).
  * (55) F. Beaudoin, J. M. Gambetta, and A. Blais, Phys. Rev. A 84, 043832 (2011).
  * (56) A Settineri, V. Macrí, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta, Phys. Rev. A 98, 053834 (2018).
  * (57) B. Li, L. Wang, and G. Casati, Appl. Phys. Lett. 88, 143501 (2006).
  * (58) N. B. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012).
  * (59) D. H. He, S. Buyukdagli, and B. Hu, Phys. Rev. B 80, 104302 (2009).
  * (60) D. H. He, B. Q. Ai, H. K. Chan, and B. Hu, Phys. Rev. E 81, 041131 (2010).
  * (61) H. K. Chan, D. H. He, and B. Hu, Phys. Rev. E 89, 052126 (2014).
  * (62) Q. H. Chen, Y. Y. Zhang, T. Liu, and K. L. Wang, Phys. Rev. A 78, 051801 (2008).
  * (63) J. Ren, J. X. Zhu, J. E. Gubernatis, C. Wang, and B. Li, Phys. Rev. B 85, 155443 (2012).
  * (64) T. Ruokola and T. Ojanen, Phys. Rev. B 83, 045417 (2011).
  * (65) C. Wang, J. Ren, and J. Cao, Sci. Rep. 5, 11787 (2015).
  * (66) D. Z. Xu and J. Cao, Front. Phys. 11, 110308 (2016).
  * (67) C. Wang, J. Ren, and J. Cao, Phys. Rev. A 95, 023610 (2017).
  * (68) T. Q. Chen, V. Balachandran, C. Guo, and D. Poletti, Phys. Rev. E 102, 012155 (2020).
  * (69) C. R. Du, C. Y. Hsieh, J. J. Liu, J. L. Wu, and J. Cao, J. Phys. Chem. Lett. 11, 4080 (2020).
  * (70) J. Ren and J. X. Zhu, Phys. Rev. B 87, 241412 (2013).
  * (71) J. Ren, Phys. Rev. B 88, 220406(R) (2013).
  * (72) J. H. Jiang, M Kulkarni, D. Segal, and Y. Imry, Phys. Rev. B 92, 045309 (2015).
  * (73) J. H. Jiang and Y. Imry, Phys. Rev. Applied 7, 064001 (2017).
  * (74) J. C. Lu, R. Q. Wang, C. Wang, and J. H. Jiang, Phys. Rev. B 102, 125405 (2020).
  * (75) J. C. Lu, J. H. Jiang, and Y. Imry, Phys. Rev. B 103, 085429 (2021).
  * (76) R. J. Glauber, Phys. Rev. 130, 2529 (1963).
  * (77) J. H. Eberly and N. E. Rehler, Phys. Rev. A 2, 1607 (1970).
  * (78) R. Bonifacio, P. Schwendimann, and Fritz Haake, Phys. Rev. A 4, 302 (1971).
  * (79) R. Bonifacio, P. Schwendimann, and Fritz Haake, Phys. Rev. A 4, 854 (1971).
  * (80) D. Meiser and M. J. Holland, Phys. Rev. A 81, 063827 (2010).
  * (81) A. Aufféves, D. Gerace, S. Portolan, A. Drezet, and M. F. Santos, New J. Phys. 13, 093020 (2011).
  * (82) E. Mascarenhas, D. Gerace, M. F. Santos, and A. Aufféves, Phys. Rev. A 88, 063825 (2013).
  * (83) A. Imamoglu, H. Schmidt, G. Woods, and M. Deutsch, Phys. Rev. Lett. 79, 1467 (1997).
  * (84) P. Rabl, Phys. Rev. Lett. 107, 063601 (2011).
  * (85) D. Roy, C. M. Wilson, and O. Firstenberg, Rev. Mod. Phys. 89, 021001 (2017).
  * (86) Q. Bin, X. Y. Lü, F. P. Laussy, F. Nori, and Y. Wu, Phys. Rev. Lett. 124, 053601 (2020).
  * (87) Q. Bin, Y. Wu, and X. Y. Lü, Phys. Rev. Lett. 127, 073602 (2021).