# Dicke superradiant heat current enhancement in circuit quantum electrodynamics

Gian Marcello Andolina, Paolo Andrea Erdman, Frank Noé, Jukka Pekola, Marco Schirò

► To cite this version:

Gian Marcello Andolina, Paolo Andrea Erdman, Frank Noé, Jukka Pekola, Marco Schirò. Dicke superradiant heat current enhancement in circuit quantum electrodynamics. Physical Review Research, 2024, 6 (4), pp.043128. 10.1103/physrevresearch.6.043128. hal-04794694

HAL Id: hal-04794694

https://hal.science/hal-04794694v1

Submitted on 21 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dicke superradiant heat current enhancement in circuit quantum electrodynamics

Gian Marcello Andolina,¹ Paolo Andrea Erdman,² Frank Noé,³,²,⁴,⁵ Jukka Pekola,⁶ and Marco Schirò¹

¹JEIP, UAR 3573 CNRS, Collège de France, PSL Research University, 11 Place Marcelin Berthelot, F-75321 Paris, France

²Freie Universität Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany

³Microsoft Research AI4Science, Karl-Liebknecht Strasse 32, 10178 Berlin, Germany

⁴Freie Universität Berlin, Department of Physics, Arnimallee 6, 14195 Berlin, Germany

⁵Rice University, Department of Chemistry, Houston, Texas 77005, USA

⁶Pico Group, QTF Centre of Excellence, Department of Applied Physics, Aalto University, P.O. Box 15100, FI-00076 Aalto, Finland

<figure><img src="image_2.png" /></figure>

(Received 5 March 2024; accepted 23 September 2024; published 12 November 2024)

Controlling the flow of energy and heat at the microscale is crucial to achieve energy-efficient quantum technologies, for on-chip thermal management, and to realize quantum heat engines and refrigerators. Yet, the efficiency of current quantum technologies is affected by thermal noise, and efficient cooling of quantum devices remains challenging in various solid-state implementations such as superconducting circuits. Collective effects, such as the Dicke superradiant emission, have been exploited to enhance the performance of quantum devices. However, the inherently transient nature of Dicke superradiant emission raises questions about its impact on steady-state properties. Here, we study how to enhance the steady-state heat current flowing between a hot bath and a cold bath through an ensemble of *N* qubits, that are collectively coupled to the thermal baths. Remarkably, we find a regime where the heat current scales quadratically with *N* in a finite-size scenario. Conversely, when approaching the thermodynamic limit, we prove that the collective scenario exhibits a parametric enhancement over the noncollective case. We then consider the presence of a third uncontrolled *parasitic* bath, interacting locally with each qubit, that models unavoidable couplings to the external environment. Despite having a nonperturbative effect on the steady-state currents, we show that the collective enhancement is robust to such an addition. Finally, we discuss the feasibility of realizing such a Dicke heat valve with superconducting circuits. Our findings indicate that in a minimal realistic experimental setting with two superconducting qubits, the collective advantage offers an enhancement of approximately 10% compared to the noncollective scenario.

DOI: 10.1103/PhysRevResearch.6.043128

I. INTRODUCTION

Quantum thermodynamics [1–6] is the study of heat and work management in quantum systems. Within the recent blooming of quantum technologies, the critical aspect of energy management becomes increasingly crucial [7] due to the natural interest in building energy-efficient quantum technologies and limiting associated energy waste. In this context, it is particularly relevant to investigate quantum heat transport in superconducting circuit quantum electrodynamics (circuit QED) [8], which is among the most promising platforms for quantum technologies and quantum computation [9]. The heat current flowing across circuit QED devices has been recently measured across various designs [10–14], and circuit QED is emerging as a platform to realize quantum heat engines and refrigerators [15–27].

Two recent circuit QED experiments have investigated the heat transport of microwave photons scattering off a single qubit, realizing a photonic heat valve [28,29]. In these experiments, a single transmon qubit was capacitively coupled to two microwave resonators, each of them in contact with a

resistance acting as a heat bath. A natural question is therefore whether the performance of such a device could be enhanced in the presence of multiple qubits interacting among each other. Indeed, the transport of heat and energy are well known to be sensitive probes of collective and many-body effects [30–35]. A notable example where collective phenomena result in an enhanced emission and superextensive scaling is provided by the Dicke model [36], where an ensemble of *N* atoms in an optical cavity collectively radiates with a superextensive intensity that scales as *N*², i.e., enhanced by a factor *N* with respect to ordinary fluorescence, where atoms emit independently. In the Dicke model, the electrical dipoles of the atoms synchronize thanks to their collective coupling to the optical cavity modes, leading to an enhanced emission which has been dubbed “superradiance” [37,38]. Superradiant emission has been observed in various systems, including Rydberg atoms in cavities [39], color centers in diamonds [40], and superconducting qubits [41].

Collective effects, including Dicke superradiance, have been proposed to improve the performance of thermometers [42–44], quantum heat engines [45–56], quantum batteries [57–61], refrigerators [62], heat transfer [63], and energy transfer between quantum devices [64] and to reduce the dissipated work in finite-time thermodynamics [65].

In this work, we investigate the role of collective superradiant coupling in enhancing the heat current flowing between

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. This figure shows an illustration of our system, where a heat current $\dot{Q}_h$ flows from a hot bath, characterized by a temperature $T_h$, to a cold bath with a temperature $T_c$. This heat current is mediated by an ensemble of $N$ qubits, that are not directly coupled. The collective coupling to the baths leads to a collective enhancement of the heat current.</figcaption></figure>

two thermal baths through an ensemble of artificial atoms (qubits), as depicted in Fig. 1. This is a nontrivial question, particularly since superradiance is typically a transient phenomenon observed in emission and absorption, whereas our focus is on a nondriven, steady-state situation. Remarkably, we find that even under these conditions, the collective coupling between the qubits and the thermal bath can lead to the emergence of superextensive scaling in the heat current for systems of finite size. While the superextensive scaling vanishes in the thermodynamic limit, we find that the heat current is nonetheless enhanced by a temperature-dependent prefactor that diverges as the temperature increases.

While many previous proposals have examined the influence of collective effects on thermal machines [42–62,64] and heat transport [63], it remains crucial to discuss the robustness of these effects to noise. Noise can hinder the coherence of the dynamics, potentially undermining the collective enhancement. Here, we show that this collective enhancement is robust to the addition of a third uncontrolled *parasitic* bath, which interacts locally with each qubit. This finding is crucial since, as we will show, the presence of an infinitesimally small local noise has a finite and nonperturbative effect on the steady-state heat currents, and removes the dependence of the steady-state heat current from the initial state preparation. Furthermore, we demonstrate the resilience of the superradiant effect in this realistic noisy environment.

Despite the intense theoretical interest in exploiting collective effects in quantum thermodynamical machines [42–58,60–64], only a few experiments have been conducted in this context [59,66]. Hence, we discuss the feasibility of measuring such a Dicke enhancement of the heat current in an experimental setup. Such a device can be realized within the framework of circuit QED [8,28,67,68], where $N$ transmons are capacitively coupled to two *RLC* circuits where the dissipative nature of a thermal bath stems from the presence of the resistive elements. Our findings indicate that using experimentally realistic parameters in the minimal case of $N = 2$ superconducting qubits, the collectively

<figure><img src="image_4.png" /><figcaption>FIG. 2. This figure shows a lumped-element circuit diagram of a circuit-QED realization of this system. The circuit features two *LC* resonators, with associated voltages $V_i$ where $i = h, c$. Each resonator consists of an inductance $L$ and a capacitance $C$, and is coupled to a resistance $R_i$. The ensemble of $N$ transmons is represented by fluxes $\phi_j$ for $j = 1, ..., N$, and their time derivatives $\dot{\phi}_j$. Each transmon is made up of a Josephson junction with an associated inductance $L_J$ and a capacitance $C_T$. The transmons are capacitively coupled to the resonators via capacitors $C_c^{(i)}$.</figcaption></figure>

enhanced heat current is approximately 10% higher than in scenarios where collective effects are absent. This highlights the potential of leveraging collective quantum behaviors in practical thermodynamic applications, offering a measurable enhancement over more traditional designs, and providing a platform to experimentally observe superradiant effects in measurable steady-state heat currents. Our proposed device also represents a many-body collective version of a heat valve, a device that has been previously implemented in the context of circuit QED [28]. A heat valve is a device designed to control the flow of heat between two baths. Indeed, the circuit pictured in Fig. 2 can function as a heat valve. The *LC* circuits in Fig. 2 (left and right elements) act as a filter; when the qubits are tuned to resonate at the *LC* frequency, there is an efficient flow of heat between the two baths. Conversely, if the qubits are detuned from the *LC* frequency, the heat flow is effectively impeded due to the frequency mismatch. This dynamic tuning capability allows for controlled manipulation of heat transfer, embodying the core function of a heat valve.

This article is organized as follows, in Sec. II, we derive the model and discuss the superextensive behavior of the heat current in the noiseless case, i.e., when the modulus of the collective spin operator is conserved. In Sec. III, we discuss the resilience of this superextensive behavior to the addition of a realistic parasitic bath, breaking the conservation of the modulus of the collective spin. In Sec. IV, we discuss the experimental feasibility of observing our findings and we propose a minimal experimental setup that exhibits the superextensive behavior of the heat current. In Sec. V, we draw our conclusions. Appendices A–D contain a series of technical details.
## II. THE DICKE SUPERRADIANT HEAT VALVE

The lumped-element circuit diagram, including capacitances, inductances, resistances, and the various variables, is depicted in Fig. 2. The quantum heat valve under study consists of N transmon qubits coupled to two resistors denoted by $R_i$ in parallel with LC resonators. These RLC circuits form two heat baths characterized by their temperature $T_i$. We denote the two baths as “hot” (h) and “cold” (c), hence $i = \{h, c\}$ and $T_h \ge T_c$. As derived in Appendix A, the total Hamiltonian of the system can be written as

$$ \hat{H}_{\text{tot}} = \hat{H}_{0} + \hat{H}_{\text{int}} + \hat{H}_{\text{RLC}}. \quad (1) $$

$\hat{H}_0$ and $\hat{H}_{\text{int}}$, respectively, describe the Hamiltonian of the qubits and their coupling to the RLC circuit. Before projecting onto the qubit subspace, the Hamiltonian of the circuit is given by

$$ \begin{aligned} \hat{H}_0 &= \frac{1}{2} \sum_{j=1}^{N} \left[ \frac{1}{C_{\text{eff}}} \hat{Q}_j^2 - E_j \cos \left( \frac{2\pi \hat{\phi}_j}{\phi_0} \right) \right], \\ \hat{H}_{\text{int}} &= - \sum_{j=1}^{N} \sum_{i=h,c} \frac{C_c^{(i)}}{C_T} \hat{Q}_j \hat{V}_i, \end{aligned} \quad (2) $$

where $C_{\text{eff}} = C_T[1 - \sum_{i=h,c} C_c^{(i)}/C_T]^{-1}$, $\hat{Q}_j$, $\hat{\phi}_j$ are the charge and flux associated with the jth transmon, $C_T$ ($E_j$) is the transmon capacitance (Josephson energy), $\phi_0 = h/(2e)$, $C_c^{(i)}$ are the coupling inductances, and $\hat{V}_i$ is the voltage of the ith RLC circuit (for $i = h, c$), as reported in the lumped-circuit diagram in Fig. 2. By projecting the Hilbert space onto the two-level subspace of each transmon, we obtain

$$ \hat{H}_0 = \sum_{j=1}^{N} \frac{\hbar \omega_0}{2} \hat{\sigma}_z^{(j)}, \quad (3) $$

where $\hat{\sigma}_z^{(j)}$ represents the Pauli z operator for the jth qubit and $\hbar\omega_0$ is the qubit energy, and

$$ \hat{H}_{\text{int}} = - \sum_{j=1}^{N} \sum_{i=h,c} \hbar G_i \frac{\hat{\sigma}_x^{(j)}}{2} \hat{V}_i, \quad (4) $$

where $\hat{\sigma}_x^{(j)}$ is the dipole operator of the jth qubit (the Pauli x operator) and $G_i$ is the coupling strength that we consider different for the two RLC circuits, but otherwise uniform across the qubits. This coupling plays a crucial role in the interaction dynamics of the qubits with the thermal baths. Finally, $\hat{H}_{\text{RLC}}$ describes the Hamiltonian of the two RLC circuits as described in Appendix A.

It is useful to introduce a collective spin operator $\hat{J}_\alpha$, with $\alpha = x, y, z$, as

$$ \hat{J}_{\alpha} = \sum_{j=1}^{N} \frac{\hat{\sigma}_{\alpha}^{(j)}}{2}, \quad (5) $$

where $\hat{\sigma}_\alpha^{(j)}$ are the Pauli operators for the jth qubit. The total Hamiltonian in Eq. (1) can thus be rewritten (up to an additive constant) as

$$ \hat{H}_{\text{tot}} = \hbar\omega_0\hat{J}_z - \sum_{i=h,c} \hbar G_i \hat{V}_i + \hat{H}_{\text{RLC}}. \quad (6) $$

This Hamiltonian conserves the norm of the collective spin operator $\hat{J}^2$, where

$$ \hat{J}^2 \equiv \sum_{\alpha=x,y,z} \hat{J}_{\alpha}^2. \quad (7) $$

