# Nonequilibrium thermal transport and photon squeezing in a quadratic qubit-
resonator system

Chen Wang1, Hua Chen1<EMAIL_ADDRESS>Jie-Qiao Liao2,
<EMAIL_ADDRESS>1Department of Physics, Zhejiang Normal University, Jinhua
321004, Zhejiang , P. R. China
2Key Laboratory of Low-Dimensional Quantum Structures and Quantum Control of
Ministry of Education, Key Laboratory for Matter Microstructure and Function
of Hunan Province, Department of Physics and Synergetic Innovation Center for
Quantum Effects and Applications, Hunan Normal University, Changsha 410081,
China

###### Abstract

We investigate steady-state thermal transport and photon statistics in a
nonequilibrium hybrid quantum system, in which a qubit is longitudinally and
quadratically coupled to an optical resonator. Our calculations are conducted
with the method of the quantum dressed master equation combined with full
counting statistics. The effect of negative differential thermal conductance
is unravelled at finite temperature bias, which stems from the suppression of
cyclic heat transitions and large mismatch of two squeezed photon modes at
weak and strong qubit-resonator hybridizations, respectively. The giant
thermal rectification is also exhibited at large temperature bias. It is found
that the intrinsically asymmetric structure of the hybrid system and negative
differential thermal conductance show the cooperative contribution. Noise
power and skewness, as typical current fluctuations, exhibit global maximum
with strong hybridization at small and large temperature bias limits,
respectively. Moreover, the effect of photon quadrature squeezing is
discovered in the strong hybridization and low-temperature regime, which shows
asymmetric response to two bath temperatures. These results would provide some
insight to thermal functional design and photon manipulation in qubit-
resonator hybrid quantum systems.

##  I Introduction

The efficient control of energy transport and information processing is of
fundamental importance in low dimensional nano-structures gchen2005book ,
which spurs the advance of interdisciplinary frontiers, ranging from quantum
thermodynamics and thermal machines gbenenti2017pr ; mesposito2019prl ,
quantum electrodynamics (QED) aronzani2018np ; jsenior2020cp ; dwang2019np  to
topological quantum optics hxu2016nature ; sbarik2018sci ; ywang2020pr . Of
particular interest is the nonequilibrium thermal transport in open small
quantum systems, where the degrees of freedom of the target system are much
fewer than those of the external baths uweiss2012book . Accordingly, the
microscopic modeling of the system-bath interaction is crucial for theoretical
characterization of the heat flow. Moreover, the strong system-bath coupling
significantly fertilizes the physical picture of thermal transport
dsegal2006prb ; dsegal2011prb ; ksaito2013prl ; mcarrega2016prl ; wdou2018prb
; anazir2020prl . Hence, it is intriguing to understand the underlying
mechanism of the nonequilibrium thermal transport in small quantum systems.

One of the representative small quantum systems to describe nonequilibrium
thermal transport probably is the nonequilibrium spin-boson model (SBM)
aleggett1987rmp , which is composed by a two-level system interacting with two
individually bosonic thermal baths. The heat flow in the nonequilibrium SBM
has been extensively investigated from weak to strong system-bath couplings by
various analytical approaches, including the Redfield equation dsegal2005prl ,
noninteracting blip approximation dsegal2006prb ; dsegal2011prb ,
nonequilibrium polaron-transformed Redfield equation cwang2015sr ;
cwang2017pra ; xfcao2021prb , nonequilibrium green function yy2014epl ;
jjliu2017pre , and reaction coordinate mapping anazir2020prl . While the rise
of quantum engineering spawns another typical small quantum system platform,
the hybrid quantum system (HQS) mwallquist2009ps ; gkurizki2015pnas ;
aclerk2020np . For example, the circuit-QED system akockum2019nrp ;
pdiaz2019rmp  and optomechanical system maspelmeyer2014rmp , they enable the
strong hybridization of the artificial qubit and mechanical phonon with
optical photon, respectively. Interestingly, the reaction coordinate mapping
tries to theoretically establish the bridge between the SBM and the coupled
qubit-resonator system anazir2014pra ; anazir2020prl ; anazir2016jcp ;
dnewman2017pre ; srestrepo2018njp ; cmc2019jcp . Recently, two-photon based
quadratic interaction in HQSs, e.g., two-photon Rabi model, has attracted
increasing attention, which has practical correspondence in the circuit-QED
platform sfelicetti2018pra1 ; sfelicetti2018pra2 , and theoretically leads to
spectral collapse qhchen2012pra ; lduan2016jpa , quantum phase transition
lgarbe2017pra ; xchen2018pra ; scui2020pra , ultrafast charging quantum
battery ac2020prb ; ad2021entropy  and enhancement of superradiant spontaneous
emission nparxiv2020 . Moreover, such nonlinear interactions play a crucial
role in modeling the electric current based on the quantum dot-vibration
hybrid model jjliu2020nl .

From the aspect of quantum thermal transport, increasing works are carried out
in HQSs. Specially, chiral thermal management aseif2017nc ; zdenis2020prl  and
Berry-phase-like heat pump wjnie2020pra  are proposed in optomechanical
systems, which result in non-reciprocal heat flow. Recently, Wu _et al._
yang2020nc  experimentally analyzed the phonon heat transport in a mechanical
resonator dimer mediated by one optical cavity mode, which directly measured
heat current and proved thermodynamics uncertainty relation. While for the
coupled qubit-resonator system, typical thermal functional devices, e.g.,
quantum thermal transistor and quantum rectifier, have been proposed by
linearly hybridizing the qubit to the optical resonator with the longitudinal
and transverse forms mmajland2020prb ; cwang2021cpl ; cwang2021cpb ;
yjzhao2015pra ; xwang2016pra ; xwang2017pra ; sricher2016prb ; nlambert2018prb
. The novel two-peak feature of linear thermal conductance by tuning the bath
temperature was unravelled, which cannot be simply observed in the SBM
tkato2020arxiv . Meanwhile, the heat flow was practically measured in the
circuit-QED setup with weak qubit-photon hybridization aronzani2018np ;
jsenior2020cp .

In this paper, we study the nonequilibrium thermal transport and photon
squeezing in a longitudinal and quadratic coupling qubit-resonator system.
Concretely, for the steady-state heat current, we clearly unravel the effect
of negative differential thermal conductance (NDTC) in a wide qubit-photon
interaction regime by tuning the bath temperature bias. We also find the giant
thermal rectification at strong qubit-photon hybridization. For the heat
current fluctuations, we analyze the noise power and skewness, which show
global maximum with strong hybridization strength at small and large
temperature biases, respectively. While for the photon squeezing, we find that
strong qubit-photon hybridization and low bath temperatures is helpful to
exhibit the photon squeezing. Moreover, the photon squeezing shows asymmetric
response to two bath temperatures.

The rest of this paper is organized as follows. In Sec. II, we introduce the
coupled qubit-resonator model with quadratic interaction, derive the dressed
master equation (DME) to obtain the steady state, and combine DME with full
counting statistics (FCS) to obtain the heat current fluctuations. In Sec.
III, we first study thermal functionalities, i.e., the NDTC and heat
amplification, by modulating both the qubit-photon hybridization strength and
the bath temperature bias. We also investigate the steady-state noise power
and skewness as the representative current fluctuations. Finally, we analyze
the effect of photon squeezing both at equilibrium and out of equilibrium. In
Sec. IV, we present a conclusion.

##  II Model and method

###  II.1 Qubit-resonator hybrid system

The nonequilibrium hybrid quantum system under consideration in Fig. 1 is
composed by an optical resonator quadratically interacting with a qubit, each
coupled to a bosonic thermal bath. The Hamiltonian of the coupled qubit-
resonator system reads ($\hbar=1$)

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

where $\hat{a}^{\dagger}~{}(\hat{a})$ denotes the creation (annihilation)
operator of the optical resonator with the resonant frequency $\omega_{a}$,
$\sigma_{x}=|1{\rangle}{\langle}0|+|0{\rangle}{\langle}1|$ and
$\sigma_{z}=|1{\rangle}{\langle}1|-|0{\rangle}{\langle}0|$ are the Pauli
operators of the qubit composed by two states $|1{\rangle}$ and $|0{\rangle}$,
and $\lambda$ is the qubit-photon interaction strength. Thereafter, we set
$\omega_{a}=1$ without losing any generality. The $u$th thermal bath is
described as noninteracting bosons
$\hat{H}^{u}_{\textrm{B}}=\sum_{k}\omega_{k,u}\hat{b}^{\dagger}_{k,u}\hat{b}_{k,u}$,
where $\hat{b}^{\dagger}_{k,u}~{}(\hat{b}_{k,u})$ is the
creation(annihilation) operator of the bosonic mode with the momentum $k$ and
the frequency $\omega_{k,u}$. The interactions of the photon and qubit with
the corresponding thermal baths are given by

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

with $g_{k,\textrm{R}(\textrm{Q})}$ being the photon (qubit)-boson interaction
strengthes. Based on the above description, it is known that the Hamiltonian
of the whole system, which both include the hybrid system and the environment,
is described by

$\displaystyle\hat{H}=\hat{H}_{\textrm{S}}+\sum_{u=\textrm{R},\textrm{Q}}(\hat{H}^{u}_{\textrm{B}}+\hat{V}_{u}).$
(3)

Figure 1: (Color online) (a) Schematic of the nonequilibrium coupled qubit-
resonator system. The harmonics and two-level system denote the optical
resonator and qubit, respectively. The pair of double black crossing curves
shows the quadratic qubit-photon interaction; the red and blue panels
represent the bosonic thermal baths, characterized as the temperatures
$T_{\textrm{R}}$ and $T_{\textrm{Q}}$, respectively.

It is interesting to point out that the system Hamiltonian given in Eq. (1)
can be exactly solved. Under the basis states $|1{\rangle}$ and $|0{\rangle}$
of the qubit, the Hamiltonian in Eq. (1) can be decomposed as
$\hat{H}_{\textrm{S}}=\sum_{\sigma=1,0}\hat{H}^{\sigma}_{s}|\sigma{\rangle}{\langle}\sigma|$,
with
$\hat{H}^{\sigma}_{\textrm{S}}=(\omega_{a}+2\lambda_{\sigma})\hat{a}^{\dagger}\hat{a}+\lambda_{\sigma}(\hat{a}^{{\dagger}2}+\hat{a}^{2})+({\varepsilon_{\sigma}}/{2}+\lambda_{\sigma})$,
where the renormalized quantities are given by
$\eta_{\sigma}=\sqrt{1-({2\lambda}/{\omega_{\sigma}})^{2}}$,
$\omega_{\sigma}=\omega_{a}+2\lambda_{\sigma}$,
$\varepsilon_{\sigma}=(-1)^{\sigma+1}\varepsilon$, and
$\lambda_{\sigma}=(-1)^{\sigma+1}\lambda$. Then, by applying the
transformation operator
$\hat{U}=\exp\\{{-\sum_{\sigma=1,0}{\alpha_{\sigma}}/{2}|\sigma{\rangle}{\langle}\sigma|{\otimes}(\hat{a}^{2}-\hat{a}^{{\dagger}2})}\\}$
to $\hat{H}_{\textrm{S}}$ with

