Nonequilibrium thermal transport and photon squeezing in a quadratic qubit-resonator 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)

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.

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

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