It is thus useful to introduce the Dicke states, given by

$$ |J, m_J\rangle = \sqrt{\frac{(J+m_J)!(J-m_J)!}{(2J)!}} (\hat{J}_+)^{J+m_J} |J, -J\rangle, \quad (8) $$

where $J(J+1) (m_J)$ are the eigenvalues of $\hat{J}^2$ ($\hat{J}_z$), and $\hat{J}_+$ ($\hat{J}_-$) is the raising (lowering) operator.

Since the total Hamiltonian given by Eq. (6) commutes with $\hat{J}^2$, the system dynamics is confined to a subspace characterized by a fixed value of J. The allowed values for J are non-negative and follow the sequence $J = N/2, (N-1)/2, (N-2)/2, \dots$. Correspondingly, $m_J$ can range from $-J$ to J in integer steps, namely, $m_J = -J, -J+1, \dots, J-1, J$. As an example, consider a system of N = 2 qubits. Here, the total collective spin J can either be 1 or 0. This means that the system can either be in a triplet state ($J=1$) with $m_J = -1, 0, 1$ or in a singlet state ($J=0$) with $m_J = 0$.

Motivated by existing experimental setups [28], we are interested in a scenario where the qubits are weakly coupled to the baths. In this regime, it is possible to invoke the Born-Markov approximation and describe the open system dynamics of the qubits with a suitable Lindblad master equation,

$$ \frac{d\hat{\rho}}{dt} = -\frac{i}{\hbar}[\hat{H}_0, \hat{\rho}] + \mathcal{D}(\hat{\rho}), \quad (9) $$

where $\hat{\rho}$ is the density operator of the qubits. A detailed derivation of this Lindblad master equation is carried out in Appendix B, while here we summarize the results needed for our analysis. The dissipator $\mathcal{D}(\hat{\rho})$ arises from integrating out the RLC circuits which are coupled to our system through the voltage $\hat{V}_i$ [see Eqs. (4) and (6)]. Within the Born-Markov weak-coupling approximation, the dynamics of the voltage is fully encoded in the voltage dynamical structure factor [69] $S_{\hat{V}_i, \hat{V}_j}(\omega)$, defined in Eq. (B10), which describe the voltage fluctuations occurring in the RLC circuits.

The interaction term in Eq. (6), coupling the qubits system with the two resonators, is proportional to $\hat{J}_x = \hat{J}_+ + \hat{J}_-$. Hence, after performing the secular approximation as needed to derive a Lindblad master equation [70], the suitable Lindblad operators for our system, describing transitions between the Dicke states induced by the exchange of energy with thermal baths, can be conveniently written using the collective spin raising and lowering operators $\hat{J}_+$ and $\hat{J}_-$. Hence, the total dissipator in Eq. (9) can be separated into contributions from the hot and cold baths,

$$ \mathcal{D}(\hat{\rho}) = \mathcal{D}_h(\hat{\rho}) + \mathcal{D}_c(\hat{\rho}), \quad (10) $$

given by

$$ \begin{aligned} \mathcal{D}_i(\hat{\rho}) ={}& \gamma_i(1 + n_i)(\hat{J}_-\hat{\rho}\hat{J}_+ - \frac{1}{2}\{\hat{J}_+\hat{J}_-, \hat{\rho}\}) \\ & + \gamma_i n_i(\hat{J}_+\hat{\rho}\hat{J}_- - \frac{1}{2}\{\hat{J}_-\hat{J}_+, \hat{\rho}\}), \end{aligned} \quad (11) $$

where $i = h, c$ and $\gamma_i$ and $\gamma_h$ are the transition rates for the cold and hot baths, respectively. The rates $\gamma_i$, that represent
the coupling strength between the baths and the qubits can be expressed in terms of the microscopic parameters of the circuit, as shown in Appendix B. The thermal occupation numbers $n_i$ are given by the Bose-Einstein distribution,

$$n_i = \frac{1}{\exp\left(\frac{\hbar\omega_0}{k_B T_i}\right) - 1}, \quad (12)$$

where $k_B$ is the Boltzmann's constant.

To summarize, after integrating out the two RLC circuits representing the two thermal baths, we have obtained an effective description where the system of N qubits is in contact with two structured thermal baths at finite temperature, each providing injection and depletion of energy through global spin-flip processes as described by Eq. (11). It is important to notice that the Lindblad master equation in Eq. (9) correctly inherits the conservation of $\hat{J}^2$ from the total Hamiltonian in Eq. (6). Specifically, this becomes a strong symmetry of the Lindbladian superoperator since $\hat{J}^2$ commutes both with $H_0$ as well as with each jump operator $\hat{J}_{\pm}$. As such, this corresponds to a conserved quantity of the Lindblad dynamics [71]. Interestingly, analogous Lindblad master equations have been derived to describe dissipative (boundary) time crystals [72], where usually one considers a single bath at zero temperature.

### A. Steady-state density matrix and effective temperature

Since we are interested in computing the steady-state heat current flowing across the device, we start discussing the steady-state solution of the Lindblad master equation (9). In the case of a single reservoir, i.e., $\mathcal{D}(\hat{\rho}) = \mathcal{D}_i(\hat{\rho})$ for $i = h$ or $i = c$, the steady-state solution corresponds to the thermal distribution (within each subspace with fixed $J$) characterized by the corresponding temperature $T_i$ since the rates satisfy the detailed balance condition,

$$\frac{\gamma_i n_i}{\gamma_i(1+n_i)} = \exp\left(-\frac{\hbar\omega_0}{k_B T_i}\right), \quad (13)$$

capturing the equilibrium relation between absorption and emission rates in each bath. This result, consistent with the laws of thermodynamics, comes from using a master equation written in terms of the global jump operators for the coupled qubits (see Refs. [73–75] for a discussion of this point). In the case of two baths, the total dissipator in Eq. (10) can be expressed as

$$\mathcal{D}(\hat{\rho}) = [\gamma_c(1+n_c) + \gamma_h(1+n_h)](\hat{J}_{-}\hat{\rho}\hat{J}_{+} - \frac{1}{2}\{\hat{J}_{+}\hat{J}_{-}, \hat{\rho}\}) \\ + [\gamma_c n_c + \gamma_h n_h](\hat{J}_{+}\hat{\rho}\hat{J}_{-} - \frac{1}{2}(\hat{J}_{-}\hat{J}_{+}, \hat{\rho})]. \quad (14)$$

For clarity in this discussion, we introduce a temperature $T_0$ associated with the qubit frequency $\omega_0$, $T_0 \equiv (\hbar\omega_0)/k_B$. The dissipator in Eq. (14) defines a detailed balance equation,

$$\frac{\gamma_c n_c + \gamma_h n_h}{\gamma_c(1+n_c) + \gamma_h(1+n_h)} = \exp\left(-\frac{T_0}{T^*}\right), \quad (15)$$

which, in turn, defines an effective temperature $T^*$, such that $T_c \le T^* \le T_h$. We can invert the previous equation to get an

explicit expression for the effective temperature,

$$T^* = \frac{T_0}{\ln\left[\frac{\gamma_c(1+n_c)+\gamma_h(1+n_h)}{\gamma_c n_c + \gamma_h n_h}\right]}. \quad (16)$$

The steady-state solution $\hat{\rho}^{(s)}$, determined setting $\mathcal{D}(\hat{\rho}^{(s)}) = 0$, is thus a thermal state at temperature $T^*$ within each subspace with fixed $J$. However, since $\hat{J}^2$ is conserved, the relative occupation of different subspaces is fixed by the initial state. Assuming that the initial state does not contain any coherence between subspaces with different $J$, the stationary state $\hat{\rho}^{(s)}$ is given by

$$\hat{\rho}^{(s)} = \sum_{J} P_J \sum_{m_J} P(m_J|J) |J, m_J\rangle\langle J, m_J|, \quad (17)$$

where $m_J$ are the eigenvalues of $\hat{J}_z$ compatible with $J$, $P_J$ is the occupation probability of subspace $J$ determined by the initial state, and $|J, m_J\rangle$ are the Dicke states, given in Eq. (8). It should be noted that the summation over $J$ in Eq. (17) incorporates the degeneracies associated with $J$, as is customary in angular momentum composition [76]. For example, considering a system composed of three qubits, each characterized by a spin of $1/2$, the resultant collective angular momentum can exhibit values of $J = 1/2, 1/2, 3/2$, with the $J = 1/2$ appearing twice. This is due to the composition rule that three spin-$1/2$ entities combine to yield collective angular momenta of $3/2$ and two instances of $1/2$. This composition rule is customarily represented as $1/2 \otimes 1/2 \otimes 1/2 = 3/2 \oplus 1/2 \oplus 1/2$. The dissipator in Eq. (14) fixes the occupation within each subspace $J$, i.e., it fixes the conditional probabilities $P(m_J|J)$ to be thermal at the effective temperature $T^*$,

$$P(m_J|J) = \frac{\exp\left(-\frac{m_J T_0}{T^*}\right)}{Z_{J,T^*}}, \quad (18)$$

where $Z_{J,T^*} = \sum_{m_J} \exp[-m_J T_0/T^*]$ is the partition function.

### B. Heat current

The general expression for the total heat current $\dot{Q}$ flowing out of the thermal baths can be defined as [5,77]

$$\dot{Q} \equiv \mathrm{Tr}\left[\hat{H}_0 \frac{d\hat{\rho}}{dt}\right] = \mathrm{Tr}[\hat{H}_0 \mathcal{D}(\hat{\rho})]. \quad (19)$$

The second equality has been obtained by enforcing the dynamics of the Lindblad master equation, given by Eq. (9). The total heat current $\dot{Q}$ can be split into two contributions corresponding to the different baths,

$$\dot{Q} = \dot{Q}_h + \dot{Q}_c, \quad (20)$$

where

$$\dot{Q}_i \equiv \mathrm{Tr}[\hat{H}_0 \mathcal{D}_i(\hat{\rho})]. \quad (21)$$

Here, we are interested only in the steady-state dynamics where $\dot{Q} = 0$ and $\dot{Q}_h = -\dot{Q}_c$.

Using the specific forms of the dissipators for the cold and hot baths [in Eq. (11)], as detailed in Appendix C, we arrive
at the following expressions for the heat currents:

$$ \dot{Q}_i = \hbar\omega_0\gamma_i[-(1 + n_i)\langle\hat{J}_+\hat{J}_-\rangle + n_i\langle\hat{J}_-\hat{J}_+\rangle], \quad (22) $$

where $\langle\hat{x}\rangle = \mathrm{Tr}[\hat{\rho}\hat{x}]$. This result is reminiscent of the well-known superradiant energy emission in an ensemble of excited qubits [37,38,78,79]. In such a context, the system is typically considered to be in contact with a single cold bath at zero temperature, corresponding to setting $\gamma_h = 0$ in our case, and is initially prepared in a highly excited state, such as the state with $J = N/2$ and $M = J$. Consequently, the superradiant heat current flowing into the cold bath is given by $\dot{Q}_c = -\hbar\omega_0\gamma_c\langle\hat{J}_+\hat{J}_-\rangle$. The term $\langle\hat{J}_+\hat{J}_-\rangle$ can be expressed in terms of the individual qubit raising and lowering operators, $\hat{\sigma}_+^{(j)}$ and $\hat{\sigma}_-^{(j)}$, as $\sum_{j,j'}\langle\hat{\sigma}_+^{(j)}\hat{\sigma}_-^{(j')}\rangle$. This summation encompasses $N$ local population terms ($j = j'$), and $\sim N^2$ nonlocal coherent terms ($j \neq j'$) that may enable a superextensive scaling effect. However, it is important to notice that this superradiant effect is only present in the transient dynamics. Over time, the qubits will eventually reach thermal equilibrium with the zero-temperature bath, resulting in the cessation of the heat current.

Going back to the case under study, our primary focus is on the steady-state heat current that flows between the two baths. Using that $\dot{Q}_c = -\dot{Q}_h$, we can express the heat current as

$$ \dot{Q}_h = \left(\frac{1}{\gamma_h} + \frac{1}{\gamma_c}\right)^{-1} \left(\frac{1}{\gamma_h}\dot{Q}_h - \frac{1}{\gamma_c}\dot{Q}_c\right) \quad (23) $$

$$ = \frac{2\gamma_h\gamma_c}{\gamma_h + \gamma_c} \hbar\omega_0\langle-\hat{J}_z\rangle(n_h - n_c), \quad (24) $$

where we used that $[\hat{J}_+, \hat{J}_-] = 2\hat{J}_z$. Notice that the case of a single qubit can be obtained from this equation setting $J = 1/2$.

Notably, Eq. (24) reveals that the steady-state heat current does not explicitly depend on the sum of the nonlocal coherent terms that are present in $\langle\hat{J}_+\hat{J}_-\rangle$ and $\langle\hat{J}_-\hat{J}_+\rangle$, which are typically associated with a superextensive behavior in the context of superradiant emission. Instead, the steady-state heat current is controlled by the operator $\hat{J}_z$, which can be represented as a sum of $N$ local terms, $\hat{J}_z = \sum_{j=1}^N \hat{\sigma}_z^{(j)}$. This would suggest the absence of a superextensive scaling due to superradiant effects in this context. Indeed, upper bounding the expectation value $\langle-\hat{J}_z\rangle$ in Eq. (24) with $N/2$, we find the following exact upper bound:

$$ \dot{Q}_h \le \dot{Q}_h^{\max} = \frac{\gamma_h \gamma_c}{\gamma_h + \gamma_c} \hbar \omega_0 N (n_h - n_c), \quad (25) $$

revealing that for large values of $N$, the heat current cannot scale more than linearly in the number of qubit. It is interesting to notice that the bound in Eq. (25), which is saturable in the thermodynamic limit $N \to \infty$ (see Sec. II C), has a “Landauer-Büttiker-like” dependence on the temperature of the baths [80,81] and is analogous to the heat current flowing through a harmonic oscillator coupled to bosonic heat baths [82]—but enhanced by a factor $N$. Intuitively, this can be expected because, in the thermodynamic limit, the collective spin operators effectively become bosonic operators [83] and should thus behave similarly to a harmonic oscillator. However, the rates here are increased by a factor $N$ since the qubits are coupled in parallel.

Surprisingly, despite the linear scaling of $\dot{Q}_h^{\max}$ with $N$, below we show that the heat current $\dot{Q}_h$ can nonetheless exhibit a superextensive scaling in the finite-size regime. Indeed, if the qubit were in a global thermal state, it can be shown from Eq. (24) that the scaling would be linear since the thermal state is a tensor product of local qubit states. However, since $\hat{J}^2$ is conserved, the steady-state density matrix $\hat{\rho}^{(s)}$, as indicated by Eq. (17), may not be purely thermal, and the presence of coherence between different qubits can enable a nonlinear dependence of the heat current on the qubit number $N$.

As a benchmark for our study, we will compare our results to the case where $N$ qubits (with equal frequency $\omega_0$) are independently coupled to a hot bath and a cold bath. This scenario serves as a reference case, allowing for a clear comparison of our collective case described by a collective Lindblad master equation, given by Eq. (9). To compute the heat current associated with a single qubit, we employ the Lindblad master equation in Eq. (9), along with the expression of the steady-state heat current given in Eq. (24), by setting $\bar{J} = 1/2$, corresponding to $N = 1$. As all qubits in our model are identical, the total heat current in this independent case, denoted as $\dot{Q}_h^{\mathrm{ind}}$, is determined by multiplying the heat current of a single qubit by the total number of qubits, $N$. This leads us to the following expression for the total heat current in the independent scenario:

$$ \dot{Q}_h^{\mathrm{ind}} = N\hbar\omega_0 \left[ \frac{\gamma_h\gamma_c}{\gamma_c(2n_c + 1) + \gamma_h(2n_h + 1)} \right] (n_h - n_c). \quad (26) $$

An analogous result for $N = 1$ was previously derived in the spin-boson model [82,84].

### C. Results for the heat current in a fixed subspace

For the remainder of this section, we focus on scenarios where the initial state has a well-defined $J = \bar{J}$, meaning that

$$ P_J = \delta_{J, \bar{J}}, \quad (27) $$

in Eq. (17). Notice that in general, $\bar{J} \le N/2$, and the case $\bar{J} = N/2$ can be implemented by initializing all the qubits in their ground states $\otimes_{j=1}^N |0\rangle_j$, corresponding to the Dicke state $|N/2, -N/2\rangle$ and hence the system's dynamics is constrained to the subspace with $J = N/2$. While this simplified case may not be directly applicable in experimental settings, it is a useful starting point for gaining insights into this problem before exploring more realistic scenarios.

In this setting, the steady-state heat current $\dot{Q}_h$, as given in Eq. (24), can be analytically derived by averaging the current operator over the thermal state at an effective temperature denoted by $T^*$ and given in Eq. (16). In the general case, the detailed calculation of this quantity is provided in Appendix C.

Analogously to the derivation of Eq. (25), we can upper bound the expectation value $\langle-\hat{J}_z\rangle$ in Eq. (24) with the conserved value of $\bar{J}$, yielding a generally tighter bound,

$$ \dot{Q}_h \le \dot{Q}_h^{\mathrm{lim}} = \frac{2\gamma_h\gamma_c}{\gamma_h + \gamma_c} \hbar\omega_0 \bar{J}(n_h - n_c). \quad (28) $$

Interestingly, as derived in Eq. (C10), the heat current saturates this bound in the limit $\bar{J}T_0 \gg T^*$. Physically, this can be interpreted as the low-temperature regime, i.e., when the
effective temperature $T^*$ is considerably lower than $\bar{J}T_0$, or as the thermodynamic limit $N \to \infty$, choosing $\bar{J} = N/2$. Indeed, in these cases, the majority of spins are in the inverted state, corresponding to $\langle \hat{J}_z \rangle \approx -\bar{J}$. Since $\bar{J}$ is, at most, linear in $N$, this reveals that the heat current cannot scale superextensively in the number of qubits both in the low-temperature regime and in the thermodynamic limit.

Nonetheless, in the thermodynamic limit $N \to \infty$, we still find a parametric advantage in the heat current. In this limit, the ratio between the heat current in the collective case $\dot{Q}_h^{\text{lim}}$ in the $\bar{J} = N/2$ case and the heat current in the independent case $\dot{Q}_h^{\text{ind}}$ is given by Eq. (C12), i.e.,

$$ \frac{\dot{Q}_h^{\text{lim}}}{\dot{Q}_h^{\text{ind}}} = \coth \left[ \frac{T_0}{2T^*} \right]. \qquad (29) $$

The ratio $\dot{Q}_h^{\text{lim}}/\dot{Q}_h^{\text{ind}}$ is always greater than 1, indicating that the collective scenario yields a strictly larger heat current than the independent case. In the limit $T_0/T^* \gg 1$, the ratio approaches one. Conversely, in the limit where $T_0/T^* \ll 1$, the ratio becomes $\dot{Q}_h^{\text{lim}}/\dot{Q}_h^{\text{ind}} \approx (2T^*/T_0)$, which can be arbitrarily large. This indicates that even in the thermodynamic limit, collective effects can provide an unbounded parametric enhancement to the heat current.

We now consider the high-temperature regime, i.e., $\bar{J}T_0 \ll T^*$. Using Eq. (C11), we find that

$$ \dot{Q}_h \approx \frac{2 \gamma_h \gamma_c T_0}{3 \gamma_h + \gamma_c T^*} (\hbar\omega_0) \bar{J}(1 + \bar{J})(n_h - n_c). \qquad (30) $$

In this scenario, the heat current exhibits a superlinear scaling with $\bar{J}$. This scaling behavior is derived from the nonlinear dependence of the thermal populations $P(m_{\bar{J}}|\bar{J})$ on $\bar{J}$, given in Eq. (18), as shown in Appendix C. In the case $\bar{J} = N/2$, this corresponds to a superextensive scaling in the number of qubits, $N$. However, it is crucial to recognize that in the thermodynamic limit ($N \to \infty$), the condition $\bar{J}T_0 \ll T^*$ is never met when $\bar{J} = N/2$ since, as $N$ approaches infinity, $\bar{J} = N/2$ eventually surpasses $T^*/T_0$ and the upper limit $\dot{Q}_h^{\text{lim}}$ in Eq. (28), linear in $\bar{J}$, becomes tight.

In Figs. 3(a) and 3(b), we plot the steady-state heat current $\dot{Q}_h$ (normalized by $\bar{J}\hbar\omega_0^2$) (14) focusing on system sizes of experimental relevance, specifically considering a system size up to $N = 16$, which corresponds to $\bar{J} = 8$. For convenience, we fix the temperature of the cold bath to be $T_c = T_0/3$ and tune the ratio $T_0/T_h$, which, according to the discussion of the steady-state's effective temperature, allows tuning the effective temperature across the relevant scale $T_0$ as discussed above. In Fig. 3(a), we see that in the low-temperature regime, $T_0 > 1.5T_h$ (corresponding to $T_0 > T^*$), the ratio $\dot{Q}_h/(\bar{J}\hbar\omega_0^2)$ does not depend on $\bar{J}$, i.e., the heat current scales linearly in $\bar{J}$, and thus in $N$. However, in the high-temperature regime, i.e., $T_0 < 1.5T_h$ (corresponding to $T_0 < T^*$), a superlinear scaling emerges. In Fig. 3(a), we also show as a dotted line the heat current obtained in the thermodynamic limit $\dot{Q}_h^{\text{lim}}$ (taking $\bar{J} = N/2$), as given in Eq. (28), which is also an upper bound to the heat current. While, as expected, a saturation to the extensive limit $\dot{Q}_h^{\text{lim}}$ is seen for low temperatures, there is no saturation for very high temperatures, $T_0/T_h \ll 1$, where the superlinear scaling persists up to $N = 16$. The scaling with the size is explicitly shown in Fig. 3(b), where the

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 3. (a) The heat current as a function of the ratio $T_0/T_h$ (where $T_0 \equiv \hbar\omega_0/k_B$) for various values of $\bar{J}$. The black dashed line indicates the independent current $\dot{Q}_h^{\text{ind}}$, whereas the red dotted curve represents the heat current in the thermodynamic limit, denoted as $\dot{Q}_h^{\text{lim}}$. (b) The heat current is plotted as a function of $\bar{J}$ for various ratios of $T_0/T_h$. The dotted lines in this panel correspond to the heat current in the thermodynamic limit, $\dot{Q}_h^{\text{lim}}$, in the case $\bar{J} = N/2$. Each of these dotted lines corresponds to the same temperature ratio $T_0/T_h$ of the data points of matching color. The temperature of the cold bath is chosen to be $T_0/T_c = 3$ in both panels. Other parameters are $\gamma_h = \gamma_c = \omega_0/10$.</figcaption></figure>

dependence of the ratio $\dot{Q}_h/(\bar{J}\hbar\omega_0^2)$ on $\bar{J}$ is plotted for different temperatures, highlighting the superlinear scaling for small values of $T_0/T_h$. Again, the heat current in the thermodynamic limit $\dot{Q}_h^{\text{lim}}$ is shown as dotted lines of correspondent color. Analogously to Fig. 3(a), for low temperatures $T_0/T_h \ge 0.7$, the heat current saturates the bound given by $\dot{Q}_h^{\text{lim}}$. Instead, in the high-temperature regime, i.e., for $T_0/T_h \le 0.3$ and for realistic system sizes, the heat current always shows a super-linear behavior, far from saturation.

To summarize, for experimentally relevant cases where $N \le 10$, a superextensive scaling of the heat current is always achievable provided the hot bath is sufficiently warm (high-temperature regime). In these scenarios, the value of $\dot{Q}_h^{\text{lim}}$ essentially serves as an upper limit. Below we show that such a temperature regime that yields a superextensive scaling of
the heat current is within experimental reach and, importantly, this behavior remains robust in the presence of noise.

### III. IMPACT OF A PARASITIC BATH

The dynamics governed by the Lindblad master equation in Eq. (9) preserves the norm of the collective spin operator $\hat{J}^2$. Consequently, the steady-state current depends on the probabilities $P_J$, which are fixed by the initial state as described in Eq. (17). However, a dependence of the steady-state heat current on the initial state preparation of the qubits is not realistic in practice. Although short-term dynamics could be influenced by the initial state preparation, the long-term steady state is likely to be dominated by factors such as noise and local dissipation. These effects would eventually erase the memory of the initial state. Indeed, as we now show, such parasitic effects have a large impact on the steady-state heat current, even in the limit of vanishing small parasitic coupling, thus a nonperturbative effect. To address this effect, we introduce a third, parasitic thermal bath that interacts locally with each qubit. Unlike the collective interactions from the primary thermal baths, this parasitic bath interacts uniformly but locally with each qubit, thereby breaking the conservation of $\hat{J}^2$. This parasitic bath is characterized by a temperature $T_p$ and a coupling strength $\gamma_p$.

The dissipator associated with the parasitic bath, $D_p(\hat{\rho})$, can be expressed as

$$D_p(\hat{\rho}) = \sum_{j=1}^{N} D_p^{(j)}(\hat{\rho}), \qquad (31)$$

and

$$D_p^{(j)}(\hat{\rho}) = \gamma_p n_p (\hat{\sigma}_+^{(j)} \rho \hat{\sigma}_-^{(j)} - \frac{1}{2} \{\hat{\sigma}_-^{(j)} \hat{\sigma}_+^{(j)}, \rho\}) \\ + \gamma_p (1 + n_p) (\hat{\sigma}_-^{(j)} \rho \hat{\sigma}_+^{(j)} - \frac{1}{2} \{\hat{\sigma}_+^{(j)} \hat{\sigma}_-^{(j)}, \rho\}). \quad (32)$$

Here, $n_p = 1/\{\exp[\hbar\omega_0/(k_B T_p)] - 1\}$ represents the mean thermal occupation number of the parasitic bath. The overall dissipative dynamic of the system influenced by all three baths is thus described by

$$\frac{d\hat{\rho}}{dt} = -\frac{i}{\hbar}[\hat{H}_0, \hat{\rho}] + D_h(\hat{\rho}) + D_c(\hat{\rho}) + D_p(\hat{\rho}). \quad (33)$$