$\displaystyle\tanh\alpha_{\sigma}=(-1)^{\sigma+1}\sqrt{\frac{1-\eta_{\sigma}}{1+\eta_{\sigma}}},$
(4)

we find that the transformed Hamiltonian
$\hat{U}\hat{H}_{\textrm{S}}\hat{U}^{\dagger}$ is fully diagonalized. From the
aspect of inverse design, we can express $\hat{H}_{\textrm{S}}$ in a concise
way as

$\displaystyle~{}\hat{H}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle\sum_{\sigma=0,1}[\eta_{\sigma}\omega_{\sigma}\hat{A}^{\dagger}_{\sigma}\hat{A}_{\sigma}+({\eta_{\sigma}-1})\omega_{a}/{2}+{\varepsilon_{\sigma}}/{2}+\lambda_{\sigma}]$
(5) $\displaystyle{\otimes}|\sigma{\rangle}{\langle}\sigma|,$

where the new bosonic operators are specified as
$\hat{A}^{\dagger}_{\sigma}=f_{\sigma}\hat{a}^{\dagger}+g_{\sigma}\hat{a}$ and
$\hat{A}_{\sigma}=f_{\sigma}\hat{a}+g_{\sigma}\hat{a}^{\dagger}$, fulfilling
the commutation relation $[\hat{A}_{\sigma},\hat{A}^{\dagger}_{\sigma}]=1$,
and the coefficients are given by

$\displaystyle~{}f_{\sigma}=\sqrt{\frac{1+\eta_{\sigma}}{2\eta_{\sigma}}},\hskip
14.22636ptg_{\sigma}=(-1)^{\sigma+1}\sqrt{\frac{1-\eta_{\sigma}}{2\eta_{\sigma}}}.$
(6)

Accordingly, the eigen-equation of $\hat{H}_{\textrm{S}}$ is given by
$\hat{H}_{\textrm{S}}|\psi^{\sigma}_{n}{\rangle}=E^{\sigma}_{n}|\psi^{\sigma}_{n}{\rangle}$.
The $n$th eigenvalue is given by
$E^{\sigma}_{n}=\eta_{\sigma}\omega_{\sigma}{n}+[{(\eta_{\sigma}-1)}\omega_{a}/{2}+{\varepsilon_{\sigma}}/{2}+\lambda_{\sigma}]$,
and the corresponding eigenstate is
$|\psi^{\sigma}_{n}{\rangle}=(1/{\sqrt{n!}}){\hat{A}^{{{\dagger}}n}_{{\sigma}}}|\psi^{\sigma}_{0}{\rangle}$,
where the squeezed vacuum state is given by

$\displaystyle|\psi^{\sigma}_{0}{\rangle}=|\sigma{\rangle}{\otimes}(1-|\phi_{\sigma}|^{2})^{1/4}\sum^{\infty}_{n=0}\frac{\phi^{n}_{\sigma}}{2^{n}n!}(\hat{a}^{\dagger})^{2n}|0{\rangle}_{a},$
(7)

with the ratio $\phi_{\sigma}=-g_{\sigma}/f_{\sigma}$ ($|\phi_{\sigma}|<1$)
and the bare vacuum state $\hat{a}|0{\rangle}_{a}=0$. In the squeezing basis,
the frequencies of two squeezing bosonic modes are renormlized by the factor
$\eta_{\sigma}$ and displaced by $\pm 2\lambda_{\sigma}$, respectively. Hence,
in the strongly coupled qubit-photon case, it is straightforward to find that
$\eta_{1}\omega_{1}{\gg}\eta_{0}\omega_{0}$ and the hybridization strength is
bounded by $\lambda<0.25\omega_{a}$.

Recently, quantum heat transport in qubit-resonator hybrid devices have been
theoretically proposed tkato2020arxiv ; mmajland2020prb ; cwang2021cpl ;
cwang2021cpb  and experimentally established aronzani2018np ; jsenior2020cp
based on the circuit-QED setups, where the superconducting qubit can both
longitudinally and transversely interact with the photon resonator
yjzhao2015pra ; xwang2016pra ; xwang2017pra ; sricher2016prb ; nlambert2018prb
. While for the present model Eq. (1) with quadratic interaction, it could
also be realized in the superconducting quantum circuit, e.g., no dc current
biasing the SQUID sfelicetti2018pra1 , where the qubit shows longitudinal
quadratic coupling with the resonator. Moreover, the bosonic thermal baths are
simulated by or the LC circuit coupled to a resistor mmajland2020prb ;
aaclerk2010rmp . Hence, we may be able to analyze heat energy in the circuit-
QED system under the temperature bias.

###  II.2 Quantum master equation

When considering weak system-bath interactions, we are able to separately
perturb the photon-bath and qubit-bath interactions. Under the Born-Markov
approximation, the total density operator can be decomposed as
$\hat{\rho}_{\textrm{tot}}{\approx}\hat{\rho}_{\textrm{S}}{\otimes}\hat{\rho}_{\textrm{B}}$,
where $\hat{\rho}_{\textrm{S}}$ is the density operator of the reduced qubit-
resonator system, and
$\hat{\rho}_{\textrm{B}}=\frac{1}{Z}\Pi_{u=\textrm{R},\textrm{Q}}\exp[-\hat{H}^{u}_{\textrm{B}}/(k_{B}T_{u})]$
is the density operator of thermal baths at equilibrium, with
$Z=\textrm{Tr}_{\textrm{B}}\\{\Pi_{u=\textrm{R},\textrm{Q}}\exp[-\hat{H}^{u}_{\textrm{B}}/(k_{B}T_{u})]\\}$
the partition function, $T_{u}$ the temperature of the $u$th bath, and $k_{B}$
is the Boltzmann constant ($k_{B}$ is set to $1$ for convenience). Then, the
quantum master equation is obtained as (see Appendix)

$\displaystyle~{}\frac{d}{dt}\hat{\rho}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle
i[\hat{\rho}_{\textrm{S}},\hat{H}_{\textrm{S}}]+\frac{1}{2}\sum_{\omega,\omega^{\prime};u=\textrm{R},\textrm{Q}}\\{\kappa_{u}(\omega^{\prime})[\hat{P}_{u}(\omega^{\prime})\hat{\rho}_{s},\hat{P}_{u}(\omega)]$
(8) $\displaystyle+\mathrm{h.c.}\\},$

where the rate is given by
$\kappa_{u}(\omega)=\gamma_{u}(\omega)n_{u}(\omega)$, with
$\gamma_{u}(\omega)=2\pi\sum_{k}|g_{k,u}|^{2}\delta(\omega-\omega_{k})$ the
spectral function and $n_{u}(\omega)=1/[\exp(\omega_{k}/T_{u})-1]$ the Bose-
Einstein distribution function. The spectral function is selected to be super-
Ohmic case
$\gamma_{u}(\omega)=\pi\alpha_{u}({\omega^{3}}/{\omega^{2}_{c}})e^{-\omega/\omega_{c}}\theta(\omega)$,
where $\alpha_{u}$ is the coupling strength, $\omega_{c}$ is the cut-off
frequency, and the heaviside step function is $\theta(x)=1$ for $x{\geq}0$ and
$\theta(x)=0$ for $x<0$. The projecting operator of the resonator originates
from
$[\hat{a}^{\dagger}(-\tau)+\hat{a}(-\tau)]=\sum_{\omega}\hat{P}_{\textrm{R}}(\omega)e^{-i\omega\tau}$,
which is specified as
$\hat{P}_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})=\sum_{\sigma}(f_{\sigma}-g_{\sigma})|\sigma{\rangle}{\langle}\sigma|\hat{A}^{\dagger}_{\sigma}$
and
$\hat{P}_{\textrm{R}}(-\eta_{\sigma}\omega_{\sigma})=\hat{P}^{\dagger}_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})$.
While the projecting operator of the qubit, based on the relation
$\hat{\sigma}_{x}(-\tau)=\sum_{\omega}\hat{P}_{\textrm{Q}}(\omega)e^{-i\omega\tau}$,
is given by
$\hat{P}_{\textrm{Q}}(\omega)=\sum_{n,m,\sigma}{\langle}\phi^{\sigma}_{n}|\hat{\sigma}_{x}|\phi^{\bar{\sigma}}_{m}{\rangle}|\phi^{\sigma}_{n}{\rangle}{\langle}\phi^{\bar{\sigma}}_{m}|\delta(\omega-E^{n,\sigma}_{m,\bar{\sigma}})$,
with the energy gap
$E^{n,\sigma}_{m,\bar{\sigma}}=E_{n,\sigma}-E_{m,\bar{\sigma}}$.

It should be noted that after a long-time evolution, the off-diagonal density
matrix elements of the coupled qubit-resonator system are confirmed to be
negligible. Then, the populations are decoupled from the off-diagonal
elements, which restricts the pair of projectors $\hat{P}_{\mu}(\omega)$ and
$\hat{P}_{\mu}(\omega^{\prime})$ in Eq. (8) to
$\hat{P}_{\mu}(\omega=E^{\sigma}_{n}-E^{\sigma^{\prime}}_{n^{\prime}})={\langle}\psi^{\sigma}_{n}|\hat{A}_{\mu}|\psi^{\sigma^{\prime}}_{n^{\prime}}{\rangle}|\psi^{\sigma}_{n}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{n^{\prime}}|$
and
$\hat{P}_{\mu}(\omega^{\prime}=E^{\sigma^{\prime}}_{n^{\prime}}-E^{\sigma}_{n})={\langle}\psi^{\sigma^{\prime}}_{n^{\prime}}|\hat{A}_{\mu}|\psi^{\sigma}_{n}{\rangle}|\psi^{\sigma^{\prime}}_{n^{\prime}}{\rangle}{\langle}\psi^{\sigma}_{n}|$.
Consequently, quantum master equation can be simplified to the DME
asettineri2018pra

$\displaystyle~{}\frac{d}{dt}\hat{\rho}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{\textrm{S}},\hat{\rho}_{\textrm{S}}]$ (9)
$\displaystyle+\sum_{m,\sigma}\\{\Gamma^{\textrm{R},+}_{m,\sigma}\mathcal{\hat{D}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma}_{m-1}|]\hat{\rho}_{\textrm{S}}$
$\displaystyle+\Gamma^{\textrm{R},-}_{m,\sigma}\mathcal{\hat{D}}[|\psi^{\sigma}_{m-1}{\rangle}{\langle}\psi^{\sigma}_{m}|]\hat{\rho}_{\textrm{S}}\\}$
$\displaystyle+\sum_{m,m^{\prime},\sigma}\\{\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}\mathcal{\hat{D}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\overline{\sigma}}_{m^{\prime}}|]\hat{\rho}_{\textrm{S}}$
$\displaystyle+\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma}\mathcal{\hat{D}}[|\psi^{\overline{\sigma}}_{m^{\prime}}{\rangle}{\langle}\psi^{\sigma}_{m}|]\hat{\rho}_{\textrm{S}}\\},$