Although we assume that this parasitic bath is weakly coupled to the system compared to the primary baths ($\gamma_p \ll \gamma_h, \gamma_c$), its main effect is breaking the conservation of $\hat{J}^2$. Therefore, while the coupling strength $\gamma_p$ can be neglected when calculating the population within a given subspace, its influence is nonperturbative in determining the steady-state probabilities $P_J$ of different subspaces. Hence, once we extract the values of $P_J$ from the numerics, we can still use the results of Sec. II C to determine the occupations $P(m_J|J)$ inside a given subspace.

In Figs. 4(a) and 4(c), we show how the heat current $\dot{Q}_h$, divided by $N$, is influenced by the parasitic bath. Figure 4(a) focuses on the low-temperature regime, i.e., when both primary thermal baths have temperatures lower than $T_0$, with $T_0/T_h = 1$ and $T_0/T_c = 3$. This corresponds to the system being effectively colder than $T_0$, $T_0/T^* \approx 1.42$. As we can see, a hot parasitic bath $T_0/T_p \le 1$ (blue and orange dots) suppresses the heat current even below the independent case in Eq. (26)

(plotted as a black dashed line). In the particular case where $T_p \approx T^* \approx 1.42T_0$ (green dots), the heat current coincides with the value derived for independent qubits, according to Eq. (26). Under this condition, the parasitic bath has the only effect of suppressing the coherence between different qubits, rendering them effectively independent. However, when the parasitic bath is sufficiently cold, $T_0/T_p \ge 10$ (red and purple dots), the heat current $\dot{Q}_h^{\text{ind}}$ coincides with the one observed when the dynamics are constrained to a single subspace with $J = N/2$ (represented by black triangles). Even in this single subspace, a superextensive scaling in $N$ is absent. This absence is attributed to the low-temperature regime of the baths, which maintains the system close to its ground state, leading to $\langle -\hat{J}_z \rangle \approx N/2$ in Eq. (24), which corresponds to a linear behavior.

Figure 4(c) illustrates the high-temperature scenario where $T_0/T_h = 1/3$ and $T_0/T_c = 3$. In this case, the system has an effective temperature bigger than $T_0$, $T_0/T^* \approx 0.57$. Similar to the case in Fig. 4(a), a hot parasitic bath (blue dots) negatively impacts the heat current. If $T_p \approx T^* \approx 0.57T_0$ (orange dots), the heat current reverts to the value computed for independent atoms using Eq. (26). However, when the parasitic bath is sufficiently cold ($T_0/T_p \ge 3$, red and purple dots), the heat current exhibits superextensive scaling with $N$.

Despite this, a discernible gap remains between the heat current $\dot{Q}_h$ calculated in the presence of the parasitic bath and the ideal one without the parasitic bath calculated assuming Eq. (27) with $\dot{J} = N/2$ (represented by black triangles). This gap exists because the parasitic bath redistributes the population across various subspaces characterized by different values of $J$, whereas in Eq. (27) we assume that a single value of $\dot{J} = N/2$ is occupied.

To further interpret our results, in Figs. 4(b) and 4(d), we show the dependence of the total angular momentum $\langle \hat{J}^2 \rangle$, divided by $\sim N^2$, on the temperature of the parasitic bath. In both panels, we see that a hot parasitic bath (blue, orange, and green dots) corresponds to low values of $\langle \hat{J}^2 \rangle$, which is consistent with low values of the heat currents shown in the corresponding Figs 4(a) and 4(c). This can be interpreted in the following way. Let us consider the high-temperature regime of the main baths, which is the one that can exhibit the superlinear scaling of the heat current (see Sec. II C). In the limit of infinitely hot parasitic bath, all states would be equally populated. As $N$ increases, the states with $J \sim 0$ exponentially outnumber the other states, leading to the observed decrease of $\hat{J}^2$. This, in turn, implies a vanishing $\langle \hat{J}_z \rangle$, which, using Eq. (24), implies a vanishing heat current even below the independent qubits case—as observed, e.g., in the orange and blue dots in Figs. 4(a) and 4(c).

Therefore, to achieve superextensive scaling, a cold parasitic bath $T_p \ll T_0$ is required to ensure that the subspace with maximum collective spin $J \approx N/2$ is occupied, corresponding to $\langle \hat{J}^2 \rangle \approx N/2(N/2 + 1)$ in Figs. 4(b) and 4(d). However, this requirement appears to conflict with the need for a high effective temperature $T^*$ to have the superlinear scaling emerging from Eq. (30). Nevertheless, given that the parasitic bath is weakly coupled to the system, it does not significantly influence the system's effective temperature $T^*$, which is dominated by $T_h$ and $T_c$. As a result, a cold parasitic
<figure><img src="image_7.png" /><figcaption>FIG. 4. (a),(c) The heat current, divided by $N$, as a function of the number of qubits, $N$, for various values of the ratio $T_0/T_p$. The independent case $\dot{Q}_h^{\text{ind}}$, given by Eq. (26), is shown as a black dashed line. The ideal case without the parasitic bath, and calculated assuming Eq. (27) with $\bar{J} = N/2$, is shown as black triangles. (b),(d) The total angular momentum $\dot{J}^2$ [normalized by $N/2(N/2+1)$] as a function of $N$ for various values of the ratio $T_0/T_p$. (a),(b) $T_h = T_0$; (c),(d) $T_h = 3T_0$. Other parameters are $T_c = T_0/3$, $\gamma_h = \gamma_c = \omega_0/10$, $\gamma_p = 0.001\omega_0$.</figcaption></figure>

bath $T_p \ll T_0$ does not necessarily force the system toward its ground state, allowing for a finite window where superex- tensive scaling of the heat current can be observed when $T^* \gg T_0$, as shown in Fig. 4(b).

IV. EXPERIMENTAL PROPOSAL

In this section, we discuss an experimental proposal to observe the enhanced heat current in a circuit-QED platform. In Fig. 5(a), we plot the dependence of the heat current $\dot{Q}_h$ on the number of qubits for different temperatures of the parasitic bath, as in Fig. 4, but using physical parameters— reported in Table I—that have been used to describe a previous experiment [28]. This plot demonstrates that the collective enhancement of the heat current is quite pronounced when using parameters that are feasible in experimental settings, as long as the temperature of the parasitic bath is roughly below 150 mK.

Furthermore, we propose a minimal single-device experi- ment to demonstrate the collective advantage in a setup with

$N=2$ transmon qubits. Ideally, we would like to compare the collective current $\dot{Q}_h$ with the independent case $\dot{Q}_h^{\text{ind}}$ in Eq. (26), where each qubit is independently coupled with its bath. However, working with two different devices introduces variability in Hamiltonian parameters due to the fabrication differences, which makes a direct comparison difficult. To circumvent this problem, we devise a protocol where we exploit the possibility of controlling the frequency of each transmon by adjusting the external magnetic flux threading through it. By controlling the detuning between the qubits, we can switch between the collective and the independent coupling scheme on the same device.

Let us assume that the two transmon qubits have tunable frequencies, $\omega_0^{(1)}$ and $\omega_0^{(2)}$. When both of these frequencies are resonant with the LC resonator frequency ($\omega_0^{(k)} = \omega_{LC}$ for $k=1, 2$), the transmons are collectively coupled to the same modes of the baths. Under this condition, the heat current is the collective one $\dot{Q}_h$ determined by Eq. (22) and discussed in previous sections. However, if we detune, e.g., the second qubit setting $\omega_0^{(2)} = \omega_{LC}/2$, while keeping the first qubit on
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5. (a) The heat current $\dot{Q}_h$ (divided by $N$) and expressed in femtoWatt (fW) as a function of $N$, for different values of $T_p$ in milliKelvin (mK). Here, we choose $\omega_0/(2\pi) = 4.0$ GHz, $\gamma_h = \gamma_c = 0.5$ GHz. (b) The heat current $\dot{Q}_h$ in femtoWatt (fW) as a function of $T_p$ in milliKelvin (mK) for $N = 2$ qubits. In this panel, the collective heat current $\dot{Q}_h$ is depicted by a blue continuous line, the heat current in the independent scenario $\dot{Q}_h^{\text{ind}}$ is represented by a red dashed line, and the heat current in the detuned scenario is indicated by a green dotted line. For the detuned scenario, the parameters are set as $\omega_0^{(1)}/(2\pi) = 4.0$ GHz and $\omega_0^{(2)}/(2\pi) = 2.0$ GHz, with $\gamma_1 = 0.5$ GHz. The rate $\gamma_2$ is calculated in accordance with Eq. (D5), assuming $Q_i = 20$. In both panels, additional parameters are $\nu_{T_h} = 8.0$ GHz, $\nu_{T_c} = 1.0$ GHz (corresponding to bath temperatures $T_h = 384.0$ mK and $T_c = 96.0$ mK), and a parasitic bath coupling rate of $\gamma_p = 0.01$ GHz.</figcaption></figure>

resonance, $\omega_0^{(1)} = \omega_{LC}$, we can effectively “disconnect” the second qubit. We label the heat current in this configuration as $\dot{Q}_h^{\text{det}}$. In this detuned scenario, the qubits interact with independent modes of the baths and behave independently. Additionally, due to the significant detuning of the second qubit, it negligibly contributes to the heat current. Therefore, only the first qubit significantly contributes to the heat current, and hence we expect that $\dot{Q}_h^{\text{det}} \approx \dot{Q}_h^{\text{ind}}/2$. In Appendix D, we detail the calculation of the heat current $\dot{Q}_h^{\text{det}}$ in the detuned scenario. This derivation allows us to quantitatively support the previous analysis. In Fig. 5(b), we show the heat current

<table><caption>TABLE I. Physical parameters compatible with the experiment in Ref. [28]. In the first line, we reported the frequencies $\nu_T$ associated with the relative temperatures $T$ as $\nu_T \equiv k_B T/h$.</caption><thead><tr><th></th><th>ω<sub>0</sub>/(2π)</th><th>T<sub>c</sub></th><th>T<sub>h</sub></th><th>T<sub>p</sub></th><th>γ<sub>c</sub></th><th>γ<sub>h</sub></th><th>γ<sub>p</sub></th></tr></thead><tbody><tr><td>GHz</td><td>4.0</td><td>2.0</td><td>8.0</td><td>1.04–10.4</td><td>0.5</td><td>0.5</td><td>0.01</td></tr><tr><td>mK</td><td></td><td>96.0</td><td>384.0</td><td>50–500</td><td></td><td></td><td></td></tr></tbody></table>

as a function of the temperature of the parasitic bath for the minimal single-device system comprising $N = 2$ transmon qubits. The continuous blue line represents the collective heat current $\dot{Q}_h$, while the red dashed line indicates the heat current in the independent case $\dot{Q}_h^{\text{ind}}$, calculated as per Eq. (26). The green dotted line depicts twice the heat current in the scenario where one qubit is in resonance and the other is detuned, namely, $2\dot{Q}_h^{\text{det}}$, and determined by Eq. (D4). As we can see in the figure, the detuned case can be used to accurately estimate the independent scenario as $2\dot{Q}_h^{\text{det}} \approx \dot{Q}_h^{\text{ind}}$, and thus to experimentally validate and quantify the heat current enhancement in the collective case.

Significantly, Fig. 5(b) reveals that at a parasitic bath temperature of $T_p = 50$ mK, the collective $\dot{Q}_h$ exhibits an enhancement of approximately 13% compared to the independent scenario. In contrast, at $T_p = 450$ mK, the collective heat current shows a reduction of about 4% relative to the independent case. Concerning the experimental implementation, fabricating a multiqubit circuit is not significantly more demanding than making a single qubit. We believe that the main experimental challenge in realizing such a minimal single-device experiment is ensuring that the two transmons have the same frequency $\omega_0$ and are coupled to the same bath with the same intensity to guarantee the conservation of $\mathcal{J}^2$. Small imperfections in the physical implementation would inevitably result in breaking the conservation of $\mathcal{J}^2$. However, heuristically, these effects would be similar to the presence of the parasitic baths that we studied earlier.

V. DISCUSSION AND CONCLUSIONS

In this work, we have illustrated how the collective interactions between $N$ qubits and two thermal baths can enhance the heat transport across the device, compared to a scenario where each qubit is independently coupled with the baths.

The device can be implemented within current quantum technologies, and we focus on its circuit-QED implementation. After revealing the physical origin of the collective enhancement of the heat current, we analyze its resilience to a third, parasitic thermal bath that acts locally on each qubit, modeling unavoidable single-qubit noise. We then assess the experimental feasibility of our proposal by choosing system parameters that have been measured in a previous experiment [28], and we propose a minimal experimental device based on two superconducting qubits. Our findings indicate that this collective advantage is not only robust to local noise, but also observable under experimental conditions.

It is important to note that the circuit depicted in Fig. 2 is a simplified model. In real-world scenarios, additional factors may come into play. For example, unintended capacitive
coupling between qubits is likely to occur in practical im-
plementations, leading to effective dipole-dipole interactions
between qubits. However, this coupling would primarily affect
the conservation of $\hat{J}^2$, leading to effects qualitatively analo-
gous to those introduced by the parasitic bath.

Another important limitation of our predictions stems from
the use of a Markovian Lindblad master equation, which
is well justified only for small transition rates compared to
the qubit frequency $\omega_0$. This is primarily required by the
secular approximation and the weak-coupling approximation
(commonly referred to as the Born-Markov approximation),
as detailed in the derivation of our master equation in Ap-
pendix B. Given that we chose relatively small rates, i.e.,
approximately 10% of the bare frequency $\omega_0$, the use of the
Lindblad master equation is reasonably justified. Additionally,
it has been shown that nonsecular terms tend to be more
relevant during transient dynamics, while their influence is
significantly reduced at the steady state [85]. Nonetheless, it
would be interesting to explore the influence of nonsecular
terms and non-Markovian effects that emerge at strong cou-
pling and can play a significant role in an experiment [86].

In the future, it will be interesting to investigate similar
collective effects when time-dependent driving is present, al-
lowing the device to also function as a heat engine [87,88] or
as a refrigerator [89,90], depending on the specific physical
parameters.

ACKNOWLEDGMENTS

Numerical work has been performed by using the PYTHON
toolbox QuTiP2 [91] and the PIQS library [92]. We wish
to thank R. Fazio, F. Campaioli, D. De Bernardis, L. Gia-
comelli, M. Brunelli, D. Ferraro, and F. Minganti for useful
discussions. G.M.A. and M.S. acknowledge funding from
the European Research Council (ERC) under the Euro-
pean Union's Horizon 2020 research and innovation program
(Grant Agreement No. 101002955 – CONQUER). P.A.E.
gratefully acknowledges funding by the Berlin Mathematics
Center MATH+ (Grant No. AA2-18). F.N. gratefully ac-
knowledges funding by the BMBF (Berlin Institute for the
Foundations of Learning and Data – BIFOLD), the European
Research Commission (Grant No. ERC CoG 772230) and the
Berlin Mathematics Center MATH+ (Grants No. AA1-6 and
No. AA2-8). This work was supported by Research Council
of Finland Grant No. 312057 (Centre of Excellence program).

APPENDIX A: DERIVATION OF THE HAMILTONIAN (1)

In this Appendix, we provided a detailed derivation of the
Hamiltonian in Eq. (1). We begin by revisiting the quantiza-
tion processes for two key components of the circuits drawn
in Fig. 2(b): the LC circuit and the transmon qubit. The
derivation of these Hamiltonians is detailed in Appendices A 1
and A 2, respectively. Building upon these concepts, we then
proceed to derive the Hamiltonian in Eq. (1) from the circuit
drawn in Fig. 2.

1. Quantization of the LC superconducting circuit

The Hamiltonian for an *LC* circuit, where *L* is the induc-
tance and *C* is the capacitance, is given by

$$
\hat{H}_{LC} = \frac{\hat{Q}^2}{2C} + \frac{\hat{\phi}^2}{2L}, \qquad (A1)
$$

where $\hat{Q}$ denotes the charge on the capacitor and $\hat{\phi}$ represents the magnetic flux through the inductor. These two variables fulfill the canonical commutation relation $[\hat{Q}, \hat{\phi}] = i\hbar$. Since Eq. (A1) is a quadratic Hamiltonian, it is useful to express $\hat{Q}$ and $\hat{\phi}$ in terms of ladder operators $\hat{a}, \hat{a}^{\dagger}$ as

$$
\hat{Q} = \sqrt{\frac{\hbar}{2Z_{LC}}} (\hat{a} + \hat{a}^{\dagger}), \qquad (A2)
$$

$$
\hat{\phi} = i \sqrt{\frac{\hbar Z_{LC}}{2}} (\hat{a}^{\dagger} - \hat{a}), \qquad (A3)
$$

where ladder operator satisfy $[\hat{a}, \hat{a}^{\dagger}] = 1$ and $Z_{LC} = \sqrt{L/C}$ is the impedance associated with the LC circuit. The Hamiltonian can be thus diagonalized in the form

$$
\hat{H}_{LC} = \hbar\omega_{LC}(\hat{a}^{\dagger}\hat{a} + \frac{1}{2}), \qquad (A4)
$$

where the resonant frequency $\omega_{LC}$ of the LC circuit reads
$\omega_{LC} = 1/\sqrt{LC}$.

2. Quantization for a transmon qubit

A transmon qubit can be modeled as a nonlinear circuit
composed of a shunt capacitance and a Josephson junction
in parallel. The Hamiltonian for the transmon qubit can be
written as

$$
\hat{H}_{\mathrm{T}}=\frac{1}{2 C_{\mathrm{T}}} \hat{Q}_{\mathrm{T}}^{2}-E_{\mathrm{J}} \cos \left(\frac{2 \pi \hat{\phi}_{\mathrm{T}}}{\phi_{0}}\right),
\qquad(\mathrm{A} 5)
$$

where $C_T$ is the transmon capacitor, $\phi_0 = h/(2e)$, and $E_J$ represents the Josephson energy. The variables $\hat{Q}_T$ and $\hat{\phi}_T$ are quantum operators representing the charge of excess Cooper pairs and the phase difference across the Josephson junction, respectively. Again, they satisfy the canonical commutation relation,

$$
[\hat{Q}_T, \hat{\phi}_T] = i\hbar. \qquad (A6)
$$

Transmon qubits operate in the regime $E_J \gg E_C$—where $E_C = e^2/(2C)$ denotes the charging energy—to minimize sensitivity to charge noise. When the potential of the transmon is examined in this regime, it can be approximated as a quasi-harmonic potential. Linearizing around the minimum of this potential, the cosine term in the Hamiltonian can be expanded, leading to an almost harmonic behavior with small anharmonic corrections.

Neglecting the nonlinearity stemming from the Josephson
energy, Eq. (A5) can be approximated as

$$
\hat{H}'_{T} = \frac{1}{2C_{T}} \hat{Q}_{T}^{2} + \frac{1}{2L_{J}} \hat{\phi}_{T}^{2}, \qquad (A7)
$$

where $L_J = (\hbar/2e)^2(1/E_J)$.
This is analogous to the Hamiltonian of a quantum harmonic oscillator in Eq. (A1). Hence, $\hat{Q}_T$ and $\hat{\phi}_T$ can be expressed in terms of ladder operators $\hat{b}$, $\hat{b}^\dagger$ (fulfilling the commutation relation $[\hat{b}, \hat{b}^\dagger] = 1$) as

$$ \hat{Q}_T = \sqrt{\frac{\hbar}{2Z_T}} (\hat{b} + \hat{b}^\dagger), \quad (A8) $$

$$ \hat{\phi}_T = i \sqrt{\frac{\hbar Z_T}{2}} (\hat{b}^\dagger - \hat{b}), \quad (A9) $$

where $Z_T = \sqrt{L_J/C_T}$ is the impedance associated with the equivalent LC circuit. Substituting Eqs. (A8) and (A9) into the linearized Hamiltonian $\hat{H}'_T$, we have

$$ \hat{H}'_T = E'_T \hat{b}^\dagger \hat{b}, \quad (A10) $$

where $E'_T = \sqrt{8E_J E_C}$ is the effective frequency of the transmon. This is the harmonic spectrum of the transmon, at the lowest order in $E_C/E_J$. However, the presence of anharmonic terms is crucial for the transmon to function as a qubit. Specifically, the anharmonicity ensures that the energy levels of the transmon are not evenly spaced, allowing for selective addressing of specific energy transitions.

The quartic anharmonic term in the potential, which arises from the Taylor expansion of the cosine function in Eq. (A5), can be expressed as

$$ \delta \hat{H}_T = - \frac{E_J}{4!} \left( \frac{2\pi \hat{\phi}_T}{\phi_0} \right)^4, \quad (A11) $$

where the transmon Hamiltonian is approximated by $\hat{H}_T \approx \hat{H}'_T + \delta\hat{H}_T$. Given the relationship between $\hat{\phi}$ and $\hat{b}$ and $\hat{b}^\dagger$ in Eq. (A9), the quartic term can be reexpressed in terms of these ladder operators. Expanding $\hat{\phi}^4$ and neglecting off-diagonal terms that do not conserve the number of excitations, the perturbation $\delta\hat{H}_T$ becomes

$$ \delta \hat{H}_T = -\frac{E_C}{2} (\hat{b}^\dagger \hat{b}^\dagger \hat{b} \hat{b} + 2 \hat{b}^\dagger \hat{b}). \quad (A12) $$

The term $\hat{b}^\dagger \hat{b}^\dagger \hat{b} \hat{b}$ in the potential introduces an anharmonicity, leading to unevenly spaced energy levels. This anharmonicity—proportional to $E_C$—is crucial for the operation of the transmon as a qubit, allowing for selective addressing of its states. The term $\hat{b}^\dagger \hat{b}$ results in a small shift in the effective energy $E_T$ of the transmon,

$$ E_T \approx \sqrt{8E_J E_C} - E_C. \quad (A13) $$

This equation incorporates the anharmonic corrections from the Josephson potential.

Due to the introduced anharmonicity, our focus can be limited to the subspace spanned by the unexcited state $|0\rangle_T$ (where $\hat{b}|0\rangle_T = 0$) and the excited state $|1\rangle_T \equiv \hat{b}^\dagger|0\rangle_T$. In this subspace, we can truncate the ladder operator to $\hat{b} \approx \hat{\sigma}_-$, $\hat{b}^\dagger \approx \hat{\sigma}_+$, and subsequently $\hat{b}^\dagger \hat{b} \approx (\hat{\sigma}_z + 1)/2$. Hence, the transmon's observables can be expressed in the truncated

basis as

$$ \hat{H}_T = \frac{E_T}{2} (\hat{\sigma}_z + 1), \quad (A14) $$

$$ \hat{Q}_T = \sqrt{\frac{\hbar}{2Z_T}} \hat{\sigma}_x, \quad (A15) $$

$$ \hat{\phi}_T = -\sqrt{\frac{\hbar Z_T}{2}} \hat{\sigma}_y. \quad (A16) $$

### 3. Derivation of the Hamiltonian of the system

Here we derive the Hamiltonian of the system introduced in the main text. The lumped-element circuit diagram, including capacitances, inductances, resistances, and the various variables, is depicted in Fig. 2. This circuit corresponds to the following Lagrangian for the transmons:

$$ \mathcal{L} = \frac{1}{2} \sum_{j=1}^{N} \left[ C_T \dot{\phi}_j^2 - \frac{1}{L_J} \phi_j^2 + \sum_{i=\text{h,c}} C_c^{(i)} (\dot{\phi}_j - V_i)^2 \right], \quad (A17) $$

where $V_i$ is the voltage drop occurring in the resistance $R_i$ and $C_c^{(i)}$ corresponds to the capacitive coupling between the transmons and the RLC circuits. In the following, the capacitive couplings are assumed to be small with respect to the transmon capacities, $C_c^{(i)} \ll C_T$. Under this approximation, in the following, we will trace out the variables $V_i$ and obtain a Markovian master equation for the transmons. The Lagrangian gives rise to canonical momenta $Q_j$ (the charges), which are the derivatives of the Lagrangian with respect to $\dot{\phi}_j$, $Q_j \equiv \partial \mathcal{L} / \partial \dot{\phi}_j$. In terms of $\dot{\phi}_j$, the charges $Q_j$ can be expressed as

$$ Q_j = \left( C_T + \sum_{i=\text{h,c}} C_c^{(i)} \right) \dot{\phi}_j - \sum_{i=\text{h,c}} C_c^{(i)} V_i. \quad (A18) $$

The previous equation can be inverted to obtain the voltages $\dot{\phi}_j$. This relation, at the first order in $C_c^{(i)}/C_T$, reads

$$ \dot{\phi}_j = Q_j \left( \frac{1}{C_T} - \frac{\sum_{i=\text{h,c}} C_c^{(i)}}{C_T^2} \right) + \sum_{i=\text{h,c}} C_c^{(i)} V_i. \quad (A19) $$

The Hamiltonian of the system can thus be obtained as the Legendre transform of the total Lagrangian $\mathcal{L}$ in Eq. (A17) as $H \equiv \sum_{j=1}^N Q_j \dot{\phi}_j - \mathcal{L}$. At this stage, it is now possible to quantize the system Hamiltonian by promoting classical variables to quantum ones as follows:

$$ \begin{aligned} \dot{\phi}_j &\rightarrow \hat{\phi}_j, \\ Q_j &\rightarrow \hat{Q}_j, \\ V_i &\rightarrow \hat{V}_i. \end{aligned} \quad (A20) $$

The total quantized Hamiltonian can be expressed at the first order $C_c^{(i)}/C_T$ as

$$ \hat{H} = \sum_{j=1}^{N} \frac{1}{2} \left[ \frac{1}{C_T} \left( 1 - \frac{\sum_{i=\text{h,c}} C_c^{(i)}}{C_T} \right) \hat{Q}_j^2 + \frac{1}{L_J} \hat{\phi}_j^2 \right] \\ - \sum_{j=1}^{N} \sum_{i=\text{h,c}} \frac{C_c^{(i)}}{C_T} \hat{Q}_j \hat{V}_i, \quad (A21) $$
where we neglected terms $\hat{V}_i^2$ that do not act on the transmons. By truncating this Hamiltonian on the transmon qubit Hilbert space, as in Eq. (A14), we have

$$ \hat{H} \approx \sum_{j=1}^{N} \frac{E_T}{2} (\hat{\sigma}_z^{(j)} + 1) - \sum_{j=1}^{N} \sum_{i=h,c} \hbar G_i \frac{\hat{\sigma}_x^{(j)}}{2} \hat{V}_i, \quad (A22) $$