where the detailed procedures are exhibited in Appendix. The dressed state
dissipator is given by
$\mathcal{\hat{D}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|]\hat{\rho}_{\textrm{S}}=|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|\hat{\rho}_{\textrm{S}}|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma}_{m}|-\frac{1}{2}(|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|\hat{\rho}_{\textrm{S}}+\hat{\rho}_{\textrm{S}}|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|)$,
where the populations are naturally decoupled from the off-diagonal elements
of the density operator. The transition rates induced by the photon-bath
interaction in Eq. (2a) are

$\displaystyle\Gamma^{\textrm{R},+}_{m,\sigma}=$ $\displaystyle
m(f_{\sigma}-g_{\sigma})^{2}\gamma_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})n_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma}),~{}$
(10a) $\displaystyle\Gamma^{\textrm{R},-}_{m,\sigma}=$ $\displaystyle
m(f_{\sigma}-g_{\sigma})^{2}\gamma_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})[1+n_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})]~{}.$
(10b)

Physically,
$\Gamma^{\textrm{R},+}_{m,\sigma}~{}(\Gamma^{\textrm{R},-}_{m,\sigma})$
describes the excitation (relaxation) of the hybrid quantum system from the
eigenstate $|\psi^{\sigma}_{m-1(m)}{\rangle}$ to
$|\psi^{\sigma}_{m(m-1)}{\rangle}$ by absorbing (releasing) the energy
$\eta_{\sigma}\omega_{\sigma}$ from(to) the Rth bath. While the transition
rates from the qubit-bath interaction in Eq. (2b) are

$\displaystyle\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}=$
$\displaystyle\theta(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})G^{m^{\prime},\overline{\sigma}}_{m,\sigma}G^{m,\sigma}_{m^{\prime},\overline{\sigma}}\gamma_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})n_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})~{},$
(11a) $\displaystyle\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma}=$
$\displaystyle\theta(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})G^{m^{\prime},\overline{\sigma}}_{m,\sigma}G^{m,\sigma}_{m^{\prime},\overline{\sigma}}\gamma_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})[1+n_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})]~{},$
(11b)

where the positive energy gap is
$E^{m,\sigma}_{m^{\prime},\overline{\sigma}}=E_{m,{\sigma}}-E_{m^{\prime},{\overline{\sigma}}}$,
and the squeezing state overlap coefficient is defined as
$G^{n,\sigma}_{m,\bar{\sigma}}\equiv{\langle}\psi^{\sigma}_{n}|\hat{\sigma}_{x}|\psi^{\bar{\sigma}}_{m}{\rangle}=\;_{a}\\!{\langle}n|\exp\\{-(1/2l){(\alpha_{\sigma}-\alpha_{\overline{\sigma}})}(\hat{a}^{2}-\hat{a}^{{\dagger}^{2}})\\}|m{\rangle}_{a}$,
which is specified as pkral1990jmo

$\displaystyle~{}G^{m,\sigma}_{m^{\prime},\bar{\sigma}}$ $\displaystyle=$
$\displaystyle\frac{({v_{\sigma}}/{2u_{\sigma}})^{m/2}}{(m!m^{\prime}!u_{\sigma})^{1/2}}\sum^{\min\\{m,m^{\prime}\\}}_{l=0}\frac{m^{\prime}!H_{m^{\prime}-l}}{l!(m^{\prime}-l)!}\frac{m!H_{m-l}}{(m-l)!}$
(12)
$\displaystyle{\times}\left(\frac{2}{u_{\sigma}v_{\sigma}}\right)^{l/2}\left(\frac{-v^{*}_{\sigma}}{2u_{\sigma}}\right)^{(m^{\prime}-l)/2},$

with the coefficients $H_{m}=(-1)^{m/2}m!/(m/2)!$ for an even $m$ and
$H_{m}=0$ for an odd $m$,
$u_{\sigma}=\cosh(\alpha_{\sigma}-\alpha_{\overline{\sigma}})$, and
$v_{\sigma}=-\sinh(\alpha_{\sigma}-\alpha_{\overline{\sigma}})$. The rate
$\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}~{}(\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma})$
describes the process that the hybrid quantum system is assisted to transit
$|\psi^{\overline{\sigma}}_{m^{\prime}}{\rangle}~{}(|\psi^{{\sigma}}_{m}{\rangle})$
up (down) to
$|\psi^{{\sigma}}_{m}{\rangle}~{}(|\psi^{\overline{\sigma}}_{m^{\prime}}{\rangle})$
by exchanging energy $E^{m,\sigma}_{m^{\prime},\overline{\sigma}}$ with the
Qth bath. It need note that due to the specific selection rule of the
transition factor $G^{m,\sigma}_{m^{\prime},\bar{\sigma}}$ in Eq. (12),
$G^{m,\sigma}_{m^{\prime},\bar{\sigma}}{\neq}0$ tightly relies on the
condition $|m-m^{\prime}|$ to be even. Hence, the exchange between the qubit
and the Qth bath, quantified by
$\Gamma^{\textrm{Q},\pm}_{m,m^{\prime},\sigma}$, occurs only between the
states $|\psi^{\sigma}_{m}{\rangle}$ and $|\psi^{\sigma}_{m{\pm}2l}{\rangle}$.

###  II.3 Quantum master equation with FCS

We apply the FCS to individually count the energy current into the $u$th
bosonic thermal bath with the counting parameter $\chi_{\mu}$ after a long-
time evolution, which is based on the Markovian dressed master equation (14).
The technique of FCS was initially proposed by Levitov _et al._
levitov1992jetp ; levitov1996jmp  to investigate the electron flow and
fluctuations. Technically, based on the two-time nondemolition measurement
scheme, the generating function of the energy transfer is obtained as
mesp2009rmp

$\displaystyle~{}\mathcal{Z}(\mathrm{X},t)$ $\displaystyle{\equiv}$
$\displaystyle\textrm{Tr}\\{e^{i\sum_{\mu}\chi_{\mu}\hat{H}^{\mu}_{\textrm{B}}(0)}e^{-i\sum_{\mu}\chi_{\mu}\hat{H}^{\mu}_{\textrm{B}}(t)}\hat{\rho}_{\textrm{S}}(0)\\}$
(13) $\displaystyle=$
$\displaystyle\textrm{Tr}\\{\hat{\rho}_{\mathrm{X}}(t)\\},$

where the counting parameter set is $\mathrm{X}=\\{\chi_{\mu}\\}$, the
modified density operator
$\hat{\rho}_{\mathrm{X}}(t)=\hat{U}_{-\mathrm{X}}(t)\hat{\rho}_{\textrm{S}}(0)\hat{U}^{\dagger}_{\mathrm{X}}(t)$,
the unitary evolution operator is
$\hat{U}_{\mathrm{X}}(t)=\exp(-i\hat{H}_{\mathrm{X}}t)$, the initial density
operator of the hybrid system is $\hat{\rho}_{\textrm{S}}(0)$, and the
modified Hamiltonian is
$\hat{H}_{\mathrm{X}}{\equiv}\exp(i\sum_{\mu}\chi_{\mu}\hat{H}^{\mu}_{\textrm{B}}/2)\hat{H}\exp(-i\sum_{\mu}\chi_{\mu}\hat{H}^{\mu}_{\textrm{B}}/2)$,
which is specified as
$\hat{H}_{\mathrm{X}}=\hat{H}_{\textrm{S}}+\sum_{\mu}[\hat{H}^{\mu}_{\textrm{B}}+\hat{V}_{\mu}(\chi_{\mu})]$,
with
$\hat{V}_{\textrm{R}}(\chi_{\textrm{R}})=(\hat{a}^{\dagger}+\hat{a})\sum_{k}(g_{k,\textrm{R}}e^{i\omega_{k}\chi_{\textrm{R}}/2}\hat{b}^{\dagger}_{k,\textrm{R}}+e^{-i\omega_{k}\chi_{\textrm{R}}/2}g^{*}_{k,\textrm{R}}\hat{b}_{k,\textrm{R}})$
and
$\hat{V}_{\textrm{Q}}(\chi_{\textrm{Q}})=\hat{\sigma}_{x}\sum_{k}(g_{k,\textrm{Q}}e^{i\omega_{k}\chi_{\textrm{Q}}/2}\hat{b}^{\dagger}_{k,\textrm{Q}}+g^{*}_{k,\textrm{Q}}e^{-i\omega_{k}\chi_{\textrm{Q}}/2}\hat{b}_{k,\textrm{Q}})$.
Then, employing the Born-Markov approximation, the density operator is
decomposed as
$\hat{\rho}_{\mathrm{X}}(t)=\hat{\rho}^{\textrm{S}}_{\mathrm{X}}(t){\otimes}\hat{\rho}_{\textrm{B}}$.
Within the framework of dressed master equation, we perturb
$\hat{V}_{\mu}(\chi_{\mu})$ up to the second order and obtain

$\displaystyle~{}\frac{d}{dt}\hat{\rho}^{\textrm{S}}_{\textrm{X}}$
$\displaystyle=$
$\displaystyle-i[\hat{H}_{\textrm{S}},\hat{\rho}^{\textrm{S}}_{\textrm{X}}]$
(14)
$\displaystyle+\sum_{m,\sigma}\\{\Gamma^{\textrm{R},+}_{m,\sigma}\mathcal{\hat{D}}_{\chi_{\textrm{R}}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma}_{m-1}|]\hat{\rho}^{\textrm{S}}_{\textrm{X}}$
$\displaystyle+\Gamma^{\textrm{R},-}_{m,\sigma}\mathcal{\hat{D}}_{-\chi_{\textrm{R}}}[|\psi^{\sigma}_{m-1}{\rangle}{\langle}\psi^{\sigma}_{m}|]\hat{\rho}^{\textrm{S}}_{\textrm{X}}\\}$
$\displaystyle+\sum_{m,m^{\prime},\sigma}\\{\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}\mathcal{\hat{D}}_{\chi_{\textrm{Q}}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\overline{\sigma}}_{m^{\prime}}|]\hat{\rho}^{\textrm{S}}_{\textrm{X}}$
$\displaystyle+\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma}\mathcal{\hat{D}}_{-\chi_{\textrm{Q}}}[|\psi^{\overline{\sigma}}_{m^{\prime}}{\rangle}{\langle}\psi^{\sigma}_{m}|]\hat{\rho}^{\textrm{S}}_{\textrm{X}}\\},$

where the modified dressed dissipator is given by
$\mathcal{\hat{D}}_{\chi}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|]\hat{\rho}^{\textrm{S}}_{\textrm{X}}=e^{-i{\chi}E^{m,\sigma}_{l,\sigma^{\prime}}}|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|\hat{\rho}^{\textrm{S}}_{\textrm{X}}|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma}_{m}|-\frac{1}{2}(|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|\hat{\rho}^{\textrm{S}}_{\textrm{X}}+\hat{\rho}^{\textrm{S}}_{\textrm{X}}|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|)$.
Hence, we obtain the cumulant generating function at steady state as