where coupling between the qubits and the voltage $\hat{V}_i$ is denoted by $G_i$, which is defined as $G_i = 2(C_c^{(i)}/C_T)(2/\hbar Z_T)^{1/2}$. In the main text, we have used the notation $\hbar\omega_0 = E_T$ to represent the qubit energy. With these definitions and notations in place, we have completed the derivation of the Hamiltonian of the system, given in Eq. (6).

## APPENDIX B: DERIVATION OF LINDBLAD MASTER EQUATION

In this Appendix, we derive the Lindblad master equation for the collection of transmons by tracing out the RLC degrees of freedom. The complete system is described by the total Hamiltonian $\hat{H}_{tot}$,

$$ \hat{H}_{tot} = \hat{H} + \hat{H}_{RLC}. \quad (B1) $$

Here, $\hat{H}_{RLC}$ represents the Hamiltonian for RLC circuits, while $\hat{H}$ is the transmon's Hamiltonian defined in Eq. (6). As $\hat{H}_{RLC}$ contains a resistive part, writing an explicit expression involves the coupling with a collection of infinite LC elements, as discussed, for example, in Ref. [69]. As we will show later, here we do not need to write the full RLC Hamiltonian, $\hat{H}_{RLC}$. Instead, thanks to the fluctuation-dissipation theorem, it will be sufficient to know the classical equation of motion of the RLC circuit. We further define the interaction Hamiltonian $\hat{H}_{int}$ as

$$ \hat{H}_{int} = \hat{H}_{tot} - \hat{H} - \hat{H}_{RLC} \quad (B2) $$

$$ = - \sum_{i=h,c} \hbar G_i \hat{J}_x \hat{V}_i. \quad (B3) $$

We employ the Born-Markov approximation, which assumes that the qubits are weakly coupled to the modes of the baths, to trace out the field and obtain a master equation governing the evolution of the qubits' density matrix $\rho(t)$. The calculations are facilitated by working in the interaction picture. For any generic operator $\hat{O}$, its corresponding form in the interaction picture is $\hat{O}^1(t) = U_0^\dagger(t) \hat{O} U_0(t)$, where $U_0(t) = \exp[-i\hat{H}_{bare}t]$ and $\hat{H}_{bare} = \hat{H}_0 + \hat{H}_r$. Meanwhile, the density matrix in the interaction picture evolves according to $\hat{\rho}^1(t) = U_0(t)\hat{\rho}(t)U_0^\dagger(t)$.

The evolution of the total density matrix in the interaction picture, $\hat{\rho}_{tot}^1(t)$, is dictated by the following master equation [70]:

$$ \partial_t \hat{\rho}_{tot}^1(t) = -\frac{1}{\hbar^2} \int_0^\infty d\tau \{ \hat{H}_{int}^1(t), [\hat{H}_{int}^1(t-\tau), \hat{\rho}_{tot}^1(t)] \}. \quad (B4) $$

In reaching this equation, we apply the Markov approximation, assuming that the environment's correlation time $\tau_{RLC}$ is much shorter than the timescale of the qubit's dynamics, $\omega_0^{-1}$, namely, $\omega_0 \tau_{RLC} \ll 1$. This allows us to approximate $\hat{\rho}_{tot}^1(t-\tau) \approx \hat{\rho}_{tot}^1(t)$ under the time integral and to extend the integration limit to $\tau \to \infty$.

Upon tracing out the field subsystem, the master equation for the matter density matrix $\hat{\rho}^1(t)$ is given by

$$ \partial_t \hat{\rho}^1(t) = -\frac{1}{\hbar^2} \int_0^\infty d\tau \{ \text{tr}_r [\hat{H}_{\text{int}}^1(t) \hat{H}_{\text{int}}^1(t-\tau) \hat{\rho}_r^1(t)] \hat{\rho}^1(t) \\ - \text{tr}_r [\hat{H}_{\text{int}}^1(t) \hat{\rho}_r^1(t) \hat{\rho}^1(t) \hat{H}_{\text{int}}^1(t-\tau)] + \text{H.c.} \}, \quad (B5) $$

where $\text{tr}_r[\cdot]$ denotes the trace on the RLC degrees of freedom. Here, we employed the Born approximation, which assumes that the RLC circuits serve as large memoryless reservoirs. This lets us approximate $\hat{\rho}_{tot}^1(t) \approx \hat{\rho}^1(t)\hat{\rho}_r^1(t)$.

In the interaction picture, the interaction Hamiltonian defined in Eq. (B2) reads

$$ \hat{H}_{\text{int}}^{1}(t) = - \sum_{i=h,c} \hbar G_i [\hat{J}_{-}^{1}(t) + \hat{J}_{+}^{1}(t)] \hat{V}_{i}^{1}(t). \quad (B6) $$

Here, $\hat{V}_i^1(t)$ are the voltage operators associated with the RLC circuits, and $\hat{J}_+^1(t)$ and $\hat{J}_-^1(t)$ are the qubits collective raising and lowering operators. The qubit raising and lowering operators in the interaction picture are given by

$$ \hat{J}_{+}^{1}(t) = e^{i\omega_{0}t} \hat{J}_{+}, \quad (B7) $$

$$ \hat{J}_{-}^{1}(t) = e^{-i\omega_{0}t} \hat{J}_{-}. \quad (B8) $$

The phase factors arise due to the transformation to the interaction picture, where $\omega_0$ is the bare frequency of the qubits. By substituting ladder operators in the interaction picture into the master equation [Eq. (B5)] and performing the secular approximation, we select operators that induce transitions between eigenstates. The secular approximation involves neglecting rapidly oscillating terms in the interaction picture, which is valid when $\omega_0$ is large compared to the system-bath couplings. In this regime, we obtain

$$ \partial_t \hat{\rho}(t) = -\frac{i}{\hbar} [\hat{H}_0, \hat{\rho}(t)] - \sum_{i,j} G_i G_j \{ S_{\hat{\nu}_i, \hat{\nu}_j}(-\omega_0) [\hat{J}_+ \hat{J}_- \hat{\rho}(t) \\ - \hat{J}_+ \hat{\rho}(t) \hat{J}_-] + S_{\hat{\nu}_i, \hat{\nu}_j}(\omega_0) [\hat{J}_- \hat{J}_+ \hat{\rho}(t) \\ - \hat{J}_- \hat{\rho}(t) \hat{J}_+] + \text{H.c.} \}, \quad (B9) $$

where $S_{\hat{\nu}_i, \hat{\nu}_j}(\omega)$ is the voltage dynamical structure factor of the RLC circuits,

$$ S_{\hat{\nu}_i, \hat{\nu}_j}(\omega) = \int_0^\infty d\tau e^{-i\omega\tau} \langle \hat{V}_i^1(\tau) \hat{V}_j \rangle_r. \quad (B10) $$

The voltage dynamical structure factor can be calculated by employing the fluctuation-dissipation theorem [69],

$$ S_{\hat{\nu}_i, \hat{\nu}_j}(\omega) = 2\hbar\omega\text{Re}[Z_{\text{tot},i}(\omega)]\mathcal{N}_i(\omega)\delta_{i,j}, \quad (B11) $$

where $Z_{\text{tot},i}(\omega)$ is the total impedance of the RLC circuit and $\mathcal{N}_i(\omega) = [n_i(\omega) + 1]\Theta(\omega) - n_i(-\omega)\Theta(-\omega)$ is a thermal occupation factor. The total impedance of the RLC circuits in Fig. 2 can be obtained by considering that the LC branch is in parallel with the resistor branch, hence,

$$ Z_{\text{tot},i}(\omega) = \left( \frac{1}{Z_{\text{LC}}(\omega)} + \frac{1}{Z_{\text{R}_i}} \right)^{-1}, \quad (B12) $$
where $Z_R = R_i$ is the impedance of the ith resistor and $Z_{LC}(\omega)$ is the impedance of the LC circuit,

$$Z_{LC}(\omega) = \frac{i}{C} \frac{\omega}{\omega^2 - \omega_{LC}^2}. \qquad (B13)$$

By placing the explicit form of $Z_{LC}(\omega)$ and $Z_R$ in Eq. (B12), we obtain

$$\mathrm{Re}[Z_{\mathrm{tot},i}(\omega)] = \frac{\omega^2 R_i}{\omega^2 + R_i^2 C^2 (\omega^2 - \omega_{LC}^2)^2}. \qquad (B14)$$

This expression can be further simplified by introducing the quality factor $Q_i$ as the ratio of the frequency of the RLC circuit divided by its line width, $\kappa_i$, $Q_i = \omega_{LC}/\kappa_i$. The linewidth of the RLC circuit resonance is given by $\kappa_i = R_i$. Hence, we have

$$\mathrm{Re}[Z_{\mathrm{tot},i}(\omega)] = \frac{R_i}{1 + Q_i^2 \left( \frac{\omega}{\omega_{LC}} - \frac{\omega_{LC}}{\omega} \right)^2}. \qquad (B15)$$

By plugging this equation into Eq. (B11), we obtain an explicit expression for the dynamical structure factor,

$$S_{\hat{\nu}_i, \hat{\nu}_j}(\omega) = \left[ \frac{2\hbar\omega R_i N_i(\omega)}{1 + Q_i^2 \left( \frac{\omega}{\omega_{LC}} - \frac{\omega_{LC}}{\omega} \right)^2} \right] \delta_{ij}. \qquad (B16)$$

By plugging the explicit form of the dynamical structure factor in Eq. (B9), we can finally arrive at a master equation that can be written in the form of Eq. (9),

$$\partial_t \hat{\rho}(t) = -\frac{i}{\hbar} [\hat{H}_0, \hat{\rho}(t)] + \mathcal{D}[\hat{\rho}(t)]. \qquad (B17)$$

Here, the dissipator $\mathcal{D}[\rho(t)]$ is given by

$$\mathcal{D}[\rho(t)] = \sum_{i=\mathrm{h},\mathrm{c}} \gamma_i \left\{ (1+n_i) \left[ \hat{J}_{-} \hat{\rho} \hat{J}_{+} - \frac{1}{2} (\hat{J}_{+} \hat{J}_{-}, \hat{\rho}) \right] + n_i \left[ \hat{J}_{+} \hat{\rho} \hat{J}_{-} - \frac{1}{2} (\hat{J}_{-} \hat{J}_{+}, \hat{\rho}) \right] \right\}, \qquad (B18)$$

where the rates read

$$\gamma_i = \frac{4\hbar\omega_0 R_i G_i^2}{1 + Q_i^2 \left( \frac{\omega_0}{\omega_{LC}} - \frac{\omega_{LC}}{\omega_0} \right)^2}. \qquad (B19)$$

## APPENDIX C: DETAILS OF THE HEAT CURRENT

In this Appendix, we provide detailed mathematical calculations related to the heat current. Appendix C 1 is dedicated to deriving an explicit expression for the heat current, specifically, Eq. (22). Following this, in Appendix C 2, we discuss various limits within which a simpler analytical expression for the heat current can be obtained.

### 1. Derivation of Eq. (22)

Starting from Eq. (21), we can express the heat current contribution from a specific bath. We will focus on deriving the heat current from the cold bath, $\dot{Q}_h$.

The heat current from the cold bath is given by

$$\dot{Q}_i = \hbar\omega_0 \mathrm{Tr}[\hat{J}_z \mathcal{D}_i(\hat{\rho})].$$

By inserting the expression for $\mathcal{D}_i(\hat{\rho})$, using the cyclic property of the trace, and introducing the notation $\langle \hat{x} \rangle \equiv \mathrm{Tr}[\hat{\rho}\hat{x}]$, we can rewrite this expression as

$$\dot{Q}_i = \hbar\omega_0\gamma_i[(1+n_i)(\langle\hat{J}_+\hat{J}_z\hat{J}_-\rangle - \frac{1}{2}\langle\langle\hat{J}_+\hat{J}_-, \hat{J}_z\rangle\rangle) + n_i(\langle\langle\hat{J}_-\hat{J}_+\rangle - \frac{1}{2}\langle\langle\hat{J}_-\hat{J}_+, \hat{J}_z\rangle\rangle)].$$

Using the fact that $\hat{J}_z$ commutes with $\hat{J}_-\hat{J}_+$, we arrive at

$$\dot{Q}_i = \hbar\omega_0\gamma_i[(1 + n_i)\langle\hat{J}_+ [\hat{J}_z, \hat{J}_-]\rangle + n_i\langle\hat{J}_- [\hat{J}_z, \hat{J}_+]\rangle].$$

Using the commutation relation $[\hat{J}_z, \hat{J}_\pm] = \pm \hat{J}_\pm$, we arrive at a desired final expression for the heat current $\dot{Q}_i$,

$$\dot{Q}_i = \hbar\omega_0\gamma_i[-(1 + n_i)\langle\hat{J}_+\hat{J}_-\rangle + n_i\langle\hat{J}_-\hat{J}_+\rangle],$$

which corresponds to Eq. (22).

### 2. Analytical calculation of the heat current

In this Appendix, we detail calculations of the heat current in Eq. (24). Essentially, one needs to calculate the average value of the collective spin inversion, $\langle -\hat{J}_z \rangle$. Given a generic operator $\hat{x}$, its average can be computed on the steady-state density matrix $\hat{\rho}^{(s)}$ in Eq. (17) as