$\displaystyle\mathcal{G}(\textrm{X})=\lim_{t{\rightarrow}\infty}\frac{1}{t}\ln[\mathcal{Z}(\textrm{X},t)].$
(15)

Then, the $n$th cumulant into the $u$th thermal bath is given by
$J^{(n)}_{u}=[{\partial}^{n}\mathcal{G}(\textrm{X})/{\partial}(i{\chi_{u}})^{n}]{\Big{|}}_{\textrm{X}=0}$.
In the following, we focus on the current fluctuations into the Qth thermal
bath. Accordingly, the steady-state heat current, noise power, and skewness
are given by

$\displaystyle J_{ss}=$
$\displaystyle[{\partial}\mathcal{G}(\textrm{X})/{\partial}(i{\chi_{\textrm{Q}}})]{\Big{|}}_{\textrm{X}=0}~{},$
(16a) $\displaystyle J^{(2)}_{ss}=$
$\displaystyle[{\partial}^{2}\mathcal{G}(\textrm{X})/{\partial}(i{\chi_{\textrm{Q}}})^{2}]{\Big{|}}_{\textrm{X}=0}~{},$
(16b) $\displaystyle J^{(3)}_{ss}=$
$\displaystyle[{\partial}^{3}\mathcal{G}(\textrm{X})/{\partial}(i{\chi_{\textrm{Q}}})^{3}]{\Big{|}}_{\textrm{X}=0}~{}.$
(16c)

In particular, by analyzing the dynamical transition processes in Eq. (14),
the current can also be alternatively expressed as

$\displaystyle~{}J_{ss}=\sum_{m,m^{\prime},\sigma}[\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma}P_{m,{\sigma}}-\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}P_{m^{\prime},\overline{\sigma}}],$
(17)

where $P_{m,{\sigma}}$ is the steady-state population of the eigenstate
$|\psi^{\sigma}_{m}{\rangle}$.

Figure 2: (Color online) (a) Influence of the temperature bias ${\Delta}T$ on
the rescaled steady-state heat current $J_{ss}/\lambda^{2}$ from weak to
strong qubit-photon hybridizations, e.g., $\lambda=0.001,0.01$, and $0.1$; (b)
bird view of the heat current $J_{ss}$ by both modulating $\Delta{T}$ and
$\lambda$; (c) and (d) are cyclic heat transitions to depict the first term of
$I_{m,1}$ and $I_{m,0}$, respectively; (e) effect of the temperature bias on
heat current component $I_{m,\sigma}$ given by Eqs. (20a-20b) in the weak
hybridization regime to describe the origin of the NDTC; (f) two different
types of heat transfer processes between eigenstates at strong qubit-photon
hybridization, i.e.,
$|\psi^{1}_{m+2k}{\rangle}{\rightarrow}|\psi^{0}_{m}{\rangle}$ and
$|\psi^{0}_{m+2l}{\rangle}{\rightarrow}|\psi^{1}_{m}{\rangle}$
($m,k,l{\in}\mathrm{N}$), where the thickness of the blue arrowed lines
denoting the transition rate strengthes. The other parameters are given by
$T_{\mathrm{R}}=T_{0}+{\Delta}T/2$, $T_{\mathrm{Q}}=T_{0}-{\Delta}T/2$,
$T_{0}=1$, $\omega_{a}=1$, $\varepsilon=1$,
$\alpha_{\mathrm{R}}=\alpha_{\mathrm{Q}}=10^{-3}$, and $\omega_{c}=10$.

##  III Nonequilibrium thermal transport

The NDTC is considered as one representative kind of nonlinear effects in
nonequilibrium thermal transport bli2006apl . It was initially proposed by Li
_et al._ nbli2012rmp  in nonlinear phononic lattices, which is characterized
as the crucial component of the thermal transistor to efficiently manage heat
transport. The NDTC describes the effect that the heat current is anomalously
suppressed by increasing the bath temperature bias within two terminal setup.
In recent years, the NDTC has been extensively investigated in small quantum
systems, e.g., hybrid quantum system cwang2021cpl ; jren2013prb ;
afornieri2016prb  and SBM dsegal2006prb . However, the influence of the
quadratic qubit-photon interaction on the NDTC is still illusive.

To investigate the NDTC in our system, we first show the steady-state heat
current $J_{ss}/\lambda^{2}$ (17) by modulating the temperature bias
${\Delta}T$ ($T_{\textrm{R}}=T_{0}+{\Delta}T/2$ and
$T_{\textrm{Q}}=T_{0}-{\Delta}T/2$) in Fig. 2(a). In the negative temperature
bias regime, the heat current shows trivial monotonic enhancement by
increasing $|{\Delta}T|$ from weak to strong qubit-photon hybridization
strengthes. In sharp contrast, it is interesting to find that for positive
${\Delta}T$, $J_{ss}$ always exhibits nontrivial nonmonotic behavior by
enlarging ${\Delta}T$, i.e., $J_{ss}$ enhanced with small ${\Delta}T$ and
suppressed at large ${\Delta}T$. This clearly demonstrates the emergence of
the NDTC. Moreover, such novel behavior of the NDTC is bounded by
$\lambda<0.25$. Hence, we focus on the positive temperature bias in the
following. In Fig. 2(b), we also show a comprehensive view of the cooperative
effect of both the hybridization strength and the temperature bias. The heat
current forms a global peak at $\lambda{\approx}0.15$ and
${\Delta}T{\approx}0.6$. Though not shown in the paper, it is confirmed that
the NDTC persists at finite energy bias regime, i.e.,
$\varepsilon{\neq}\omega_{a}$. It should be noted that the robustness of the
NDTC in the current model over the wide qubit-photon hybridization regime is
nontrivial, which is different from previous results in the limiting
conditions, e.g., linearly qubit-phonon hybrid system at weak hybridization
cwang2021cpl , nonequilibrium SBM at strong system-bath dissipation
dsegal2006prb , and electron-magnon scattering interface with weak tunneling
jren2013prb . Below, we try to explore the underlying mechanism for the
appearance of the NDTC in representative qubit-photon hybridization regimes.

We first analyze the expression of heat current in the weak qubit-photon
hybridization case. From dynamical equation in Eq. (14) and the expression of
heat current in Eq. (17), it is known that the transition rates are
fundamental to establish the dynamical energy exchange processes, which leads
to the steady-state populations and heat current. In particular, the
renormalization factor is reduced to
$\eta_{\sigma}{\approx}1-2\lambda^{2}/\omega^{2}_{a}$, the coefficients become
$f_{\sigma}{\approx}1+\lambda^{2}/\omega^{2}_{a}$,
$g_{\sigma}{\approx}(-1)^{\sigma+1}\lambda/\omega_{a}$, and
$\alpha_{\sigma}{\approx}(-1)^{\sigma+1}\lambda/\omega_{a}$. Thus, the
squeezing state overlap coefficient defined in Eq. (12) is simplified as
$G^{m,\sigma}_{m^{\prime},\overline{\sigma}}{\approx}\delta_{m,m^{\prime}}+(-1)^{\sigma}({\lambda}/{\omega_{a}})\sqrt{m^{\prime}(m^{\prime}-1)}\delta_{m,m^{\prime}-2}-(-1)^{\sigma}({\lambda}/{\omega_{a}})\sqrt{m(m-1)}\delta_{m,m^{\prime}+2}$.
Consequently, the transition rates associated with the $\mathrm{Q}$th
reservoir is given by

$\displaystyle\Gamma^{\mathrm{Q},\pm}_{m,m^{\prime},\sigma}$
$\displaystyle{\approx}$
$\displaystyle\Gamma^{\mathrm{Q},\pm}_{m,m,1}\delta_{m,m^{\prime}}\delta_{\sigma,1}+\Gamma^{\mathrm{Q},\pm}_{m,m-2,1}\delta_{m,m^{\prime}+2}\delta_{\sigma,1}$
(18)
$\displaystyle+\Gamma^{\mathrm{Q},\pm}_{m,m+2,0}\delta_{m,m^{\prime}-2}\delta_{\sigma,0}.$

The first rate component
$\Gamma^{Q,\pm}_{m,m,1}{\approx}\gamma_{Q}(\pm\varepsilon)n_{Q}(\pm\varepsilon)$
describes the individual process that the qubit flips up (down) assisted by
the Qth bath without changing the photon state. While the second and third
components
$\Gamma^{Q,\pm}_{m,m-2,1}{\approx}m(m-1)({\lambda}/{\omega_{a}})^{2}\gamma_{Q}({\pm}2\omega_{a}\pm\varepsilon)n_{Q}(\pm
2\omega_{a}\pm\varepsilon)$ and
$\Gamma^{Q,\pm}_{m-2,m,0}{\approx}m(m-1)({\lambda}/{\omega_{a}})^{2}\gamma_{Q}(\pm
2\omega_{a}\mp\varepsilon)n_{Q}(\pm 2\omega_{a}\mp\varepsilon)$ explicitly
capture the contribution of cooperative qubit-photon scattering process, i.e.,
the qubit flip is simultaneously accompanied by two photons emission
(absorption), under the influence of both the Rth and Qth thermal baths.
Hence, $\Gamma^{Q,\pm}_{m,m-2,1}$ and $\Gamma^{Q,\pm}_{m-2,m,0}$ are crucial
to analytically obtain the expression of heat current. Consequently, the
current defined in Eq. (17) with the leading order of
$(\lambda/\omega_{a})^{2}$ is expressed as

$\displaystyle~{}J_{ss}{\approx}\frac{2\lambda^{2}}{\omega_{a}}\sum^{\infty}_{m=2}m(m-1)(I_{m,1}+I_{m,0}),$
(19)

where the current components are

$\displaystyle I_{m,1}=$
$\displaystyle\frac{\gamma_{Q}(2\omega_{a}+\varepsilon)e^{-\beta_{\mathrm{R}}m\omega_{a}}}{[1+2n_{\mathrm{Q}}(\varepsilon)][1+n_{\mathrm{R}}(\omega_{a})]n_{\mathrm{R}}(2\omega_{a})}~{}$
(20a)
$\displaystyle{\times}\\{[1+n_{\mathrm{Q}}(2\omega_{a}+\varepsilon)]n_{\mathrm{Q}}(\varepsilon)n_{\mathrm{R}}(2\omega_{a})$
$\displaystyle-
n_{\mathrm{Q}}(2\omega_{a}+\varepsilon)[1+n_{\mathrm{Q}}(\varepsilon)][1+n_{\mathrm{R}}(2\omega_{a})]\\},$
$\displaystyle I_{m,0}=$
$\displaystyle\frac{\gamma_{Q}(2\omega_{a}-\varepsilon)e^{-\beta_{\mathrm{R}}m\omega_{a}}}{[1+2n_{\mathrm{Q}}(\varepsilon)][1+n_{\mathrm{R}}(\omega_{a})]n_{\mathrm{R}}(2\omega_{a})}~{}$
(20b)
$\displaystyle{\times}\\{[1+n_{\mathrm{Q}}(2\omega_{a}-\varepsilon)][1+n_{\mathrm{Q}}(\varepsilon)]n_{\mathrm{R}}(2\omega_{a})$
$\displaystyle-
n_{\mathrm{Q}}(2\omega_{a}-\varepsilon)n_{\mathrm{Q}}(\varepsilon)[1+n_{\mathrm{R}}(2\omega_{a})]\\}.$

Specifically, the component $I_{m,1}$ is contributed by two competing cyclic
transitions, where the first term characterized as
$[1+n_{\mathrm{Q}}(2\omega_{a}+\varepsilon)]n_{\mathrm{Q}}(\varepsilon)n_{\mathrm{R}}(2\omega_{a})$
is specified by the transition paths

$\displaystyle|\psi^{1}_{m}{\rangle}{\rightarrow}|\psi^{0}_{m-2}{\rangle}{\rightarrow}|\psi^{1}_{m-2}{\rangle}{\rightarrow}|\psi^{1}_{m}{\rangle},$
(21)
$\displaystyle|\psi^{1}_{m}{\rangle}{\rightarrow}|\psi^{0}_{m-2}{\rangle}{\rightarrow}|\psi^{0}_{m}{\rangle}{\rightarrow}|\psi^{1}_{m}{\rangle}.$

And the microscopic physical pictures are also graphed in Fig. 2(c).
Accordingly, the second term can be straightforwardly obtained by reversing
the transition directions of Fig. 2(c). In analogy, $I_{m,0}$ is also composed
by opposite cyclic components, where the first term proportional to
$[1+n_{\mathrm{Q}}(2\omega_{a}-\varepsilon)][1+n_{\mathrm{Q}}(\varepsilon)]n_{\mathrm{R}}(2\omega_{a})$
is contributed by transition paths

$\displaystyle|\psi^{0}_{m}{\rangle}{\rightarrow}|\psi^{1}_{m-2}{\rangle}{\rightarrow}|\psi^{0}_{m-2}{\rangle}{\rightarrow}|\psi^{0}_{m}{\rangle},$
(22)
$\displaystyle|\psi^{0}_{m}{\rangle}{\rightarrow}|\psi^{1}_{m-2}{\rangle}{\rightarrow}|\psi^{1}_{m}{\rangle}{\rightarrow}|\psi^{0}_{m}{\rangle},$

and physically depicted in Fig. 2(d).

Then, we analyze the mechanism for emergence of the NDTC with weak qubit-
resonator hybridization. At weak temperature bias limit, the heat current Eq.
(19) behaves $J_{ss}{\propto}{\Delta}T$, which denotes that $J_{ss}$ shows
linear enhancement by increasing ${\Delta}T$. While in the large ${\Delta}T$
regime, e.g., $T_{\textrm{R}}{\approx}2T_{0}$ and $T_{\textrm{Q}}{\approx}0$,
the directional transitions, i.e.,
$|\psi^{0}_{m-2}{\rangle}{\rightarrow}|\psi^{1}_{m-2}{\rangle}$,
$|\psi^{0}_{m-2}{\rangle}{\rightarrow}|\psi^{0}_{m}{\rangle}$ and
$|\psi^{0}_{m-2}{\rangle}{\rightarrow}|\psi^{1}_{m}{\rangle}$, are dramatic
suppressed. Then, all cyclic transitions to compose $I_{m,1}$ break down,
which naturally leads to $I_{m,1}{\approx}0$. Hence, $I_{m,1}$ clearly
exhibits the nonmonotonic behavior by increasing the bias ${\Delta}T$.
Moreover, the magnitude of $I_{m,1}$ is larger than $I_{m,0}$ over a wide
temperature bias regime ($0<{\Delta}T{<}1.3$), with $m=2,3$ exemplified in
Fig. 2(e). Therefore, we conclude that the effect of the NDTC originates from
the suppression of cyclic components of $I_{m,1}$. However, the finite residue
of the heat current at large temperature bias limit results from the
robustness of $I_{m,0}$, of which all cyclic components show monotonic
enhancement by increasing ${\Delta}T$.

Next, we turn to investigate the mechanism of the NDTC beyond weak qubit-
photon hybridization. Due to multi-photon transfer processes, it is rather
difficult to obtain the analytical expression of heat current. However, we may
resort to transition rates $\Gamma^{\pm}_{m,m^{\prime},\sigma}$ for the
physical insight of energy exchange processes. In particular, in large
hybridization regime the coefficient
$G^{m,\sigma}_{m^{\prime},\overline{\sigma}}$ Eq. (12) will be strongly
suppressed with the increase of number difference $|m-m^{\prime}|$. Combining
the dramatic frequency mismatch of bosonic squeezing modes Eq. (5), i.e.,
$\eta_{1}\omega_{1}{\gg}\eta_{0}\omega_{0}$, the directional transition from
$|\psi^{0}_{m+2l}{\rangle}$ and $|\psi^{1}_{m}{\rangle}$ is almost inhibited
between squeezing states $|\psi^{0}_{m}{\rangle}$ and
$|\psi^{1}_{m^{\prime}}{\rangle}$ [shown in Fig. 2(f)], which is restricted by
$l{\gg}1$ to overcome the positive energy bias $E^{0}_{m+2l}{>}E^{1}_{m}$. On
the contrary, the transition from $|\psi^{1}_{m+2k}{\rangle}$ and
$|\psi^{0}_{m}{\rangle}$ is generally allowed under the condition
$E^{1}_{m+2k}{>}E^{0}_{m}$. Hence, the transition rate
$\Gamma^{\textrm{Q},\pm}_{m+2l,m,0}$ becomes much smaller than
$\Gamma^{\textrm{Q},\pm}_{m+2k,m,1}$ with $l{\gg}k$ over a wide temperature
bias regime, as schemed in Fig. 2(f). Consequently, at large temperature limit
(e.g., ${\Delta}T{\approx}2T_{0}$) the qubit is almost localized to the ground
state $|0{\rangle}$. Such localization of the qubit apparently blocks the
energy exchange between the qubit and the Qth thermal bath, which finally
results in the emergence of the NDTC.

Figure 3: (Color online) Thermal rectification factor $\mathcal{R}$ defined in
Eq. (23) versus the temperature bias ${\Delta}T$ and the qubit-photon
hybridization strength $\lambda$. The other parameters are the same as those
in Fig. 2.

The novel thermal transport effect in this system can also be characterized by
the thermal rectification effect, which in analogy with the electronic diode,
was theoretically proposed by Li _et al._ bli2004prl  to study the heat flux
asymmetry in nonlinear phononic lattices. In recent years, the concept of
thermal rectification has been extensively investigated in various hybrid
quantum setups cwang2021cpl ; dsegal2005prl ; mjmperez2015nn . The thermal
rectification is applied to characterize nonreciprocity of the nonequilibrium
system, that within the two terminal setup heat current is larger in one
direction than the opposite counterpart by only exchanging two bath
temperatures. Quantitatively, the thermal rectification is defined by the
factor dsegal2005prl ; lfzhang2009prb

$\displaystyle~{}\mathcal{R}=\frac{|J_{ss}({\Delta}T)+J_{ss}(-{\Delta}T)|}{\max\\{|J_{ss}({\Delta}T)|,|J_{ss}(-{\Delta}T)|\\}},$
(23)

where $J_{ss}(\pm{\Delta}T)$ denote the current under the conditions
$T_{\mathrm{R}}=T_{0}\pm\Delta{T}/2$ and $T_{\mathrm{Q}}=T_{0}\mp\Delta{T}/2$.
From the definition in Eq. (23), it is known that the rectification factor
shows symmetric behavior as $\mathcal{R}(\Delta{T})=\mathcal{R}(-\Delta{T})$.
Moreover, the thermal rectification becomes significant in the upper limit
$\mathcal{R}=1$, where the extreme nonreciprocity occurs as
$|J_{ss}({\Delta}T)|{\gg}|J_{ss}(-{\Delta}T)|$. While the rectification
disappears as $\mathcal{R}=0$, which denotes the complete reciprocity
behavior, i.e., $J_{ss}(-{\Delta}T){=}-J_{ss}({\Delta}T)$.

We plot Fig. 3 to investigate the behavior of thermal rectification by tuning
the temperature bias ${\Delta}T$. In small temperature bias regime, the
rectification factor $\mathcal{R}$ generally keeps finite from weak to strong
qubit-photon hybridizations. It results from the intrinsic nonreciprocity of
the qubit-resonator hybrid system, owning the spin and boson nature
respectively. Quantitatively, we may obtain some insight from Eq. (19) at weak
qubit-resonator hybridization. By exchanging thermal bath temperatures, the
cyclic transition weight of current components are significantly changed. For
instance, the exchanged cyclic kernel in $I_{m,1}$ is given by
$\\{[1+n_{\mathrm{R}}(2\omega_{a}+\varepsilon)]n_{\mathrm{R}}(\varepsilon)n_{\mathrm{Q}}(2\omega_{a})-n_{\mathrm{R}}(2\omega_{a}+\varepsilon)[1+n_{\mathrm{R}}(\varepsilon)][1+n_{\mathrm{Q}}(2\omega_{a})]\\}$,
which is distinct from the counterpart
$\\{[1+n_{\mathrm{Q}}(2\omega_{a}+\varepsilon)]n_{\mathrm{Q}}(\varepsilon)n_{\mathrm{R}}(2\omega_{a})-n_{\mathrm{Q}}(2\omega_{a}+\varepsilon)[1+n_{\mathrm{Q}}(\varepsilon)][1+n_{\mathrm{R}}(2\omega_{a})]\\}$
in Eq. (20a). This may partially explains the existence of finite thermal
rectification. While in large temperature bias regime, the rectification
factor is further enhanced. In particular at strong hybridization (e.g.,
$\lambda{=}0.1$), the factor approaches the unit
($\mathcal{R}{\rightarrow}1$), and keeps robust. Such enhancement results from
emergence of the NDTC only at positive ${\Delta}T$, i.e.,
$-J_{ss}(\Delta{T}=-2){\gg}J_{ss}(\Delta{T}=2)$, as illustrated in Fig. 2(a).
Therefore, the qubit-resonator hybrid system with quadratic interaction may
provide a potential platform to realize the efficient thermal diode.

Figure 4: (Color online) (a) Noise power $J^{(2)}_{ss}$ and (b) skewness
$J^{(3)}_{ss}$ defined in Eqs. (16b-16c) by tuning the bath temperature bias
${\Delta}T$ and qubit-photon hybridization strength $\lambda$. The other
parameters are the same as those given in Fig. 2.