$$\langle \hat{x} \rangle = \mathrm{Tr}[\hat{\rho}^{(s)} \hat{x}] \qquad (C1)$$

$$= \sum_{J} P_{J} \langle \hat{x} \rangle_{J,T}, \qquad (C2)$$

where $\langle \hat{x} \rangle_{J,T}$ denote averages within a subspace with a given $J$, over the thermal occupations $P(m_J|J) = e^{-m\hbar\omega_0/k_B T}/Z_{J,T}$ depending on $J$ and on the temperature $T$,

$$\langle \hat{x} \rangle_{J,T} = \sum_{m=-J}^{J} \frac{\exp\left(-\frac{m\hbar\omega_0}{k_B T}\right)}{Z_{J,T}} \langle J, m|\hat{x}|J, m \rangle, \qquad (C3)$$

and $Z_{J,T}$ is the partition function,

$$Z_{J,T} = \sum_{m=-J}^{J} \exp\left(-\frac{m\hbar\omega_0}{k_B T}\right). \qquad (C4)$$

Then our goal reduces to calculating the average value $\langle \hat{J}_z \rangle$ in a subspace with a fixed $J$,

$$\langle \hat{J}_z \rangle_{J,T} = \frac{\sum_{m=-J}^{J} m \exp \left( -\frac{m\hbar\omega_0}{k_B T} \right)}{Z_{J,T}}. \qquad (C5)$$

This sum can be performed by noticing that $Z_{J,T}$ is a geometric series. The sum of this geometric series is given by

$$Z_{J,T} = \frac{\exp\left(\frac{J\hbar\omega_0}{k_B T}\right) - \exp\left(-\frac{(J+1)\hbar\omega_0}{k_B T}\right)}{1 - \exp\left(-\frac{\hbar\omega_0}{k_B T}\right)}. \qquad (C6)$$

Given the relation

$$\frac{\partial}{\partial\left(\frac{1}{k_B T}\right)} \exp\left(-\frac{m\hbar\omega_0}{k_B T}\right) = -\hbar\omega_0 m \exp\left(-\frac{m\hbar\omega_0}{k_B T}\right), \qquad (C7)$$

we can rewrite $\langle \hat{J}_z \rangle_{J,T}$ as

$$\langle \hat{J}_z \rangle_{J,T} = -\frac{1}{\hbar\omega_0} \frac{\partial \ln(Z_{J,T})}{\partial\left(\frac{1}{k_B T}\right)}. \qquad (C8)$$

This expression provides a direct way to calculate the expectation value $\langle \hat{J}_z \rangle_{J,T}$ by differentiating the logarithm of the
<figure><img src="image_10.png" /><figcaption>FIG. 6. This figure shows the average $\langle \hat{J}_z \rangle_{J,T}$ in a given subspace with $J = \bar{J}$ fixed as a function of the temperature $T$.</figcaption></figure>

partition function with respect to $1/(k_B T)$. By doing this derivative in Eq. (C8), we obtain an explicit expression for $\langle \hat{J}_z \rangle_{J,T}$,

$$ \langle \hat{J}_z \rangle_{J,T} = -J + \frac{2J + 1}{1 - \exp\left(\frac{(2J+1)\hbar\omega_0}{k_B T}\right)} + \frac{1}{\exp\left(\frac{\hbar\omega_0}{k_B T}\right) - 1}. \quad (C9) $$

Figure 6 depicts the variation of the average $\langle \hat{J}_z \rangle_{J,T}$ as a function of temperature $T$. The graph illustrates a linear behavior for small effective temperatures, where $J\hbar\omega_0 \gg k_B T$, and transitions to a superlinear scaling in the opposite regime, $J\hbar\omega_0 \ll k_B T$. This behavior is indicative of the different thermal regimes influencing the collective angular momentum's $z$ component.

It is insightful to consider the specific limits of Eq. (C9). In the regime of low temperatures, characterized by $\hbar\omega_0 J \gg k_B T$, the average $\langle \hat{J}_z \rangle_{J,T}$ simplifies to

$$ \langle \hat{J}_z \rangle_{J,T} \approx -J. \quad (C10) $$

This expression proves that the bound in Eq. (28) is tight in the low-temperature regime.

On the other hand, in the high-temperature regime, where $\hbar\omega_0 J \ll k_B T$, we obtain

$$ \langle \hat{J}_z \rangle_{J,T} \approx -\left(\frac{\hbar\omega_0}{3k_B T}\right)J(1+J). \quad (C11) $$

This equation is utilized in the main text to derive Eq. (30). Additionally, evaluating the ratio $\langle \hat{J}_z \rangle_{J,T} / (J \langle \hat{J}_z \rangle_{1/2,T})$ in the limit $J \to \infty$ yields

$$ \frac{\langle \hat{J}_z \rangle_{J,T}}{J \langle \hat{J}_z \rangle_{1/2,T}} \approx 2 \coth \left[ \frac{\hbar\omega_0}{2k_B T} \right]. \quad (C12) $$

This expression is applied in the main text to derive Eq. (28).

## APPENDIX D: CALCULATION OF $\hat{Q}_h^{\text{det}}$

In this Appendix, we detail the calculation of the heat current $\hat{Q}_h^{\text{det}}$ in the detuned scenario. In this case, each qubit interacts independently with distinct modes of the bath. The total heat current is thus given by the sum of the independent contributions of each qubit, each one determined by Eq. (22) setting $N=1$, with the qubit's state determined by the master equation in Eq. (33). We assume equal strength in the couplings with both the hot and cold baths. Therefore, in the detuned scenario, the total heat current in Eq. (22) can be expressed as

$$ \hat{Q}_h^{\text{det}} = \sum_{k=1,2} \hbar\omega_0^{(k)} \gamma_k \left[ -\left(1 + n_h^{(k)}\right) \langle \hat{\sigma}_+^{(k)} \hat{\sigma}_-^{(k)} \rangle + n_h^{(k)} \langle \hat{\sigma}_-^{(k)} \hat{\sigma}_+^{(k)} \rangle \right], $$

where $\gamma_k$ represents the decay rate of the $k$th qubit, $\hat{\sigma}_+^{(k)}(\hat{\sigma}_-^{(k)})$ are the creation (destruction) Pauli operators acting on the $k$th qubit, and the thermal occupations $n_i^{(k)}$ of the $i$th bath depend on the bath temperature $T_i$ and the frequency of the $k$th qubit,

$$ n_i^{(k)} = \frac{1}{\exp\left(\frac{\hbar\omega_0^{(k)}}{k_B T_i}\right) - 1}. \quad (D1) $$

At the steady state, the average values $\langle \hat{\sigma}_+^{(k)} \hat{\sigma}_-^{(k)} \rangle$, $\langle \hat{\sigma}_-^{(k)} \hat{\sigma}_+^{(k)} \rangle$ reduce to the thermal populations $p_0^{(k)}$, $p_1^{(k)}$, as

$$ \begin{aligned} \langle \hat{\sigma}_+^{(k)} \hat{\sigma}_-^{(k)} \rangle &= p_1^{(k)}, \\ \langle \hat{\sigma}_-^{(k)} \hat{\sigma}_+^{(k)} \rangle &= p_0^{(k)}, \end{aligned} \quad (D2) $$

where the qubit's populations $p_0^{(k)}$, $p_1^{(k)}$ are thermally occupied as dictated by Eq. (33): $p_0^{(k)} = 1/[1 + \exp(-\frac{\hbar\omega_0^{(k)}}{k_B T_k^{*}})]$ and $p_1^{(k)} = 1 - p_0^{(k)}$.

The effective temperature $T_k^*$ is calculated using Eq. (33), analogously to the effective temperature in the absence of the parasitic bath as given in Eq. (15),

$$ \frac{\gamma_k (n_h^{(k)} + n_c^{(k)}) + \gamma_p n_p^{(k)}}{\gamma_k (n_h^{(k)} + n_c^{(k)} + 2) + \gamma_p (n_p^{(k)} + 1)} = \exp\left(-\frac{\hbar\omega_0^{(k)}}{k_B T_k^*}\right). \quad (D3) $$

Hence, Eq. (D1) can be recast as

$$ \hat{Q}_h^{\text{det}} = \sum_{k=1,2} \hbar\omega_0^{(k)} \gamma_k \left[ -\left(1 + n_h^{(k)}\right) p_1^{(k)} + n_h^{(k)} p_0^{(k)} \right]. \quad (D4) $$

The rate for the first qubit is set at $\gamma_1 = 1$ GHz, as per Table I. Detuning the second qubit from the LC resonance significantly reduces its rate, based on the microscopic expression in Eq. (B19). Hence, for the second qubit's rate $\gamma_2$, we use

$$ \gamma_2 = \frac{\gamma_1}{1 + Q_i^2 (\omega_0 / \omega_{LC} - \omega_{LC} / \omega_0)^2}, \quad (D5) $$

with $Q_i$ being the quality factors of the resonators.

[1] R. Kosloff, Quantum thermodynamics: A dynamical viewpoint, Entropy 15, 2100 (2013).

[2] G. Benenti, G. Casati, K. Saito, and R. S. Whitney, Fundamental aspects of steady-state conversion of heat to work at the nanoscale, Phys. Rep. 694, 1 (2017).

[3] S. Vinjanampathy and J. Anders, Quantum thermodynamics, Contemp. Phys. 57, 545 (2016).

[4] J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, The role of quantum information in thermodynamics—A topical review, J. Phys. A: Math. Theor. 49, 143001 (2016).
[5] F. Binder, L. Correa, C. Gogolin, J. Anders, and G. Adesso, *Thermodynamics in the Quantum Regime Fundamental Aspects and New Directions*, Fundamental Theories of Physics (Springer International Publishing, Cham, 2018).

[6] M. Lostaglio, An introductory review of the resource theory approach to thermodynamics, *Rep. Prog. Phys.* **82**, 114001 (2019).

[7] A. Auffèves, Quantum technologies need a quantum energy initiative, *PRX Quantum* **3**, 020101 (2022).

[8] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, *Rev. Mod. Phys.* **93**, 025005 (2021).

[9] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler et al., Quantum supremacy using a programmable superconducting processor, *Nature (London)* **574**, 505 (2019).

[10] A. Fornieri and F. Giazotto, Towards phase-coherent caloritronics in superconducting circuits, *Nat. Nanotechnol.* **12**, 944 (2017).

[11] O. Maillet, D. Subero, J. T. Peltonen, D. S. Golubev, and J. P. Pekola, Electric field control of radiative heat transfer in a superconducting circuit, *Nat. Commun.* **11**, 4326 (2020).

[12] N. Ligato, F. Paolucci, E. Strambini, and F. Giazotto, Thermal superconducting quantum interference proximity transistor, *Nat. Phys.* **18**, 627 (2022).

[13] A. Gubaydullin, G. Thomas, D. S. Golubev, D. Lvov, J. T. Peltonen, and J. P. Pekola, Photonic heat transport in three terminal superconducting circuit, *Nat. Commun.* **13**, 1552 (2022).

[14] E. Gümüş, D. Majidi, D. Nikolić, P. Raif, B. Karimi, J. T. Peltonen, E. Scheer, J. P. Pekola, H. Courtois, W. Belzig, and C. B. Winkelmann, Calorimetry of a phase slip in a Josephson junction, *Nat. Phys.* **19**, 196 (2023).

[15] M. Campisi, R. Blattmann, S. Kohler, D. Zueco, and P. Hänggi, Employing circuit QED to measure nonequilibrium work fluctuations, *New J. Phys.* **15**, 105028 (2013).

[16] B. Karimi and J. P. Pekola, Otto refrigerator based on a superconducting qubit: Classical and quantum performance, *Phys. Rev. B* **94**, 184503 (2016).

[17] P. P. Hofer, M. Perarnau-Llobet, J. B. Brask, R. Silva, M. Huber, and N. Brunner, Autonomous quantum refrigerator in a circuit QED architecture based on a Josephson junction, *Phys. Rev. B* **94**, 235420 (2016).

[18] P. P. Hofer, J.-R. Souquet, and A. A. Clerk, Quantum heat engine based on photon-assisted Cooper pair tunneling, *Phys. Rev. B* **93**, 041418(R) (2016).

[19] R. Sánchez, Correlation-induced refrigeration with superconducting single-electron transistors, *Appl. Phys. Lett.* **111**, 223103 (2017).

[20] K. Y. Tan, M. Partanen, R. E. Lake, J. Govenius, S. Masuda, and M. Möttönen, Quantum-circuit refrigerator, *Nat. Commun.* **8**, 15189 (2017).

[21] A. U. C. Hardal, N. Aslan, C. M. Wilson, and O. E. Müstecaplioğlu, Quantum heat engine with coupled superconducting resonators, *Phys. Rev. E* **96**, 062120 (2017).

[22] K. Funo, N. Lambert, B. Karimi, J. P. Pekola, Y. Masuyama, and F. Nori, Speeding up a quantum refrigerator via counterdiabatic driving, *Phys. Rev. B* **100**, 035407 (2019).