The higher-order thermal transport effect can be characterized by
investigating the steady-state current fluctuations, which are typically
quantified by the noise power and skewness. The noise power $J^{(2)}_{ss}$ in
Eq. (16b) characterizes zero-frequency power spectrum of the heat noise
generated by the stochastic heat exchange between the hybrid system and
$\mathrm{Q}$th bath, which is generally not directly related with the heat
current at small quantum systems levitov1992jetp ; mesp2009rmp ; fzhan2011prb
. Thus, the noise power is considered as one indispensable factor to quantify
the current fluctuation. It is interesting to find that $J^{(2)}_{ss}$ shows
global peak with small temperature bias at strong hybridization
($\lambda{\approx}0.18$), whereas it is generally suppressed as the
temperature bias increases, as shown in Fig. 4(a). Hence, the monotonic
suppression behavior of the noise power by increasing ${\Delta}T$ is distinct
from the steady-state heat current. While the finite skewness $J^{(3)}_{ss}$
defined in Eq. (16c) characterizes asymmetric probability distribution of the
transferred heat energy around the mean value, whereas the vanishing skewness
corresponds to a symmetric Gaussian distribution rbelousov2016pre . The
dramatic positive skewness can be found at strong qubit-photon hybridization,
which exhibits that the probability distribution of the transferred heat
deviates from the Gaussian case. Furthermore, it is significantly enhanced at
large bias limit (${\Delta}T=2$) , as shown in Fig. 4(b). Therefore, the
current fluctuations exhibit nontrivial features compared to the heat current,
which fertilize heat transport pictures of the hybrid quantum system.

Via the Fourier transformation the generating function Eq. (13) to count heat
flow into the $\mu$th thermal bath, can be re-expressed as nasinitsyn2007epl ;
jren2010prl
$\mathcal{Z}_{t}(\chi_{\mu})=\int^{\infty}_{-\infty}d{Q}P_{t}(Q)e^{iQ\chi_{\mu}}$,
where $P_{t}(Q)$ denotes the conditional probability to transport heat energy
$Q$ into the $\mu$th bath by the time $t$. Then, the cumulant generating
function at steady state is obtained as
$\mathcal{G}(\chi_{\mu})=\lim_{t{\rightarrow}\infty}\frac{1}{t}\ln[\mathcal{Z}_{t}(\chi_{\mu})]$.
Consequently, the noise power and skewness are given by
$J^{(2)}_{ss}=\lim_{t{\rightarrow}\infty}\frac{1}{t}[{\langle}(Q-{\langle}Q{\rangle})^{2}{\rangle}]$
and
$J^{(3)}_{ss}=\lim_{t{\rightarrow}\infty}\frac{1}{t}[{\langle}(Q-{\langle}Q{\rangle})^{3}{\rangle}]$,
where
${\langle}Q^{n}{\rangle}=\int^{\infty}_{-\infty}d{Q}Q^{n}P_{t}(Q)/\mathcal{Z}_{t}(0)$.
From the practical view, $Q=\int^{t}_{0}{J(\tau)d\tau}$ could be
experimentally detected via the trajectory measurement yang2020nc , with
$J(\tau)$ the instantaneous heat current. By collecting many trajectories of
quantum heat transport after long-time evolution, the distribution
$P_{t{\rightarrow}\infty}(Q)$ can be approximately obtained. Under such
measurement scheme, noise power and skewness could be experimentally measured.

##  IV Stationary Photon quadrature squeezing

Figure 5: (Color online) Photon quadrature squeezing factor
$\xi^{2}_{\textrm{B}}$ defined in Eq. (24) (a) at equilibrium
$T_{\mathrm{R}}=T_{\mathrm{Q}}=T$, with the inset showing the angle
distribution of $({\Delta}\hat{X}_{\theta})^{2}$ at $\lambda=0.2$, and (b)
far-from-equilibrium by tuning two bath temperatures $T_{\mathrm{R}}$ and
$T_{\mathrm{Q}}$ individually, with $\lambda=0.2$. The other parameters are
the same as those given in Fig. 2.

Quadrature squeezing is considered as one important nonclassical effect, where
the noise of a single bosonic quadrature breaks the Heisenberg fluctuation
limit at vacuum dstoler1970prd ; dstoler1971prd ; ula2016ps . Typically, from
the dynamical view the noise is suppressed at some momentum accompanied by the
enhancement of noise at other times jrk1988pra ; ax2020np . While from the
stationary view (i.e., ground state or steady state), the noise of one bosonic
quadrature component is below the zero-point level, whereas the other
quadrature one surpasses the value in the vacuum state as2013prl ; ak2013pra ;
eew2015science . Generally, the quadrature squeezing factor can be defined as
jma2011pr

$\displaystyle~{}\xi^{2}_{\textrm{B}}=\min_{\theta{\in}[0,2\pi)}({\Delta}\hat{X}_{\theta})^{2}$
(24)

to quantify the photon squeezing at steady state, where
$\hat{X}_{\theta}=\hat{X}\cos\theta+\hat{P}\sin\theta$, with the canonical
position operator $\hat{X}=\hat{a}^{\dagger}+\hat{a}$ and the canonical
momentum operator $\hat{P}=i(\hat{a}^{\dagger}-\hat{a})$. If
$\xi^{2}_{\textrm{B}}<1$, the photon is squeezed, whereas the squeezing effect
vanishes as $\xi^{2}_{\textrm{B}}>1$. While for the photon mode at the
coherent state, the factor becomes the unit. The concept of quadrature
squeezing has been extensively applied in various hybrid quantum setups, e.g.,
optomechanics ak2013pra ; eew2015science ; jqliao2011pra ; xylu2015pra  and
quantum electrodynamics systems (QEDs) mm2008prl ; jkxie2020pra . Here, we
investigate the steady-state photon quadrature squeezing under the influence
of two bath temperatures and the qubit-photon hybridization strength.

We first study the photon squeezing at equilibrium in Fig. 5(a), i.e.,
$T_{\mathrm{R}}=T_{\mathrm{Q}}=T$. It is found that the squeezing effect
($\xi^{2}_{\textrm{B}}<1$) occurs in the low-temperature regime, and is
gradually enhanced as the qubit-photon hybridization strength increases. By
analyzing the angle distribution of $\xi^{2}_{\textrm{B}}<1$ in the squeezing
regime [e.g., $\lambda=0.2$ and $T=0$ in the inset of Fig. 5(a)], we find that
the noise of one quadrature component $\hat{P}$ is significantly reduced,
whereas the other component $\hat{X}$ fails to squeeze. Moreover, we
investigate the effect of bath temperatures bias on the squeezing factor at
strong hybridization ($\lambda=0.2$) in Fig. 5(b). The squeezing factor is
clearly shown asymmetrically driven by two bath temperatures. Particularly, in
the low temperature of $T_{\mathrm{R}}$, e.g., $T_{\mathrm{R}}=0$, the
squeezing effect persists as $T_{\mathrm{Q}}<0.9$. Therefore, we conclude that
the photon squeezing effect favors strong qubit-photon hybridization, and
persists even at comparatively high temperature of $T_{\mathrm{Q}}$.

##  V Conclusion

In conclusion, we studied the steady-state heat flow and photon quadrature
squeezing in an open qubit-resonator quantum system, where the optical
resonator is quadratically coupled to the qubit, with each individually
interacting with the bosonic thermal bath. We applied the dressed master
equation combined with the full counting statistics to quantify the heat
current and current fluctuations, which is applicable after long time
evolution with negligible off-diagonal elements of the system density matrix
in the squeezed states basis $\\{|\psi^{\sigma}_{m}{\rangle}\\}$. The
transition rates $\Gamma^{\mathrm{R},\pm}_{m,\sigma}$ and
$\Gamma^{\mathrm{Q},\pm}_{m,m^{\prime},\sigma}$, characterizing the energy
exchange between the hybrid quantum system components (optical resonator and
qubit) and the corresponding thermal baths ($\mathrm{R}$th and
$\mathrm{Q}$th), are expressed at Eqs. (10a-10b) and Eqs. (11a-11b),
respectively. In particular, the rates
$\Gamma^{\mathrm{Q},\pm}_{m,m^{\prime},\sigma}$ are tightly related with the
squeezed state overlap coefficient
$G^{m,\sigma}_{m^{\prime},\overline{\sigma}}$ in Eq. (12), which clarifies the
specific transition rule between squeezed states with opposite spin branches.

We investigated the steady-state heat current by both tuning the bath
temperature bias and the qubit-resonator hybridization strength. It is
interesting to find that the heat current always shows nonmonotonic behavior
by increasing the temperature bias, i.e., the NDTC, from weak to strong
hybridizations. We also exploited the underlying mechanism to exhibit the
NDTC. At weak hybridization, the heat current is analytically contributed by
two components $I_{m,\sigma}$ in Eqs. (20a-20b), which are both described by
the cyclic heat exchange processes. And the suppression of cyclic heat
transitions of $I_{m,1}$ dominates the emergence of the NDTC. While at strong
hybridization, the almost directional transitions from the squeezed states
branch $\\{|\psi^{1}_{m}{\rangle}\\}$ to $\\{|\psi^{0}_{m}{\rangle}\\}$ mainly
contribute to the heat transfer at finite temperature bias, due to the large
frequency mismatch of two squeezed photon modes in Eq. (5). Then, we
investigated the thermal rectification effect in Eq. (23), which becomes
highly efficient at strong qubit-resonator hybridization and large temperature
bias. The intrinsic nonreciprocity of the qubit-resonator hybrid system and
the emergence of the NDTC contribute to such highly thermal rectification.
Moreover, we analyzed steady-state noise power and skewness, which show global
peaks with large hybridization strength at small and large temperature bias
limits, respectively. These novel fluctuation features fertilize heat
transport properties of the hybrid system.

We also investigated the photon quadrature squeezing effect quantified by the
generalized factor $\xi^{2}_{\mathrm{B}}$ defined in Eq. (24). At equilibrium
($T_{\mathrm{R}}=T_{\mathrm{Q}}=T$), the squeezing behavior occurs in the
comparative low-temperature regime, and becomes enhanced with the increase of
the qubit-resonator hybridization strength. The variance of the momentum
operator $\hat{P}$, corresponding to $\theta=\pi/2$, is dramatically
suppressed, whereas the other quadrature component $\hat{X}$($\theta=0$)
becomes noisy. While out of equilibrium ($T_{\mathrm{R}}{\neq}T_{\mathrm{Q}}$)
and strong hybridization, the squeezing factor shows asymmetric behavior by
tuning two bath temperatures.

We hope that the quadratically coupled qubit-resonator hybrid quantum system
can be considered as one potential hybrid platform to realize efficient
quantum thermal devices, e.g., the NDTC and thermal rectifier, and
nonclassical behavior of photons, e.g., quadrature squeezing.

## ACKNOWLEDGEMENTS

C.W. is 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. H.C. acknowledges
the National Natural Science Foundation of China under Grant No. 11704338.
J.-Q.L. is supported in part by National Natural Science Foundation of China
(Grants No. 11822501, No. 11774087, and No. 11935006), Hunan Science and
Technology Plan Project (Grant No. 2017XK2018), and the Science and Technology
Innovation Program of Hunan Province (Grant No. 2020RC4047).

## Appendix: Derivation of the quantum dressed master equation

Based on the Born approximation, the density operator of the total system is
decomposed by
$\hat{\rho}_{\textrm{tot}}{\approx}\hat{\rho}_{\textrm{S}}{\otimes}\hat{\rho}_{\textrm{B}}$,
where $\hat{\rho}_{\textrm{S}}$ is the density operator of the qubit-resonator
hybrid system, and
$\hat{\rho}_{\textrm{B}}=\Pi_{u=\textrm{R},\textrm{Q}}\exp[-\hat{H}^{u}_{\textrm{B}}/(k_{B}T_{u})]/Z$
is the density operator of thermal baths at equilibrium, with
$Z=\textrm{Tr}_{\textrm{B}}\\{\Pi_{u=\textrm{R},\textrm{Q}}\exp[-\hat{H}^{u}_{\textrm{B}}/(k_{B}T_{u})]\\}$
the partition function. Then, we perturb system-bath interactions, i.e.,
$\hat{V}_{\textrm{R}}$ and $\hat{V}_{\textrm{Q}}$, to obtain the quantum
master equation as

$\displaystyle\frac{d}{dt}\hat{\rho}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle i[\hat{\rho}_{\textrm{S}},\hat{H}_{\textrm{S}}]$
$\displaystyle+\sum_{\mu=\textrm{R},\textrm{Q}}\int^{\infty}_{0}{d\tau}\textrm{Tr}_{\textrm{B}}\\{[\hat{V}_{\mu},[\hat{V}_{\mu}(-\tau),\hat{\rho}_{\textrm{S}}{\otimes}\hat{\rho}_{\textrm{B}}]]\\},$

where the commutation relation is
$[\hat{X},\hat{Y}]=\hat{X}\hat{Y}-\hat{Y}\hat{X}$, and
$\textrm{Tr}_{\textrm{B}}$ denotes the trace over the degrees of freedom of
thermal baths. Moreover, since the interaction term is expressed as
$\hat{V}_{\mu}=\hat{A}_{\mu}{\otimes}\hat{B}_{\mu}$, with
$\hat{A}_{\textrm{R}}=(\hat{a}^{\dagger}+\hat{a})$,
$\hat{A}_{\textrm{Q}}=\hat{\sigma}_{x}$, and
$\hat{B}_{\mu}=\sum_{k}(g_{k,\mu}\hat{b}^{\dagger}_{k,\mu}+g^{*}_{k,\mu}\hat{b}_{k,\mu})~{}(\mu=\textrm{R},\textrm{S})$,
quantum master equation is specified as

$\displaystyle\frac{d}{dt}\hat{\rho}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle
i[\hat{\rho}_{\textrm{S}},\hat{H}_{\textrm{S}}]+\sum_{\mu=\textrm{R},\textrm{Q}}\int^{\infty}_{0}{d\tau}{\times}$
$\displaystyle\\{C_{\mu}(\tau)[\hat{A}_{\mu}(-\tau)\hat{\rho}_{\textrm{S}},\hat{A}_{\mu}]+\mathrm{h.c.}\\},$

where the correlation function is
$C_{\mu}(\tau)\equiv\langle\hat{B}_{\mu}(\tau)\hat{B}_{\mu}(0)\rangle=\sum_{k}|g_{k,\mu}|^{2}[n_{k,\mu}e^{i\omega_{k,\mu}\tau}+(1+n_{k,\mu})e^{-i\omega_{k,\mu}\tau}]$,
and the Bose-Einstein distribution function is
$n_{k,\mu}=1/[\exp(\omega_{k,\mu}/k_{B}T_{\mu})-1]$. Then, under the
eigenbasis of $\hat{H}_{\textrm{S}}$, i.e.
$\hat{H}_{\textrm{S}}|\psi^{\sigma}_{n}{\rangle}=E^{\sigma}_{n}|\psi^{\sigma}_{n}{\rangle}$,
we expand the operator

$\displaystyle\hat{A}_{\mu}(\tau)=\sum_{n,n^{\prime},\sigma,\sigma^{\prime}}e^{i(E^{\sigma}_{n}-E^{\sigma^{\prime}}_{n^{\prime}})\tau}{\langle}\psi^{\sigma}_{n}|\hat{A}_{\mu}|\psi^{\sigma^{\prime}}_{n^{\prime}}{\rangle}|\psi^{\sigma}_{n}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{n^{\prime}}|.$
(27)

If we reexpress
$\hat{A}_{\mu}(\tau)=\sum_{\omega=E^{\sigma}_{n}-E^{\sigma^{\prime}}_{n^{\prime}}}\hat{P}_{\mu}(\omega)e^{i\omega\tau}$,
the quantum master equation can be obtained as

$\displaystyle~{}\frac{d}{dt}\hat{\rho}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle
i[\hat{\rho}_{\textrm{S}},\hat{H}_{\textrm{S}}]+\frac{1}{2}\sum_{\omega,\omega^{\prime};u=\textrm{R},\textrm{Q}}\\{\kappa_{u}(\omega^{\prime})[\hat{P}_{u}(\omega^{\prime})\hat{\rho}_{s},\hat{P}_{u}(\omega)]$
(28) $\displaystyle+\mathrm{h.c.}\\},$

where the rates are
$\kappa_{u}(\omega^{\prime})=\gamma_{u}(\omega^{\prime})n_{u}(\omega^{\prime})$.
The Eq. (28) is identical with Eq. (8).

From Eq. (28), it is known that the populations are generally coupled with the
off-diagonal elements of $\hat{\rho}_{\textrm{S}}$ at finite-time dynamics.
However, after long-time evolution the off-diagonal elements become negligible
due to the full thermalization of the qubit-resonator hybrid system. Hence,
the populations dynamics are approximately decoupled from the off-diagonal
elements. Specifically, the projectors $\hat{P}_{\mu}(\omega)$ and
$\hat{P}_{\mu}(\omega^{\prime})$ in Eq. (28) are restricted to
$\hat{P}_{\mu}(\omega=E^{\sigma}_{n}-E^{\sigma^{\prime}}_{n^{\prime}})={\langle}\psi^{\sigma}_{n}|\hat{A}_{\mu}|\psi^{\sigma^{\prime}}_{n^{\prime}}{\rangle}|\psi^{\sigma}_{n}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{n^{\prime}}|$
and
$\hat{P}_{\mu}(\omega^{\prime}=E^{\sigma^{\prime}}_{n^{\prime}}-E^{\sigma}_{n})={\langle}\psi^{\sigma^{\prime}}_{n^{\prime}}|\hat{A}_{\mu}|\psi^{\sigma}_{n}{\rangle}|\psi^{\sigma^{\prime}}_{n^{\prime}}{\rangle}{\langle}\psi^{\sigma}_{n}|$.
Consequently, the quantum dressed master equation can be obtained as

$\displaystyle~{}\frac{d}{dt}\hat{\rho}_{\textrm{S}}$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{\textrm{S}},\hat{\rho}_{\textrm{S}}]$ (29)
$\displaystyle+\sum_{m,\sigma}\\{\Gamma^{\textrm{R},+}_{m,\sigma}\mathcal{\hat{D}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma}_{m-1}|]\hat{\rho}_{\textrm{S}}$
$\displaystyle+\Gamma^{\textrm{R},-}_{m,\sigma}\mathcal{\hat{D}}[|\psi^{\sigma}_{m-1}{\rangle}{\langle}\psi^{\sigma}_{m}|]\hat{\rho}_{\textrm{S}}\\}$
$\displaystyle+\sum_{m,m^{\prime},\sigma}\\{\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}\mathcal{\hat{D}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\overline{\sigma}}_{m^{\prime}}|]\hat{\rho}_{\textrm{S}}$
$\displaystyle+\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma}\mathcal{\hat{D}}[|\psi^{\overline{\sigma}}_{m^{\prime}}{\rangle}{\langle}\psi^{\sigma}_{m}|]\hat{\rho}_{\textrm{S}}\\},$

which is the same as Eq. (14). The dressed state dissipator is given by
$\mathcal{\hat{D}}[|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|]\hat{\rho}_{\textrm{S}}=|\psi^{\sigma}_{m}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|\hat{\rho}_{\textrm{S}}|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma}_{m}|-\frac{1}{2}(|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|\hat{\rho}_{\textrm{S}}+\hat{\rho}_{\textrm{S}}|\psi^{\sigma^{\prime}}_{l}{\rangle}{\langle}\psi^{\sigma^{\prime}}_{l}|)$.
The transition rates involved with the resonator-bath interaction are
described as

$\displaystyle\Gamma^{\textrm{R},+}_{m,\sigma}=$ $\displaystyle
m(f_{\sigma}-g_{\sigma})^{2}\gamma_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})n_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma}),$
(30a) $\displaystyle\Gamma^{\textrm{R},-}_{m,\sigma}=$ $\displaystyle
m(f_{\sigma}-g_{\sigma})^{2}\gamma_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})[1+n_{\textrm{R}}(\eta_{\sigma}\omega_{\sigma})].$
(30b)

And the transition rates related with the qubit-bath interaction are given by

$\displaystyle\Gamma^{\textrm{Q},+}_{m,m^{\prime},\sigma}=$
$\displaystyle\theta(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})G^{m^{\prime},\overline{\sigma}}_{m,\sigma}G^{m,\sigma}_{m^{\prime},\overline{\sigma}}\gamma_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})n_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}}),$
(31a) $\displaystyle\Gamma^{\textrm{Q},-}_{m,m^{\prime},\sigma}=$
$\displaystyle\theta(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})G^{m^{\prime},\overline{\sigma}}_{m,\sigma}G^{m,\sigma}_{m^{\prime},\overline{\sigma}}\gamma_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})[1+n_{\textrm{Q}}(E^{m,\sigma}_{m^{\prime},\overline{\sigma}})],$
(31b)

where the positive energy gap is
$E^{m,\sigma}_{m^{\prime},\overline{\sigma}}=E_{m,{\sigma}}-E_{m^{\prime},{\overline{\sigma}}}$,
and the squeezing state overlap coefficient is defined as
$G^{n,\sigma}_{m,\bar{\sigma}}\equiv{\langle}\psi^{\sigma}_{n}|\hat{\sigma}_{x}|\psi^{\bar{\sigma}}_{m}{\rangle}$,
which is specified as

$\displaystyle G^{m,\sigma}_{m^{\prime},\bar{\sigma}}$ $\displaystyle=$
$\displaystyle\frac{({v_{\sigma}}/{2u_{\sigma}})^{m/2}}{(m!m^{\prime}!u_{\sigma})^{1/2}}\sum^{\min\\{m,m^{\prime}\\}}_{l=0}\frac{m^{\prime}!H_{m^{\prime}-l}}{l!(m^{\prime}-l)!}\frac{m!H_{m-l}}{(m-l)!}$
(32)
$\displaystyle{\times}\left(\frac{2}{u_{\sigma}v_{\sigma}}\right)^{l/2}\left(\frac{-v^{*}_{\sigma}}{2u_{\sigma}}\right)^{(m^{\prime}-l)/2},$