[23] S. M. Tabatabaei, D. Sánchez, A. L. Yeyati, and R. Sánchez, Nonlocal quantum heat engines made of hybrid superconducting devices, *Phys. Rev. B* **106**, 115419 (2022).

[24] P. A. Erdman and F. Noé, Identifying optimal cycles in quantum thermal machines with reinforcement-learning, *npj Quantum Inf.* **8**, 1 (2022).

[25] A. Guthrie, C. D. Satrya, Y.-C. Chang, P. Menczel, F. Nori, and J. P. Pekola, Cooper-pair box coupled to two resonators: An architecture for a quantum refrigerator, *Phys. Rev. Appl.* **17**, 064022 (2022).

[26] M. A. Aamir, P. J. Suria, J. A. M. Guzmán, C. Castillo-Moreno, J. M. Epstein, N. Y. Halpern, and S. Gasparinetti, Thermally driven quantum refrigerator autonomously resets superconducting qubit, *arXiv:2305.16710*.

[27] G. Manzano and R. López, Quantum-enhanced performance in superconducting Andreev reflection engines, *Phys. Rev. Res.* **5**, 043041 (2023).

[28] A. Ronzani, B. Karimi, J. Senior, Y.-C. Chang, J. T. Peltonen, C. Chen, and J. P. Pekola, Tunable photonic heat transport in a quantum heat valve, *Nat. Phys.* **14**, 991 (2018).

[29] J. Senior, A. Gubaydullin, B. Karimi, J. T. Peltonen, J. Ankerhold, and J. P. Pekola, Heat rectification via a superconducting artificial atom, *Commun. Phys.* **3**, 40 (2020).

[30] J. P. Pekola and B. Karimi, Colloquium: Quantum heat transport in condensed matter systems, *Rev. Mod. Phys.* **93**, 041001 (2021).

[31] R. Fazio, F. W. J. Hekking, and D. E. Khmelnitskii, Anomalous thermal transport in quantum wires, *Phys. Rev. Lett.* **80**, 5611 (1998).

[32] C. L. Kane and M. P. A. Fisher, Thermal transport in a Luttinger liquid, *Phys. Rev. Lett.* **76**, 3192 (1996).

[33] S. Jezouin, F. D. Parmentier, A. Anthore, U. Gennser, A. Cavanna, Y. Jin, and F. Pierre, Quantum limit of heat flow across a single electronic channel, *Science* **342**, 601 (2013).

[34] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Finite-temperature transport in one-dimensional quantum lattice models, *Rev. Mod. Phys.* **93**, 025003 (2021).

[35] A. Larzul and M. Schirò, Energy transport across two interacting quantum baths without quasiparticles, *Phys. Rev. B* **108**, 115120 (2023).

[36] R. H. Dicke, Coherence in spontaneous radiation processes, *Phys. Rev.* **93**, 99 (1954).

[37] J. M. Raimond, M. Brune, and S. Haroche, Manipulating quantum entanglement with atoms and photons in a cavity, *Rev. Mod. Phys.* **73**, 565 (2001).

[38] M. Gross and S. Haroche, Superradiance: An essay on the theory of collective spontaneous emission, *Phys. Rep.* **93**, 301 (1982).

[39] Y. Kaluzny, P. Goy, M. Gross, J. M. Raimond, and S. Haroche, Observation of self-induced rabi oscillations in two-level atoms excited inside a resonant cavity: The ringing regime of superradiance, *Phys. Rev. Lett.* **51**, 1175 (1983).

[40] A. Angerer, K. Streltsov, T. Astner, S. Putz, H. Sumiya, S. Onoda, J. Isoya, W. J. Munro, K. Nemoto, J. Schmiedmayer, and J. Majer, Superradiant emission from colour centers in diamond, *Nat. Phys.* **14**, 1168 (2018).

[41] Z. Wang, H. Li, W. Feng, X. Song, C. Song, W. Liu, Q. Guo, X. Zhang, H. Dong, D. Zheng, H. Wang, and D.-W. Wang, Controlable switching between superradiant and subradiant states in a
10-qubit superconducting circuit, Phys. Rev. Lett. 124, 013601 (2020).

[42] L. A. Correa, M. Mehboudi, G. Adesso, and A. Sanpera, In-
dividual quantum probes for optimal thermometry, Phys. Rev.
Lett. 114, 220405 (2015).

[43] M. Mehboudi, M. R. Jørgensen, S. Seah, J. B. Brask, J.
Kołodyński, and M. Perarnau-Llobet, Fundamental limits in
Bayesian thermometry and attainability via adaptive strategies,
Phys. Rev. Lett. 128, 130502 (2022).

[44] P. Abiuso, P. A. Erdman, M. Ronen, F. Noé, G. Haack, and M. Perarnau-Llobet, Optimal thermometers with spin networks, Quantum Sci. Technol. 9, 035008 (2024).

[45] N. Jaseem, S. Vinjanampathy, and V. Mukherjee, Quadratic enhancement in the reliability of collective quantum engines, Phys. Rev. A 107, L040202 (2023).

[46] M. Campisi and R. Fazio, The power of a critical heat engine,
Nat. Commun. 7, 11895 (2016).

[47] L. Fusco, M. Paternostro, and G. De Chiara, Work extraction and energy storage in the Dicke model, Phys. Rev. E 94, 052122 (2016).

[48] W. Niedenzu and G. Kurizki, Cooperative many-body enhance-
ment of quantum thermal machine power, New J. Phys. 20,
113038 (2018).

[49] G. Manzano, G.-L. Giorgi, R. Fazio, and R. Zambrini, Boosting the performance of small autonomous refrigerators via common environmental effects, New J. Phys. 21, 123026 (2019).

[50] M. T. Mitchison, Quantum thermal absorption machines: Re-
frigerators, engines and clocks, Contemp. Phys. 60, 164 (2019).

[51] P. Abiuso and M. Perarnau-Llobet, Optimal cycles for low-
dissipation heat engines, Phys. Rev. Lett. 124, 110606 (2020).

[52] L. d. S. Souza, G. Manzano, R. Fazio, and F. Iemini, Collec-
tive effects on the performance and stability of quantum heat
engines, Phys. Rev. E 106, 014143 (2022).

[53] Y.-q. Liu, Y.-j. Yang, T.-t. Ma, and C.-s. Yu, Quantum heat valve and entanglement in superconducting *LC* resonators, Appl. Phys. Lett. 123, 144002 (2023).

[54] S. Kamimura, K. Yoshida, Y. Tokura, and Y. Matsuzaki, Univer-
sal scaling bounds on a quantum heat current, Phys. Rev. Lett.
131, 090401 (2023).

[55] M. Carrega, L. Razzoli, P. A. Erdman, F. Cavaliere, G. Benenti, and M. Sassetti, Dissipation-induced collective advantage of a quantum thermal machine, AVS Quantum Sci. 6, 025001 (2024).

[56] H. Tajima and K. Funo, Superconducting-like heat current: Ef-
fective cancellation of current-dissipation trade-off by quantum
coherence, Phys. Rev. Lett. 127, 190604 (2021).

[57] D. Ferraro, M. Campisi, G. M. Andolina, V. Pellegrini, and M. Polini, High-power collective charging of a solid-state quantum battery, Phys. Rev. Lett. 120, 117702 (2018).

[58] G. M. Andolina, M. Keck, A. Mari, M. Campisi, V. Giovannetti, and M. Polini, Extractable work, the role of correlations, and asymptotic freedom in quantum batteries, Phys. Rev. Lett. 122, 047702 (2019).

[59] J. Q. Quach, K. E. McGhee, L. Ganzer, D. M. Rouse, B. W. Lovett, E. M. Gauger, J. Keeling, G. Cerullo, D. G. Lidzey, and T. Virgili, Superabsorption in an organic microcavity: Toward a quantum battery, Sci. Adv. 8, eabk3160 (2022).

[60] F. Campaioli, S. Gherardini, J. Q. Quach, M. Polini, and G. M. Andolina, Colloquium: Quantum batteries, Rev. Mod. Phys. 96, 031001 (2024).

[61] F. Campaioli, F. A. Pollock, F. C. Binder, L. Céleri, J. Goold, S. Vinjanampathy, and K. Modi, Enhancing the charging power of quantum batteries, Phys. Rev. Lett. 118, 150601 (2017).

[62] D. Kolisnyk and G. Schaller, Performance boost of a collective qutrit refrigerator, Phys. Rev. Appl. 19, 034023 (2023).

[63] M. Vogl, G. Schaller, and T. Brandes, Counting statistics of collective photon transmissions, Ann. Phys. 326, 2827 (2011).

[64] A. Crescente, D. Ferraro, M. Carrega, and M. Sassetti, Enhanc-
ing coherent energy transfer between quantum devices via a
mediator, Phys. Rev. Res. 4, 033216 (2022).

[65] A. Rolandi, P. Abiuso, and M. Perarnau-Llobet, Collective ad-
vantages in finite-time thermodynamics, Phys. Rev. Lett. 131,
210401 (2023).

[66] J. Kim, S.-h. Oh, D. Yang, J. Kim, M. Lee, and K. An, A photonic quantum engine driven by superradiance, Nat. Photon. 16, 707 (2022).

[67] J. P. Pekola, Towards quantum thermodynamics in electronic circuits, Nat. Phys. 11, 118 (2015).

[68] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer's guide to superconducting qubits, Appl. Phys. Rev. 6, 021318 (2019).

[69] M. Cattaneo and G. S. Paraoanu, Engineering dissipation with resistive elements in circuit quantum electrodynamics, Adv. Quantum Technol. 4, 2100054 (2021).

[70] H. P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Oxford University Press, Oxford 2002).

[71] V. V. Albert and L. Jiang, Symmetries and conserved quantities in Lindblad master equations, Phys. Rev. A 89, 022118 (2014).

[72] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Boundary time crystals, Phys. Rev. Lett. 121, 035301 (2018).

[73] P. P. Hofer, M. Perarnau-Llobet, L. D. M. Miranda, G. Haack, R. Silva, J. B. Brask, and N. Brunner, Markovian master equations for quantum thermal machines: local versus global approach, New J. Phys. 19, 123037 (2017).

[74] M. Cattaneo, G. L. Giorgi, S. Maniscalco, and R. Zambrini, Local versus global master equation with common and separate baths: Superiority of the global approach in partial secular approximation, New J. Phys. 21, 113045 (2019).

[75] G. D. Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia, and M. Antezza, Reconciliation of quantum local master equations with thermodynamics, New J. Phys. 20, 113024 (2018).

[76] J. J. Sakurai, *Modern Quantum Mechanics* (Addison-Wesley,
Reading, MA, 1994).

[77] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A: Math. Gen. 12, L103 (1979).

[78] Z. Yan, J. Ho, Y.-H. Lu, S. J. Masson, A. Asenjo-Garcia, and D. M. Stamper-Kurn, Superradiant and subradiant cavity scattering by atom arrays, Phys. Rev. Lett. 131, 253603 (2023).

[79] S. J. Masson, J. P. Covey, S. Will, and A. Asenjo-Garcia, Dicke superradiance in ordered arrays of multilevel atoms, PRX Quantum 5, 010344 (2024).

[80] M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Generalized many-channel conductance formula with application to small rings, Phys. Rev. B 31, 6207 (1985).

[81] R. Landauer, Spatial variation of currents and fields due to localized scatterers in metallic conduction, IBM J. Res. Dev. 1, 223 (1957).
[82] D. Segal and A. Nitzan, Spin-boson thermal rectifier, Phys. Rev. Lett. **94**, 034301 (2005).

[83] E. Ressayre and A. Tallet, Holstein-Primakoff transformation for the study of cooperative emission of radiation, Phys. Rev. A **11**, 981 (1975).

[84] B. Bhandari, P. A. Erdman, R. Fazio, E. Paladino, and F. Taddei, Thermal rectification through a nonlinear quantum resonator, Phys. Rev. B **103**, 155434 (2021).

[85] D. Farina, G. De Filippis, V. Cataudella, M. Polini, and V. Giovannetti, Going beyond local and global approaches for localized thermal dissipation, Phys. Rev. A **102**, 052208 (2020).

[86] Meng Xu, J. T. Stockburger, and J. Ankerhold, Heat transport through a superconducting artificial atom, Phys. Rev. B **103**, 104304 (2021).

[87] H. T. Quan, Y. X. Liu, C. P. Sun, and F. Nori, Quantum thermodynamic cycles and quantum heat engines, Phys. Rev. E **76**, 031105 (2007).

[88] P. Pietzonka and U. Seifert, Universal trade-off between power, efficiency, and constancy in steady-state heat engines, Phys. Rev. Lett. **120**, 190602 (2018).

[89] N. Linden, S. Popescu, and P. Skrzypczyk, How small can thermal machines be? The smallest possible refrigerator, Phys. Rev. Lett. **105**, 130401 (2010).

[90] A. Mari and J. Eisert, Cooling by heating: Very hot thermal light can significantly cool quantum systems, Phys. Rev. Lett. **108**, 120602 (2012).

[91] J. Johansson, P. Nation, and F. Nori, QUTIP 2: A PYTHON framework for the dynamics of open quantum systems, Comput. Phys. Commun. **184**, 1234 (2013).

[92] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Open quantum systems with local and collective incoherent processes: Efficient numerical simulations using permutational invariance, Phys. Rev. A **98**, 063815 (2018).