with the coefficients $H_{m}=(-1)^{m/2}m!/(m/2)!$ for an even $m$ and
$H_{m}=0$ for an odd $m$,
$u_{\sigma}=\cosh(\alpha_{\sigma}-\alpha_{\overline{\sigma}})$, and
$v_{\sigma}=-\sinh(\alpha_{\sigma}-\alpha_{\overline{\sigma}})$.

## References

  * (1) G. Chen, _Nanoscale energy transport and conversion_ (Oxford University Press, England, 2005).
  * (2) G. Benenti, G. Casati, K. Saito, and R. S. Whitney, Phys. Rep. 694, 1 (2017).
  * (3) K. Ptaszyński and M. Esposito, Phys. Rev. Lett. 122, 150603 (2019).
  * (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) J. Senior, A. Gubaydullin, B. Karimi, J. T. Peltonen, J. Ankerhold, and J. P. Pekola,Communication Physics 3, 40 (2020).
  * (6) 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).
  * (7) H. T. Xu, D. Mason, L. Y. Jiang, and J. G. E. Harris, Nature 537, 80 (2016).
  * (8) S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi, and E. Waks, Science 359, 666 (2018).
  * (9) Y. J. Wang, J. Ren, W. X. Zhang, L. He, and X. D. Zhang, Photon. Res. 8, B39 (2020).
  * (10) U. Weiss, _Quantum dissipative systems_ (World Scientific, Singapore, 2008).
  * (11) D. Segal, Phys. Rev. B 73, 205415 (2006).
  * (12) T. Chen, X. B. Wang, and J. Ren, Phys. Rev. B 87, 144303 (2013).
  * (13) K. Saito and T. Kato, Phys. Rev. Lett. 111, 214301 (2013).
  * (14) M. Carrega, P. Solinas, M. Sassetti, and U. Weiss, Phys. Rev. Lett. 116, 240403 (2016).
  * (15) W. J. Dou, M. A. Ochoa, A. Nitzan, and J. E. Subotnik, Phys. Rev. B 98, 134306 (2018).
  * (16) H. Maguire, J. lles-Smith, and A. Nazir Phys. Rev. Lett. 123, 093601 (2019).
  * (17) A. J. Leggett, S. Chakravarty, A. T. Dorsey, Matthew P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * (18) D. Segal and A. Nitzan, Phys. Rev. Lett. 94, 034301 (2005).
  * (19) C. Wang, J. Ren, and J. S. Cao, Sci. Rep. 5 11878 (2015).
  * (20) C. Wang, J. Ren, and J. S. Cao, Phys. Rev. A 95 023610 (2017).
  * (21) X. F. Cao, C. Wang, H. Zheng, and D. H. He, Phys. Rev. B 103, 075407 (2021).
  * (22) Y. Yang and C. Q. Wu, Euro. Phys. Lett. 107, 30003 (2014).
  * (23) J. J. Liu, H. Xu, B. Li, and C. Q. Wu, Phys. Rev. E 96, 012135 (2017).
  * (24) M. Wallquist, K. Hammerer, P. Rabl, M. Liukin, and P. Zoller, Phys. Scr. T137, 014001 (2009).
  * (25) G. Kurizki, P. Bertet, Y. Kubo, K. Molmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, PNAS 112, 3866 (2015)
  * (26) A. A. Clerk, K. W. Lehnert, P. Bertet, J. R. Petta, and Y. Nakamura, Nat. Phys. 16, 257 (2020).
  * (27) A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * (28) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * (29) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * (30) J. lles-Smith, N. Lambert, and A. Nazir, Phys. Rev. A 90, 032114 (2014).
  * (31) J. lles-Smith, A. G. Dijkstra, N. Lambert, and A. Nazir, J. Chem. Phys. 144, 044110 (2016).
  * (32) D. Newman, F. Mintert, and A. Nazir, Phys. Rev. E 95, 032139 (2017).
  * (33) S. Restrepo, J. Cerrillo, P. Strasberg, and G. Schaller, New J. Phys. 20, 053063 (2018).
  * (34) C. McConnell and A. Nazir, J. Chem. Phys. 151, 054104 (2019).
  * (35) S. Felicetti, D. Z. Rossatto, E. Rico, E. Solano, and P. Forn-Díaz, Phys. Rev. A 97, 013851 (2018).
  * (36) S. Felicetti, M. J. Hwang, and A. LeBoité, Phys. Rev. A 98, 053859 (2018).
  * (37) Q. H. Chen, C. Wang, S. He, T. Liu, and K. L. Wang, Phys. Rev. A 86, 023822 (2012).
  * (38) L. W. Duan, Y. F. Xie, D. Braak, and Q. H. Chen, J. Phys. A 49, 464002 (2016).
  * (39) L. Garbe, I. L. Egusquiza, E. Solano, C. Ciuti, T. Coudreau, P. Milman, and S. Felicetti, Phys. Rev. A 95, 053854 (2017).
  * (40) X. Y. Chen and Y. Y. Zhang, Phys. Rev. A 97, 053821 (2018).
  * (41) S. F. Cui, B. Grémaud, W. A. Guo, and G. G. Batrouni, Phys. Rev. A 102, 033334 (2020).
  * (42) A. Crescente, M. Carrega, M. Sassetti, and D. Ferraro, Phys. Rev. B 102, 245407 (2020).
  * (43) A. Delmonte, A. Crescente, M. Carrega, D. Ferraro, and M. Sassetti, Entropy 23, 612 (2021).
  * (44) N. Piccione, S. Felicetti, and B. Bellomo, arXiv:2007.07844.
  * (45) J. J. Liu and D. Segal, Nano. Lett. 20, 6128 (2020).
  * (46) A. Seif, W. DeGottardi, K. Esfarjani, and M. Hafezi, Nat. Comm. 9, 1207 (2017).
  * (47) Z. Denis, A. Biella, I. Favero, and C. Ciuti, Phys. Rev. Lett. 124, 083601 (2020).
  * (48) W. J. Nie, G. Y. Li, X. Y. Li, A. X. Chen, Y. H. Lan, and S. Y. Zhu, Phys. Rev. A 102, 043512 (2020).
  * (49) C. Yang, X. R. Wei, J. T. Sheng, and H. B. Wu, Nat. Comm. 11, 4656 (2020).
  * (50) M. Majland, K. S. Christensen, and N. T. Zinner, Phys. Rev. B 101, 184510 (2020).
  * (51) C. Wang, L. Q. Wang, and J. Ren, Chin. Phys. Lett. 38, 010501 (2021).
  * (52) C. Wang, L. Q. Wang, and J. Ren, Chin. Phys. B 30, 030506 (2021).
  * (53) Y. J. Zhao, Y. L. Liu, Y. X. Liu, and F. Nori, Phys. Rev. A 91, 053820 (2015).
  * (54) X. Wang, A. Miranowicz, H. R. Li, and F. Nori, Phys. Rev. A 94, 053858 (2016).
  * (55) X. Wang, A. Miranowicz, H. R. Li, and F. Nori, Phys. Rev. A 96, 063820 (2017).
  * (56) S. Richer and D. Divincenzo, Phys. Rev. B 93, 134501 (2016).
  * (57) N. Lambert, M. Cirio, M. Delbecq, G. Allison, M. Marx, S. Tarucha, and F. Nori, Phys. Rev. B 97, 125429 (2018).
  * (58) T. Yamamoto and T. Kato, arXiv:2011.06815.
  * (59) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010).
  * (60) A. Settineri, V. Macrí, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta, Phys. Rev. A 98, 053834 (2018).
  * (61) P. Kral, J. Mod. Opt. 37, 889 (1990)
  * (62) L. Levitov and G. Lesovik, JETP Lett. 55, 555 (1992).
  * (63) L. Levitov, H. Lee, and G. Lesovik, J. Math. Phys. 37, 4845 (1996).
  * (64) M. Esposito, U. Harbola, and S. Mukamel, Rev. Mod. Phys. 81, 1665 (2009).
  * (65) B. Li, L. Wang, and G. Casati, Appl. Phys. Lett. 88, 143501 (2006).
  * (66) N. B. Li, J. Ren, L. Wang, G. Zhang, P. Hanggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012).
  * (67) J. Ren and J. X. Zhu, Phys. Rev. B 87, 241412(R) (2013).
  * (68) A. Fornieri, G. Timossi, R. Bosisio, P. Solinas, and F. Giazotto, Phys. Rev. B 93, 134508 (2016).
  * (69) B. Li, L. Wang, and G. Casati, Phys. Rev. Lett. 93, 184301 (2004).
  * (70) M. J. M. Pérez, A. Fornieri, and F. Giazotto, Nat. Nanotech. 10, 303 (2015).
  * (71) L. F. Zhang, Y. H. Yan, C. Q. Wu, J. S. Wang, and B. Li, Phys. Rev. B 80, 172301 (2009).
  * (72) F. Zhan, S. Denisov, and P. Hänggi, Phys. Rev. B 84, 195117 (2011).
  * (73) R. Belousov, E. G. D. Cohen, C. S. Wong, J. A. Goree, and Y. Feng, Phys. Rev. E 93, 042125 (2016).
  * (74) N. A. Sinitsyn and I. Nemenman, Europhys. Lett. 77, 58001 (2007).
  * (75) J. Ren, P. Hänggi, and B. Li, Phys. Rev. Lett. 104, 170601 (2010).
  * (76) D. Stoler, Phys. Rev. D 1, 3217 (1970).
  * (77) D. Stoler, Phys. Rev. D 4, 1925 (1971).
  * (78) U. L. Andersen, T. Gehring, C. Marquardt, and G. Leuchs, Phys. Scr. 91, 053001 (2016).
  * (79) J. R. Kukliński and J. L. Madajczyk, Phys. Rev. A 37, 3175 (1988).
  * (80) A. Xuereb, Nat. Phys. 16, 707 (2020).
  * (81) A. Szorkovszky, G. A. Brawley, A. C. Doherty, and W. P. Bowen, Phys. Rev. Lett. 110, 184301 (2013).
  * (82) A. Kronwald, F. Marquardt, and A. A. Clerk, Phys. Rev. A 88, 063833 (2013).
  * (83) E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Science 349, 952 (2015).
  * (84) J. Ma, X. G. Wang, C. P. Sun, and F. Nori, Phys. Rep. 509, 89 (2011).
  * (85) J. Q. Liao and C. K. Law, Phys. Rev. A 83, 033820 (2011).
  * (86) X. Y. Lü, J. Q. Liao, L. Tian, and F. Nori, Phys. Rev. A 91, 013834 (2015).
  * (87) M. Marthaler, G. Schön, and A. Shnirman, Phys. Rev. Lett. 101, 147001 (2008).
  * (88) J. K. Xie, S. L. Ma, Y. L. Ren, X. K. Li, and F. L. Li, Phys. Rev. A 101, 012348 (2020).