Photonic heat transport across a Josephson junction

TL;DR

This paper investigates how photonic heat is transported across a Josephson junction in various circuit configurations, analyzing effects of anharmonicity and coupling strength to inform quantum thermodynamics experiments.

Contribution

It provides a comprehensive analysis of photonic heat transport in Josephson junctions, including effects of anharmonicity and different coupling regimes, with analytical approximations for heat flux.

Findings

Analytical expressions for photonic heat flux in different regimes

Effects of junction anharmonicity on heat transport

Potential for quantum thermodynamics experiments using Josephson circuits

Abstract

We present a detailed study of photonic heat transport across a Josephson junction coupled to two arbitrary linear circuits having different temperatures. First, we consider the linear approximation, in which a nonlinear Josephson potential is replaced by a quadratic one and the junction acts as an inductor. Afterwards, we discuss the effects of junction anharmonicity. We separately consider the weak-coupling limit, in which the Bloch band structure of the junction energy spectrum plays an important role, and the opposite strong-coupling regime. We apply our general results to two specific models: a Josephson junction coupled to two Ohmic resistors and two resonators. We derive simple analytical approximations for the photonic heat flux in many limiting cases. We demonstrate that electric circuits with embedded Josephson junctions provide a useful platform for quantum thermodynamics experiments.

Equations

$$I_C=I_C(0)\left|\cos (\pi \Phi /{\Phi }_0)\right|$$

$$J={\int }_0^{\infty }\frac{d\omega }{2\pi }\hslash \omega \tau (\omega )\left[N_2(\omega )-N_1(\omega )\right].$$

$$\hat{H}={\hat{H}}_J+{\hat{H}}_1+{\hat{H}}_2+{\hat{H}}_1^{\mathrm{int}}+{\hat{H}}_2^{\mathrm{int}}.$$

$${\hat{H}}_J=-4E_C\frac{{\partial }^2}{\partial {\phi }^2}+E_J(1-\cos \phi )$$

$${\hat{H}}_j=\sum _k\left[\frac{{\hat{P}}_{j,k}^2}{2M_{j,k}}+\frac{M_{j,k}{\omega }_{j,k}^2{\hat{X}}_{j,k}^2}{2}\right],j=1,2,$$

$${\hat{H}}_j^{\mathrm{int}}=-\sum _k\left(c_{j,k}{\hat{X}}_{j,k}\hat{\phi }-\frac{c_{j,k}^2}{2M_{j,k}{\omega }_{j,k}^2}{\hat{\phi }}^2\right),$$

$$\frac{\pi }{2}\sum _k\frac{c_{j,k}^2}{M_{j,k}{\omega }_{j,k}}\delta (\omega -{\omega }_k)=\frac{{\hslash }^2\omega }{4e^2}\mathrm{Re}\left[\frac{1}{Z_j(\omega )}\right].$$

$$J=\frac{d}{dt}⟨{\hat{H}}_1⟩=\frac{i}{\hslash }⟨[{\hat{H}}_1^{\mathrm{int}},{\hat{H}}_1]⟩=\sum _k\frac{c_{1,k}}{M_{1,k}}⟨{\hat{P}}_{1,k}\hat{\phi }⟩.$$

$$J=\frac{i}{\hslash }{\int }_{-\infty }^{t}d{t}^{\prime }\sum _k\frac{c_{1,k}}{M_{1,k}}\mathrm{tr}\left\{{\hat{P}}_{1,k}(t)\hat{\phi }(t)[\hat{\rho }(t^{\prime }),{\hat{H}}_{\mathrm{int}}(t^{\prime })]\right\}$$

$$=-\frac{i}{\hslash }{\int }_{-\infty }^{t}d{t}^{\prime }\sum _k\frac{c_{1,k}^2}{M_{1,k}}⟨[{\hat{X}}_{1,k}(t^{\prime })\hat{\phi }(t^{\prime }),{\hat{P}}_{1,k}(t)\hat{\phi }(t)]⟩.$$

$$J$$

$$S_{\phi }(\omega )=\int dt{e}^{i\omega t}⟨\hat{\phi }(t)\hat{\phi }(0)⟩.$$

$$S_Q(\omega )=\int dt{e}^{i\omega t}⟨\hat{Q}(t)\hat{Q}(0)⟩.$$

$$\frac{d\hat{\phi }}{dt}=\frac{i}{\hslash }[\hat{H},\hat{\phi }]=\frac{4E_C}{\hslash e}\hat{Q},$$

$$S_{\phi }(\omega )=\frac{16E_C^2}{e^2{\hslash }^2{\omega }^2}{S}_Q(\omega ).$$

$$J$$

$$Z_j(\omega )$$

$${\alpha }_j$$

$$\displaystyle\hat{H}=-\frac{\hbar\omega_{01}}{2}\hat{\sigma}_{z}+\sum_{j=1,2}\big{[}\hbar\omega_{j}\hat{b}^{\dagger}_{j}\hat{b}_{j}+\hbar g_{j}(\hat{b}^{\dagger}_{j}+\hat{b}_{j})\hat{\sigma}_{x}\big{]}.$$

$$g_j=\sqrt{\frac{Z_{rj}{\omega }_j{C}_j^2}{\pi C_{\Sigma }}}{\omega }_j.$$

$$Z_J(\omega )=-i\hslash \omega /2e{I}_C.$$

$$\mathrm{Re}\left[\frac{R_q}{Z_1({\omega }_J)}+\frac{R_q}{Z_2({\omega }_J)}\right]\gtrsim 1.$$

$$\tau (\omega )=\frac{4\mathrm{Re}\left[\frac{1}{Z_1(\omega )}\right]\mathrm{Re}\left[\frac{1}{Z_2(\omega )}\right]}{{\left|-i\omega C+\frac{1}{Z_1(\omega )}+\frac{1}{Z_2(\omega )}+\frac{1}{Z_J(\omega )}\right|}^2}.$$

$$\frac{g_1^2{\omega }_J^2}{{\omega }_n^2[({\omega }_J-{\omega }_n{)}^2+{\gamma }_1^2]}\ll 1$$

$$\tau (\omega )\approx \frac{{\gamma }_J^2}{(\omega -{\omega }_J{)}^2+{\gamma }_J^2}$$

$$+\sum _{n=0}^{\infty }\frac{a_n{\gamma }_1^4}{{\left[(\omega -{\omega }_n{)}^2+{\gamma }_1^2-\frac{a_n^2{\omega }_n^2}{4}\right]}^2+a_n{\gamma }_1^4}.$$

$${\gamma }_J=\frac{2g_1^2{\omega }_J^2{\gamma }_1}{{\omega }_1^2[({\omega }_J-{\omega }_{n_0}{)}^2+{\gamma }_1^2]},$$

$$a_n=\frac{4g_1^2{\omega }_n^2}{{\omega }_1^2{\gamma }_1^2|{\omega }_n^2-{\omega }_J^2|}.$$

$$J$$

$$\tau (\omega )$$

$$J$$

$$\tau (\omega )=\pi {\gamma }_J\delta (\omega -{\omega }_J)+\frac{4R_1{R}_2{C}_1^2{C}_2^2}{C_{\Sigma }^2}{\omega }^2,$$

Full text

I Introduction

Transport of heat in nanostructures is a subject of intense research Giazotto1 ; Hanggi ; Dubi ; Sai2016 . In most nanoscaled devices, the heat is transferred by either electrons, phonons, or photons. Photonic heat transport mechanism often dominates over the two other mechanisms Cleland at temperatures below 100 mK. In addition, it provides convenient way of transmitting tiny amounts of heat over macroscopic distances Timofeev ; Matti . Hence, proper understanding of photonic heat transport is essential, for instance, for development and calibration of highly sensitive low temperature radiation detectors.

Photonic heat flux can be accurately controlled by a tunable element embedded in the electric circuit. The natural choice of such an element for low-temperature superconducting circuits is a SQUID (superconducting quantum interference device) loop with the Josephson critical current adjusted by magnetic fluxMatthias ; Alberto . Such systems are promising platforms for realizing quantum thermal machines Sai2016 ; Hofer2016 . Motivated by these considerations, in this work we theoretically study photonic heat transport through a system schematically shown in Fig. 1. It contains two linear circuits, playing the role of thermal baths, with the impedances $Z_{1}(\omega)$ and $Z_{2}(\omega)$, which include dissipative elements having the temperatures $T_{1}$ and $T_{2}$, respectively. These circuits are connected via a symmetric SQUID with the critical current modulated by magnetic flux $\Phi$,

[TABLE]

where $\Phi_{0}$ is the flux quantum. The setup of Fig. 1 is analogous to the one used in the experiment Alberto , in which a SQUID has been coupled to two resonators terminated by Ohmic resistors. It also resembles the usual setup of circuit quantum electrodynamics experiments with transmon qubitsKoch , which is suitable for heat transport experiments with quantum systems. For example, the transition from quantum to classical behavior in a Josephson junction subject to thermal radiation with increasing temperature has been observed in a similar setupWallraff .

The system depicted in Fig. 1 also provides a useful test ground for theoretical predictions about quantum heat transport. Indeed, the SQUID is an example of a quantum system coupled to two thermal baths having different temperatures. It is described by a Hamiltonian of a particle moving in a one-dimensional periodic potential with the energy spectrum given by a set of Bloch bands. The Bloch band structure is very sensitive to the ratio between the Josephson energy of the SQUID $E_{J}$ and its charging energy $E_{C}$. Depending on the ratio $E_{J}/E_{C}$ and the temperatures $T_{1}$ and $T_{2}$, the model can be approximately reduced to that of weakly scattered particle, quantum oscillator or two level systemOjanen ; Saito ; Ren_2 ; Segal ; Yang ; Saito2 coupled to the two thermal baths. Important related topics are electronic cooling and heat rectification in nanostructuresGiazotto1 ; SN ; Hanggi . For example, recently a cooler based on a voltage biased Josephson junction has been proposedAnkerhold_1 and heat rectification in nonlinear Josephson circuits has been discussedAnkerhold_2 . Thus the physics of the system shown in Fig. 1 is very rich.

Approximately replacing the nonlinear Josephson junction by a linear inductor, one can express the photonic heat flux flowing from the circuit 2 into the circuit 1 in the usual form resembling Landauer formula,

[TABLE]

Here, $N_{1,2}(\omega)=1/(e^{\hbar\omega/k_{B}T_{1,2}}-1)$ are Bose functions and $\tau(\omega)$ is the photon transmission probability. It equals to the squared absolute value of the microwave transmission coefficient between the ports 1 and 2, $\tau(\omega)=|S_{12}(\omega)|^{2}$, which can be independently measured in the experiment. In Eq. (2), we use the sign convention $J>0$ for $T_{2}>T_{1}$. The transmission probability $\tau(\omega)$ can be found in a standard way by combining microwave impedances of various circuit elements Hekking , and in the linear approximation it does not depend on the temperatures $T_{1}$ and $T_{2}$.

Since the SQUID is a nonlinear, or anharmonic, quantum system, the Landauer formula (2) for the heat flux is just an approximation. There exist several methods of including the nonlinearity into the model. In the weak-coupling limit, one can describe the SQUID dynamics by rate equations, accounting for the jumps between its quantum states, and obtain the heat flux by counting the number of photons emitted to or absorbed from the thermal baths 1 and 2. Thanks to its simplicity, this method is very popular. It has been used, for example, in the analysis of experimental results of Ref. Alberto , in the theoretical description of the heat transport through a superconducting microwave cavity Tan , of thermal rectification in a quantum system with discrete energy spectrum SN , of Berry phase effects in the heat pumping Ren , etc. More general approach, which is also valid in the strong-coupling limit, has been developed by Ojanen and Heikkilä Ojanen_2 and Ojanen and Jauho Ojanen_1 , who have used the formalism of Keldysh Green’s functions, similar to the one developed earlier for the description of the electron transport through quantum dots MW . This formalism has been used to develop the theory of thermal rectification in quantum systems Ojanen_3 , to derive bounds on thermal conductanceSegal_3 , and to study heat transport through a two level system Ojanen ; Saito ; Ren_2 ; Segal ; Yang ; Saito2 . An equivalent formally exact approach to the heat transport problems is based on path integral techniques Tero ; QF ; Aurell_2 ; Aurell_3 . In principle, it allows one to go beyond the average value of the heat flux and to study the full counting statistics of the transferred heat. Yet another exact and numerically efficient method relies on stochastic Liouville-von Neumann equation with dissipation Ankerhold_2 ; Tuorila . Interestingly, in many cases one can express the exact heat flux flowing through the nonlinear system in the Landauer form (2). However, in this case the transmission probability $\tau(\omega)$ becomes dependent on the bath temperatures $T_{1,2}$. Such dependence, for example, makes the rectification of heat possible in a nonlinear and nonsymmetric system, as has been demonstrated in a recent experiment with Josephson junctions Senior2019 .

In this work, we present a detailed study of heat transport through a Josephson junction, or a SQUID, embedded in the circuit depicted in Fig. 1. We use the formalism developed in Refs. Ojanen_2 ; Ojanen_1 combining it with the well known results from the quantum theory of Josephson junctions AZL ; Zorin ; Averin ; SZ . We consider weak and strong coupling regimes, large and small values of the Josephson energy, and various frequency dependencies of the environment impedances $Z_{1,2}(\omega)$. In many limiting cases, we derive analytical expressions both for the photon transmission probability $\tau(\omega)$ and for the heat flux (2). The paper is organized as follows. In Sec. II, we introduce the model and the system Hamiltonian; in Sec. III, we discuss the linear approximation replacing the junction by an inductor; in Sec. IV, we consider the limit of weak coupling; in Sec. V, we analyze the strong-coupling regime, and in Sec. VI, we summarize the results. Some details of the calculations and auxiliary information are moved to Appendices. In Secs. IV and V, we separately discuss the limits of large, $E_{J}\gg E_{C},k_{B}T_{1,2}$, and small, $E_{J}\ll E_{C},k_{B}T_{1,2}$, Josephson energy .

II Model

We consider a symmetric SQUID with the critical current tunable by magnetic flux (1), which we will in the following call junction for simplicity, coupled to the two linear circuits, characterized by frequency-dependent impedances $Z_{1}(\omega)$ and $Z_{2}(\omega)$, as shown in Fig. 1. The linear circuits form two thermal baths with the temperatures $T_{1}$ and $T_{2}$. The system is described by the Hamiltonian

[TABLE]

Here

[TABLE]

is the Hamiltonian of the junction. We model the linear circuits 1 and 2 as oscillator baths with the Hamiltonians

[TABLE]

where $\hat{P}_{j,k}$, $\hat{X}_{j,k}$, $\omega_{j,k}$ and $M_{j,k}$ are, respectively, momentum, coordinate, frequency and mass of the $k-$th oscillator in the thermal bath $j$. The interaction Hamiltonians $\hat{H}_{j}^{\rm int}$ have the formCL ; Grabert

[TABLE]

where the coupling constants $c_{j,k}$ should be chosen in such a way that the bath spectral densities are proportional to the real parts of the inverse impedances,

[TABLE]

Formally exact expression for the heat flux $J$ in terms of Green’s functions has been derived in Refs. Ojanen_2 ; Ojanen_1 . Here we re-derive this expression for our particular setup using perturbation theory. We define the heat flux as the time derivative of the total energy stored in the environment 1 and write it in the form

[TABLE]

The angular brackets here stand for the quantum mechanical averaging with the density matrix of the whole system $\hat{\rho}$, $\langle\hat{A}\rangle=\,{\rm tr}\{\hat{A}\hat{\rho}\}$. The expression (8) can be further transformed in the interaction representation, in which the operators acquire time dependence $\hat{A}\to\hat{A}(t)=e^{i\hat{H}_{0}t/\hbar}\hat{A}e^{-i\hat{H}_{0}t/\hbar}$ with $\hat{H}_{0}=\hat{H}_{J}+\hat{H}_{1}+\hat{H}_{2}$. In this representation the system density matrix satisfies the evolution equation $i\hbar\,d\hat{\rho}(t)/dt=[\hat{H}_{\rm int}(t),\hat{\rho}(t)]$, where $\hat{H}_{\rm int}(t)=\hat{H}_{1}^{\rm int}(t)+\hat{H}_{2}^{\rm int}(t)$. Integrating this equation over time, we get $\hat{\rho}(t)=\hat{\rho}(0)-i\int_{0}^{t}dt^{\prime}[\hat{H}_{\rm int}(t^{\prime}),\hat{\rho}(t^{\prime})]/\hbar$. Substituting this expression in Eq. (8) and taking the long time limit, we find

[TABLE]

Note that the term containing $\hat{\rho}(0)$ vanishes at sufficiently long time $t$, at which the information about the initial state of the dissipative system is lost. For such times one can also extend the integration over $t^{\prime}$ as follows, $\int_{0}^{t}dt^{\prime}\to\int_{-\infty}^{t}dt^{\prime}$. Furthermore, within the framework of perturbation theory one can factorize the averages of the products of four operators into the products of the pairwise averages $\langle\hat{X}_{1,k}(t^{\prime})\hat{P}_{1,k}(t)\rangle$ and $\langle\hat{\varphi}(t)\hat{\varphi}(t^{\prime})\rangle$. The average $\langle\hat{X}_{1,k}(t^{\prime})\hat{P}_{1,k}(t)\rangle$ can be evaluated, since in the lowest order of the perturbation theory, the oscillators of the thermal baths do not interact with the junction. Performing these operations, we arrive at the expression

[TABLE]

Here we have introduced the Fourier transformed phase-phase correlation function

[TABLE]

Although we have used perturbation theory while deriving the formula (9), it is actually exact. One can prove this result by either Keldysh Green’s function technique Ojanen_2 ; Ojanen_1 , or by path integral techniqueQF ; Aurell_2 ; Aurell_3 .

One can alternatively express the heat flux in terms of the correlation function of the operators of the charge accumulated in the junction capacitance $\hat{Q}=-2ie(\partial/\partial\varphi)$,

[TABLE]

Since the operators $\hat{\varphi}(t)$ and $\hat{Q}(t)$ are related by the equation of motion for the phase,

[TABLE]

the correlation function (10) can be written in the form

[TABLE]

Accordingly, the heat flux (9) can be expressed as

[TABLE]

Equations (9) and (14) are the starting points of our analysis. Employing various approximations, we will derive approximate expressions for the heat flux $J$ without specifying particular frequency dependence of the impedances $Z_{1,2}(\omega)$. Afterwards we will consider two specific examples. One example is the case of Ohmic dissipation with $Z_{j}(\omega)=R_{j}$, which is illustrated in Fig. 2a. Another example, shown in Fig. 2b and inspired by recent heat transport experiment Alberto , is the junction capacitively coupled to the two $\lambda/4$ resonators having characteristic impedances $Z_{r1}$ and $Z_{r2}$ and terminated by Ohmic resistors $R_{1}$ and $R_{2}$. In this case, the impedances of the two linear electric circuits are

[TABLE]

Here, $C_{j}$ are the coupling capacitors [see Fig. 2(b)], $\omega_{j}$ are the frequencies of the fundamental modes of the resonators, and the parameters $\alpha_{j}$ determine their internal quality factors, $Q_{j}=\pi/4\alpha_{j}$. Here we will only consider the regime $R_{j}<Z_{rj}$. Provided $Q_{j}\gg 1$, $E_{J}\gg E_{C}$ and the temperatures $T_{1,2}$ are sufficiently low, one can keep only the two lowest energy levels formed close to the bottom of the Josephson potential well, and approximately describe the system by Rabi Hamiltonian Koch

[TABLE]

Here the transition frequency between the levels is $\hbar\omega_{01}=\hbar\omega_{J}-E_{C}$, $\omega_{J}=\sqrt{8E_{J}E_{C}}/\hbar$ is the classical frequency of small oscillations at the bottom of the potential well, and the coupling constants $g_{1},g_{2}$ are given by

[TABLE]

Here we defined the total capacitance $C_{\Sigma}=C_{1}+C_{2}+C$.

III Linearized dynamics

In this section, we approximately replace the nonlinear Josephson junction by an inductor with the impedance

[TABLE]

This simple approximation is valid in the wide range of parameters. For the beginning, it is valid for sufficiently low temperatures, $k_{B}T_{1,2}\lesssim 2E_{J}$, sufficiently high ratio $E_{J}/E_{C}\gtrsim 1$, and sufficiently strong coupling between the junction and the environment,

[TABLE]

Here we have introduced the resistance quantum $R_{q}=h/e^{2}$. The condition (20) ensures that the width of the energy levels formed at the bottom of the Josephson potential well, which is determined by the transition rates (54), exceeds the charging energy $E_{C}$. Since the latter defines the scale of anharmonicity in the system the junction may be viewed as a linear element. The linear approximation is also valid at high temperatures $k_{B}T_{1,2}\gg E_{J}$, where one can just put $I_{C}=0$, thus removing the nonlinear element from the circuit.

Once the Josephson junction has been replaced by a linear lumped element, one can exactly reduce the full quantum problem to the solution of the classical Langevin equations for the Josephson phase and currents, which contain the stochastic noises generated by the circuit elementsHekking ; Schmid . Similar Langevin equations have been used, for example, in order to describe the transport of heat by phonons in harmonic lattices Dhar . The details of the analysis are given in Appendix A. The final expression for the heat flux acquires the Landauer form (2) with the photon transmission probability

[TABLE]

Here the junction impedance is given by Eq. (19). Since $\tau(\omega)$ does not depend on the temperatures of the two baths, the heat flux (2) has the property $J(T_{1},T_{2})=-J(T_{2},T_{1})$, which implies the absence of heat rectification in the linear approximation. Absence of rectification is the well known property of linear harmonic systems SN ; Dhar1 .

As an example, in Fig. 3, we plot the transmission probability (21) versus frequency $f=\omega/2\pi$ for a junction coupled to the two identical resonators characterized by the impedances (15). We observe that $\tau(\omega)$ has peaks at frequencies corresponding to the eigen-modes of the resonators $\omega_{n}=(2n+1)\omega_{1}$, with $n=0,1,2\dots$. Besides that, there exists a peak at the frequency of phase oscillations in the junction $\omega_{J}$. For a certain value of magnetic flux the resonant condition $\omega_{J}=\omega_{1}=\omega_{2}$ is achieved. Close to this value three peaks appear in $\tau(\omega)$ due to hybridization of the modes of the resonators and of the junction. The positions of these peaks are shown in Figs. 4a and 4b. For weak coupling (Fig. 4a) the hybridization of the modes is almost invisible, while for strong coupling (Fig. 4b) it is quite strong and manifests itself in avoided crossing of the peak position lines. In Figs. 4c and 4d we plot the dependence of the heat power (2) on magnetic flux for weak and strong coupling respectively.

In certain limiting cases one can derive analytical approximations for $\tau(\omega)$ and for the heat flux $J$. First, we consider the junction coupled to identical resonators with high quality factors. We assume that the junction mode is sufficiently far detuned from the modes of the resonators and ignore the shift of the junction frequency induced by the coupling to the resonators. This approximation is valid provided the condition

[TABLE]

is satisfied for all modes of the resonators. Here $\gamma_{1}=2\alpha_{1}\omega_{1}/(\pi+2Z_{r1}C_{1}\omega_{1})$ is the total damping rate of the resonators. In this case, the transmission probability $\tau(\omega)$ is expressed as a sum of well separated peaks,

[TABLE]

The first Lorentzian peak comes from the junction mode. Its width $\gamma_{J}$ is given by

[TABLE]

where $\omega_{n_{0}}$ is the frequency of the mode closest to the Josephson frequency $\omega_{J}$. The remaining peaks correspond to the modes of the resonators. They have non-Lorentzian shape, and their heights $a_{n}$ are determined by the detuning between the junction mode and the corresponding resonator mode,

[TABLE]

The resonator modes split into the pairs of closely lying peaks if $a_{n}\omega_{n}>2\gamma_{1}$, when the effective coupling between the modes of the resonators becomes stronger than the dissipation rate. At high temperatures, $k_{B}\max\{T_{1},T_{2}\}\gtrsim\min\{\hbar\omega_{J},\hbar\omega_{1}\}$, one can approximately replace all the peaks in Eq. (23) by $\delta-$functions. The heat flux (2) then takes the form

[TABLE]

The first term in this expression describes the heat transport through the junction mode, while the second term – through the modes of the resonators. The effective coupling between the modes of different resonators depends on the value of the Josephson frequency $\omega_{J}$.

Next, we consider arbitrary resonators, which are not necessarily identical, but assume that the temperatures are low, $k_{B}\max\{T_{1},T_{2}\}\lesssim\min\{\hbar\omega_{J},\hbar\omega_{1,2}\}$. In this case the transmission probability and the heat flux are

[TABLE]

In the regime $\hbar\omega_{J}\lesssim\hbar\omega,k_{B}\max\{T_{1},T_{2}\}\ll\hbar\omega_{1,2}$ we find

[TABLE]

In Fig. 5 we schematically overview various heat transport regimes for a Josephson junction coupled to two identical resonators with high quality factors.

Let us now consider the system with suppressed Josephson current, $I_{C}=0$. The transmission probability (21) for such system is plotted in Fig. 6 for a symmetric and an asymmetric coupling. In both cases it exhibits sharp peaks at frequencies corresponding to the eigenmodes of the resonators $\omega_{j,n}=(2n+1)\omega_{j}$. For the case of identical resonators one can derive an accurate high temperature asymptotics for the integral (2) valid at $k_{B}T_{1,2}\gtrsim\hbar\omega_{1}/2\pi$. In this case, one can average the transmission probability (21) over one period $2\omega_{1}$ and replace $\tau(\omega)$ in the Eq. (2) by its averaged value

[TABLE]

where $\omega_{0}=(Z_{r1}C_{1})^{-1}\sqrt{C_{\Sigma}/C}$. The frequency $\omega_{0}$ is usually very high so that the condition $k_{B}T_{1,2}\ll\hbar\omega_{0}$ is fulfilled. In this limit one finds

[TABLE]

where the parameter

[TABLE]

characterizes the thermal conductance between capacitively coupled resonators, and

[TABLE]

is the characteristic temperature at which the crossover between two different heat transport regimes occurs. The result (30) is consistent with the expression (26) in the limit $\omega_{J}=0$ and $R_{1}\ll Z_{r1}$. In Fig. 7 we illustrate the temperature dependence of the heat flux (2) for a symmetric system with the two identical resonators and zero Josephson current, $I_{C}=0$, and compare it with the approximation (30). We find that the approximation (30) indeed works quite well at high temperature. The low temperature limit is well described by the Eq. (28). For comparison, in the same figure we have also plotted the temperature dependence of the heat flux between two resonators with different parameters.

Let us now consider the system with ohmic dissipation, $Z_{j}(\omega)=R_{j}$. The transmission probability takes the form

[TABLE]

where $\gamma=(R_{1}+R_{2})/R_{1}R_{2}C$. In the low temperature regime $k_{B}T_{1,2}\ll\min\{\hbar\omega_{J},\hbar\omega_{J}^{2}/\gamma\}$ one finds the following approximation for the heat flux

[TABLE]

where $\alpha_{j}=h/e^{2}R_{j}$ are the dimensionless conductances of the ohmic resistors. The result (34) may be interpreted as the contribution of co-tunneling, and it has the same temperature dependence as the co-tunneling contribution to the heat flux through a two level system Ojanen . At high temperatures $k_{B}T_{1,2}\gtrsim\max\{\hbar\omega_{J},\hbar\gamma\}$ one finds

[TABLE]

Interestingly, this result is independent of $E_{J}$ even for temperatures $T_{1,2}$ well below the barrier height $2E_{J}$. At these temperatures one can roughly approximate the transmission probability (33) by a $\delta$-peak,

[TABLE]

and the heat flux by the expression

[TABLE]

This simple approximation correctly captures the leading term of the high temperature asymptotics (35), but fails to reproduce the sub-leading term as well as the low temperature power law dependence of the heat flux (34). Interestingly, the $\delta$-peak approximation (36) for the transmission probability holds even for an overdamped junction with $\hbar\omega_{J}\ll\hbar\gamma\ll k_{B}T_{1,2}$ and for a system with $I_{C}=0$ provided $\hbar\gamma\ll k_{B}T_{1,2}$. However, in these cases the $\delta$-peak shifts towards zero frequency. For an overdamped junction with $\gamma\gg\omega_{J}$ and for $\hbar\omega_{J}^{2}/\gamma\lesssim k_{B}T_{1,2}\lesssim\hbar\gamma$ the heat flux (2) takes the familiar form

[TABLE]

describing the heat transfer between two resistors. In Fig. 8 we provide a schematic diagram of the three heat transport regimes discussed above.

Finally, we consider the limit of sufficiently weak coupling and sufficiently high temperatures without specifying the frequency dependence of the impedances $Z_{j}(\omega)$. In this regime, one can replace the transmission (21) by a single $\delta$-peak. The corresponding expression reads SN

[TABLE]

Afterwards, the heat flux (2) acquires a simple formSN

[TABLE]

We test the accuracy of this approximation in Figs. 4 c and d, where we plot the heat flux through a system with two identical resonators versus the magnetic flux. We find that the approximation (40) indeed works very well in the weak coupling regime (Fig. 4c), but it fails at stronger coupling (Fig. 4d).

To conclude this section, we note that replacing the nonlinear Josephson junction by a linear inductor with the impedance (19) results in a very good approximation for the heat flux provided the coupling between the junction and the environment is sufficiently strong and the temperature is sufficiently high.

IV Effects of anharmonicity in the weak-coupling limit

In this section, we consider the effect of Josephson nonlinearity on the heat flux in the weak-coupling limit

[TABLE]

This condition should be valid for all relevant frequencies $\hbar\omega\lesssim k_{B}T_{1,2}$. It ensures that the transition rates between the eigenstates of the junction Hamiltonian (4) are smaller than the charging energy $E_{C}$, which provides the scale of anharmonicity of the junction. For a setup with two high-quality resonators, we also require sufficiently strong detuning between the Josephson frequency and the nearest mode of one of the resonators $\omega_{n_{0},j}$

[TABLE]

In this case the shift of the junction frequency caused by the coupling to the resonators is much smaller than $E_{C}$.

Under this conditions, it is convenient to work in the basis of the eigenfunctions $\psi_{qn}(\varphi)$ of the junction Hamiltonian, which satisfy Schrödinger equation

[TABLE]

Here $\epsilon_{n}(q)$ is the $2e$-periodic energy of the $n-$th Bloch band, which depends on the electric charge $q$ transferred through the junctionZorin ; Averin . According to the Bloch theorem the wave functions can be expressed in the form $\psi_{qn}(\varphi)=u_{qn}(\varphi)e^{iq\varphi/2e}$, where $u_{qn}(\varphi)$ are the $2\pi$-periodic functions normalized as $\int_{-\pi}^{\pi}d\varphi u_{qm}^{*}(\varphi)u_{qn}(\varphi)=\delta_{mn}$. In order to find the transition rates between the energy bands, we need to know the matrix elements of the momentum $\hat{p}=-i(\partial/\partial\varphi)$ and of the phase between the functions $u_{qn}(\varphi)$ belonging to different bands, $m\not=n$,

[TABLE]

Equation (12) leads to the relation between them

[TABLE]

We will also use the well known result of solid state physics, which states that the operator $\hat{p}$ within one band acts on any function of $q$ by multiplying it with the combination $(e/4E_{C})(\partial\epsilon_{n}(q)/\partial q)$. More detailed information about various analytical approximations for the Bloch band energies $\epsilon_{n}(q)$ and the matrix elements $p_{mn},\varphi_{mn}$ can be found in Refs. Averin ; Koch and in Appendix B.

With the matrix elements at hand, we can express the charge-charge correlation function (11) as a sum of intraband [$S_{Q}^{\rm b}(\omega)$] and interband [$S_{Q}^{\rm ib}(\omega)$] contributions,

[TABLE]

The intraband contribution has the form

[TABLE]

while the contribution associated with the interband transitions reads

[TABLE]

Here we have defined the interband frequency $\omega_{mn}(q)=[\epsilon_{m}(q)-\epsilon_{n}(q)]/\hbar$ and the occupation probability $w_{n}(q)$ of the quantum state described by the wave function $\psi_{nq}(\varphi)$. These occupation probabilities are normalized as

[TABLE]

The weak coupling approximation has been explicitly used in deriving Eq. (49), in which the interaction between the junction and the thermal baths has been ignored. In contrast, at this stage we keep the intraband correlation function (48) in a general form.

Let us assume that the impedances $Z_{j}(\omega)$ remain finite in the low-frequency limit $\omega\to 0$. In this case the distribution function $w_{n}(q)$ satisfies the kinetic equation Zorin ; Averin

[TABLE]

valid in the lowest nonvanishing order of the perturbation theory in the interaction Hamiltonians $\hat{H}_{1,2}^{\rm int}$. Here we have introduced the shunt resistance $R_{S}^{-1}=Z_{1}^{-1}(0)+Z_{2}^{-1}(0)$ and the effective temperature of the junction

[TABLE]

The interband transition rates, appearing in the last term of Eq. (51), are given by the sum of partial contributions of the baths 1 and 2, respectively,

[TABLE]

and the latter are defined as

[TABLE]

Since the matrix elements $\varphi_{mn}$ do not exceed 1, the condition (41) ensures that $\Gamma^{(j)}_{mn}+\Gamma^{(j)}_{nm}\ll\hbar|\omega_{mn}|$ and guaranties the validity of the kinetic equation (51).

Ignoring the interaction effects, we derive the following expressions for the symmetric and antisymmetric parts of the intraband correlation function (48):

[TABLE]

The first of these equations follows from the conservation of the charge $q$ in the absence of interactions. As a result, the symmetrized correlator $\langle[\partial\epsilon_{n}(\hat{q}(t))/\partial q]\,[\partial\epsilon_{n}(\hat{q}(0))/\partial q]+[\partial\epsilon_{n}(\hat{q}(0))/\partial q]\,[\partial\epsilon_{n}(\hat{q}(t))/\partial q]\rangle$ becomes time independent and its Fourier transform reduces to a $\delta-$function peaked at zero frequency. However, in order to treat low-frequency asymptotics correctly, one should remember that this peak has small finite width. The second equation follows from the fluctuation dissipation theorem, which in equilibrium, i.e., at $T_{1}=T_{2}=T$, states that $S_{Q}^{\rm b}(\omega)-S_{Q}^{\rm b}(-\omega)=\tanh({\hbar\omega}/{2k_{B}T})[S_{Q}^{\rm b}(\omega)+S_{Q}^{\rm b}(-\omega)]$. Since the symmetrized correlator is proportional to $\delta(\omega)$, the tangent in front of it can be expanded at small frequencies. Afterwards, one can make a replacement $T\to T_{J}$ because in the low frequency Markovian limit the intensities of the noises generated by the two environments become proportional to the corresponding temperatures and also become additive.

Substituting the correlation function (55) in the general formula (14), we obtain the heat flux as a sum of the intraband and interband contributions,

[TABLE]

The intraband contribution has the form

[TABLE]

while the interband part reads

[TABLE]

It is instructive to consider the system in which the frequency dependence of the two bath spectra is the same. More precisely, let us assume that

[TABLE]

where $a$ is a frequency independent constant. We introduce the heat flux flowing in the direction $1\to 2$, which we will denote as $J_{2}$. It is given be the expressions (56,57,58) with the interchanges indexes $1\leftrightarrow 2$. Energy conservation in the stationary case implies $J+J_{2}=0$. This property is ensured by the kinetic equation (51). Using the energy conservation condition we can re-write the heat flux in the form

[TABLE]

After such symmetrization, well known in the theory of quantum dots MW ; Konig , the total heat flux $J$ acquires the Landauer form (2) with the transmission probability given by the sum of intraband in interband contributions,

[TABLE]

These contributions read

[TABLE]

Thus, in the weak coupling limit the intraband transitions produce a narrow low frequency peak in the transmission probability, while the interband transitions result a series of transmission bands at higher frequencies. Unlike the transmission probability (21) derived in linear approximation, the probabilities (62) and (63) may depend on temperatures $T_{1}$ and $T_{2}$ via the distribution function $w_{n}(q)$. Hence, in general, $J(T_{1},T_{2})\not=-J(T_{2},T_{1})$ and heat rectification becomes possible. We also note that in the limit $E_{J}=0$, in which the Bloch bands are reduced to parabolas $(q-2en)^{2}/2C$, the intraband transmission probability (62) takes the simple form (36) with $\omega_{J}=0$, while the interband contribution (63) vanishes.

Another simple case is the system with strongly asymmetric coupling. Namely, if we assume that ${\rm Re}[Z_{1}^{-1}(\omega)]\ll{\rm Re}[Z_{2}^{-1}(\omega)]$, then the stationary solution of the kinetic equation (51) can be found exactly,

[TABLE]

In this regime the junction is thermalized with the thermal bath 2. For a strongly asymmetric system the expressions for transmission probabilities (62,63) are valid for arbitrary frequency dependence of the impedances $Z_{j}(\omega)$.

Even though the kinetic equation (51) and the heat fluxes (57,58) have relatively simple form, the solution of the problem for an arbitrary ratio $E_{J}/E_{C}$ is a challenging task, which requires numerical simulations. In the next two subsections, we consider the limits $E_{J}\gg E_{C}$ and $E_{J}\ll E_{C}$, where further approximations are possible.

IV.1 Tight-binding limit $E_{J}\gg E_{C}$.

In the limit $E_{J}/E_{C}\gg 1$, the lowest Bloch bands acquire cosine dispersion typical for tight-binding models (see Fig. 9 for illustration)

[TABLE]

Here, $E_{n}$ is the position of the $n$-th energy level in an isolated potential well and $\delta_{n}$ is the half-bandwidth proportional to the hopping amplitude between the wells. Analytical approximations for both these parameters are summarized in Appendix B. The cosine dispersion (65) is valid for the lowest energy bands with $0\leq n\leq n_{\max}$, where $n_{\max}$, defined by the Eq. (130), is the number of the highest energy level lying under the barrier $2E_{J}$. An example is shown in Fig. 9.

Next, we assume that $k_{B}T_{1,2}\lesssim 2E_{J}$, and that the coupling between the junction and the environment is not too weak. Namely, we impose the condition

[TABLE]

for all bands relevant for the transport of heat, i.e., those with $E_{n}\lesssim k_{B}T_{1,2}$. In this case, the transition rates between the levels (54) exceed the half-bandwidth $\delta_{n}$ and one can put $\delta_{n}=0$. Afterwards, the kinetic equation (51) acquires a simple form

[TABLE]

where $W_{n}=\int_{-e}^{e}dq\,w_{n}(q)$ is the total population of the $n-$th band. The interband contribution to the heat flux (58) is simplified to

[TABLE]

The small intraband heat flux, which may dominate at low temperatures, will be analyzed later. For a system with baths spectra having the same frequency dependence (59), or for a strongly asymmetric system, the transmission probability associated with interband transitions (63) acquires the form

[TABLE]

Thus every possible transition between the levels results in a separate $\delta$ peak in the transmission probability at frequency $\omega_{mn}$ corresponding to the level spacing.

The energy levels close to the bottom of the Josephson potential are not equidistant due to the quartic nonlinearity of the Josephson potential [see Eq. (132)], and the splittings between the lowest neighboring levels vary as $\hbar\omega_{n+1,n}=\hbar\omega_{J}-(n+1)E_{C}$. In Fig. 10a, we plot the dependence of the frequencies $f_{n+1,n}=\omega_{n+1,n}/2\pi$ on the magnetic flux for $n=0,1,2$ for a junction coupled to two identical resonators with the impedances (15). Every time one of these frequencies crosses the frequency of the resonators, a peak in the heat flux (68) appears, as shown in Fig. 10b. The height of these peaks is determined by temperatures $T_{1,2}$ via the occupation probabilities $W_{n}$. At the lowest temperature only one peak, coming from the transitions between the levels 0 and 1, survives. Since the transmission probability (69) depends on temperatures $T_{1,2}$ due to the anharmonicity of the junction, one should expect rectification in the system provided the couplings to the two baths are different. In Fig. 11, we show that this is indeed the case. For the chosen parameters the difference between the heat fluxes $J(T_{1},T_{2})$ and $-J(T_{2},T_{1})$ reaches up to 50%. In order to establish the correspondence with the previous section, we note that one can find an exact solution of the master equation (67) for the junction with the equidistant energy levels $E_{n}=\hbar\omega_{J}(n+1/2)$ corresponding to the spectrum of a harmonic oscillator, see Appendix C. For this system the heat flux (68) takes the form (40), while all the peaks in the transmission probability (69) collapse into a single peak and it acquires the form (39). For comparison, we have plotted the harmonic approximation for the heat flux (40) in Fig. (10). Clearly, for the chosen parameters it significantly differs from the more accurate result (68).

For small, but finite, bandwidths $\delta_{n}$ one can evaluate the integral over $q$ in the general expression for the transmission probability (63) neglecting weak $q$ dependence of the occupation probabilities $w_{n}$ and matrix elements $\varphi_{mn}$. Afterwards, $\tau^{\rm ib}(\omega)$ takes the form (69) with the $\delta$ functions replaced by finite width transmission bands having a double peak shape with square root divergences,

[TABLE]

Next, we consider the intraband contribution to the heat flux and transmission probability in the tight-binding limit $E_{J}\gg E_{C}$. For simplicity, we consider only strongly asymmetric systems with ${\rm Re}\left[Z_{1}^{-1}(\omega)\right]\ll\,{\rm Re}\left[Z_{2}^{-1}(\omega)\right]$, for which an analytical solution of the problem is possible. Taking Boltzmann distribution function (64), in which $\epsilon_{n}(q)$ have cosine form (65), we evaluate the intraband heat flux (57) and arrive at the result

[TABLE]

Here $I_{n}(x)$ are modified Bessel functions. The intraband heat flux (71) may dominate over the interband one (68) at low temperatures $k_{B}T_{2}\lesssim\hbar\omega_{J}/2\pi$ if $Z_{1}^{-1}(0)$ is sufficiently large. At these temperatures, only the lowest energy band is populated, and we can approximate $J^{\rm b}$ as

[TABLE]

One can also derive relatively simple analytical expressions for the heat flux in the temperature interval $\hbar\omega_{J}/2\pi\lesssim k_{B}T_{2}\lesssim 2E_{J}$. In this case, the interband contribution can be well approximated by harmonic result (40), which can be further simplified to the form

[TABLE]

resembling the result (35) derived for the case of Ohmic dissipation. The intraband contribution (71) in this temperature interval takes the form

[TABLE]

The heat flux (74) is a contribution of thermally activated phase slips, i.e. the jumps between the neighboring potential wells of the Josephson potential, which become relevant above the quantum to classical crossover temperature of the junction $T^{*}=\hbar\omega_{J}/2\pi$. The phase slip heat flux (74) may dominate over the interband contribution (73) if ${\rm Re}[Z^{-1}(0)]\gg{\rm Re}[Z^{-1}(\omega_{J})]$.

IV.2 Weak Josephson coupling, $E_{J}\ll E_{C}$

For $E_{J}\ll E_{C}$ the energies of the two lowest Bloch bands with $n=0$ and 1 are well approximated by GR ; Jared

[TABLE]

These two lowest Bloch bands are illustrated in Fig. 12.

Higher lying Bloch bands can be replaced by parabolas $4E_{C}(n-q/2e)^{2}$ and small energy gaps between them can be ignored. Furthermore, the matrix elements of the phase operator connecting the two bands read Averin

[TABLE]

For simplicity, here we only consider a system with strongly asymmetric coupling, for which the distribution function has Boltzmann form (64) with the temperature $T_{2}$. For $E_{J}\ll E_{C}$, the contribution of the interband transitions to the the transmission probability (63) becomes

[TABLE]

where $\theta(x)$ is the Heaviside step function. The transmission probability (77) exhibits a gap $E_{J}/\hbar$ at low frequencies and gets exponentially suppressed for $k_{B}T_{2}\lesssim E_{C}$.

For temperatures exceeding the gap, $k_{B}T_{1,2}\gg E_{J}$, one can derive rather simple analytical expression for the heat flux. The intraband contribution (57) takes the form

[TABLE]

The leading term in this expression is similar to the heat flux through a capacitor (35), while the correction $\propto E_{J}$ comes from the opening of the energy gap between the two lowest Bloch bands. The interband contribution to the heat flux is small in this limit,

[TABLE]

V Effects of anharmonicity in the strong-coupling limit

In this section, we consider the strong-coupling limit (20) in more detail and derive corrections to the results of Sec. III caused by anharmonicity of the junction. As in the previous section, we will separately consider the regimes $E_{J}\gg E_{C}$ and $E_{J}\ll E_{C}$.

V.1 Tight-binding limit $E_{J}\gtrsim E_{C}$

In this section, we consider the tight-binding regime $E_{J}\gtrsim E_{C}$ and assume that the temperature is low, $k_{B}T_{1,2}\ll\hbar\omega_{J}$. We derive the correction to the Landauer formula (2) originating from the finite bandwidth of the lowest Bloch band. For this purpose, we perform the expansion in the half-bandwidth $\delta_{0}$ keeping only the lowest nonvanishing term $\propto\delta_{0}^{2}$. This correction may also be interpreted as a contribution of phase slips, or the jumps between the neighboring potential wells. In the weak-coupling limit and for $E_{J}\gg E_{C}$ the Landauer heat flux (2) crosses over to the interband contribution (68), while the phase slip correction $\propto\delta_{0}^{2}$ — to the high-temperature expansion of the intraband contribution (72).

The phase slip correction is important only at low temperatures $k_{B}T_{1,2}\ll\hbar\omega_{J}$ unless the Landauer contribution (2) is deliberately made small by specific choice of the circuit impedances $Z_{1,2}(\omega)$. At these temperatures only the lowest Bloch band is populated and the intraband charge correlation function (48) takes the form

[TABLE]

Here, we have introduced the function

[TABLE]

similar to the one appearing in the theory of the Coulomb blockade in a Josephson junction embedded in electromagnetic environment ANO . The function ${\cal P}(\omega)$ has a physical meaning of the probability of photon absorption in either of the two thermal baths during a phase slip event. This function should be evaluated at $\delta_{0}=0$ and, since the temperature is low, one can also ignore interband transitions. Afterwards, applying the standard methods IN and adapting them to the environment consisting of the two baths with different temperatures, we obtain ${\cal P}(\omega)$ in the form of the convolution of the two probabilities ${\cal P}_{1,2}(\omega)$ describing photon absorption in the baths 1 and 2,

[TABLE]

These probabilities are given by the integrals

[TABLE]

in which the functions ${\cal F}_{j}(t)$ and ${\cal K}_{j}(t)$ read

[TABLE]

For details of the derivation, we refer the reader to the Appendix D. The substitution of the correlation function (80) in the general expression (14) for the heat flux results in the phase slip correction (see Appendix D for details)

[TABLE]

This expression has a clear physical meaning: the correction to the heat flux is the net contribution of elementary events, in which a photon is absorbed by the junction from the bath 2 and then re-emitted into the bath 1. The same expression has been recently derived for a two level system in Ref. Aurell .

Interestingly, the correction (86) can be transformed to the Landauer form (2) with the aid of the detailed balance relations

[TABLE]

The corresponding contribution to the transmission probability has the form

[TABLE]

This expression provides the generalization of the weak-coupling expression for the intraband transmission probability (62). Indeed, the low frequency $\delta$-peak now acquires a finite width, which depends on the bath spectra and may also depend on temperatures $T_{1,2}$. The validity condition of the expressions (86) and (88) is $\delta_{0}\ll k_{B}T_{1,2}\ll\hbar\omega_{J}$. However, for sufficiently strong coupling, one can use them even at $T=0$.

The computation of the functions ${\cal P}_{j}(\omega)$ for arbitrary impedances $Z_{1,2}(\omega)$ is complicated. However, in the weak-coupling limit (41), one can obtain relatively simple expressions by expanding the integrals in Eqs. (83) to the first order in small functions ${\cal F}_{j}(t)$ and ${\cal K}_{j}(t)$. In this approximation, one finds

[TABLE]

Hence, the correction to the transmission probability (88) simplifies to

[TABLE]

Let us now consider the two examples, which we have introduced before. First, we consider resonant environments with the impedances (15). For high quality-factor resonators with $R_{j}\ll Z_{rj}$, one can approximate the transmission probability (90) at $\omega\ll\omega_{1,2},\omega_{J}$ as

[TABLE]

where the first term comes form the transmission probability (21) and the second from the phase slip correction (88). The total heat flux takes the form

[TABLE]

The first term comes from the harmonic approximation (2,21) and dominates at relatively high temperatures $k_{B}T_{1,2}\gtrsim\sqrt{\delta_{0}E_{J}/\pi}$, where the junction acts as an inductor. In contrast, at low temperatures $k_{B}T_{1,2}\lesssim\sqrt{\delta_{0}E_{J}/\pi}$, the phase slip contribution (86) dominates the heat transport, and the junction more resembles a capacitor.

Our second example is the case of Ohmic environments with $Z_{j}(\omega)=R_{j}$. At low frequencies and temperatures $\hbar\omega,k_{B}T_{1,2}\ll\hbar\omega_{J}$, the functions (83) acquire the form IG

[TABLE]

Here we assumed that despite strong coupling to the environments the junction remained underdamped, which is possible if $\alpha_{1}+\alpha_{2}\lesssim\pi\sqrt{{8E_{J}}/{E_{C}}}.$ In this case, one obtains the transmission probability in the form

[TABLE]

where $\tau^{\rm b}(\omega)$ is given by Eq. (88) and the interband contribution results form the low frequency expansion of the transmission (21),

[TABLE]

One can now estimate the low temperature heat conductance of the system

[TABLE]

It reads

[TABLE]

where the interband part follows from Eq. (34)

[TABLE]

and the intraband one is given by

[TABLE]

The latter formula is valid for $\alpha_{1,2}\lesssim 4$. One can verify that for $\delta_{0}\ll k_{B}T\ll\hbar\omega_{J}$ and $\alpha_{1}\ll\alpha_{2}\ll 1$ Eq. (99) matches the thermal conductance in the weak coupling limit (72). It also agrees with the results of Refs.Saito ; Saito2 , where the heat transport through a two level system has been studied, if one identifies the parameter $\alpha$, introduced there, with the combination $(\alpha_{1}+\alpha_{2})/4$.

V.2 Heat transport at $E_{J}\ll k_{B}T_{1,2}$

In this section, we assume that $E_{J}$ is smaller than the bath temperatures $T_{1,2}$. Performing the expansion in the small parameter $E_{J}/k_{B}T$, we obtain the phase-phase correlation function in the form (see Appendix E for details)

[TABLE]

Here $Z_{J}(\omega)$ is the effective impedance of the Josephson junction in presence of strong phase fluctuations Olli

[TABLE]

and

[TABLE]

is the probability for the junction to emit a photon with the frequency $\omega$ into the environment. This function resembles the function ${\cal P}(\omega)$ defined in Eq. (81), but differs from it because the charge $q$ and the phase $\varphi$ fluctuate in a different way. The absorption probabilities of the two baths, $P_{j}(\omega)$, are given by Eqs. (83) with the functions ${\cal F}_{j}(t)$, ${\cal K}_{j}(t)$ replaced by similar functions describing phase fluctuations,

[TABLE]

Here the effective spectra of the environments ${\cal R}_{j}(\omega)$ are

[TABLE]

The functions $F_{j}(t)$ and $K_{j}(t)$ are well known from the theory of environmental Coulomb blockadeIN ; Saclay ; Falci .

The substitution of the correlation function (100) in the general expression for the heat flux (9) gives the result

[TABLE]

where the elastic contribution to the heat flux $J_{1}^{\rm el}$ is given by Landauer formula (2) with the transmission probability (21) containing the modified junction impedance (101), and $J_{1}^{\rm in}$ describes the contribution coming from the inelastic scattering of photons on the junction,

[TABLE]

The inelastic heat flux (107) has the same form as the intraband contribution (86), which is the manifestation of the phase-charge duality well known in the theory of Josephson junctions.

The probabilities $P_{j}(\omega)$ also satisfy detailed balance relations (87), which allows us to express the total heat flux (106) in the Landauer form (2) with the transmission probability given by the sum of elastic and inelastic contributions,

[TABLE]

The elastic part $\tau^{\rm el}(\omega)$ is given by Eq. (21) with $Z_{J}^{-1}(\omega)$ defined in Eq. (101), and the inelastic one reads

[TABLE]

As in the previous sections, we consider two examples. First, we assume that the junction is coupled to the two resonators with the impedances (15) and that the coupling is sufficiently weak, $|{\cal R}_{j}(\omega)|\lesssim R_{q}$. In this case, one can derive the following approximations for the functions $P_{j}(\omega)$ and the junction impedance (101):

[TABLE]

Here, we have introduced the effective total impedance of the environment seen by the junction

[TABLE]

and used the approximation $F_{j}(t)=\Gamma_{j,T}|t|$ for the functions (103). The rates of phase diffusion are

[TABLE]

Equation (110) contains the total phase diffusion rate

[TABLE]

where the effective temperature of the junction is

[TABLE]

The approximation (110) is valid if

[TABLE]

In Fig. 13, we show the dependence of the heat flow (106) on magnetic flux for a certain choice of system parameters. For these parameters the effect of the junction impedance (101) on the transmission probability (21) is insignificant, and the cosine flux dependence of the heat power predominantly comes from the pre-factor $E_{J}^{2}\propto\cos^{2}(\pi\Phi/\Phi_{0})$ in front of the inelastic term (107). We have chosen an asymmetric coupling, that is why the heat flux (107) exhibits weak rectification. At low temperatures $k_{B}T_{1,2}\ll\hbar\omega_{1,2}$, the inelastic contribution to the heat flux (107) can be estimated as

[TABLE]

Figure 13 also shows that magnetic field independent background heat flux has strong temperature dependence. This background comes from the heat transport between two capacitively coupled resonators (28,30).

Our second example is again the case of Ohmic dissipation. In this case, one can derive analytic expressions for the photon absorption probabilities IG $P_{j}(\omega)$ and the effective impedance of the junction (101),

[TABLE]

where $\alpha_{\Sigma}=\alpha_{1}+\alpha_{2}$. The expression (117) is valid provided $\pi k_{B}T_{1,2}\ll\alpha_{\Sigma}E_{C}$ and the approximation for the impedance (118) requires the temperatures to be in the interval $E_{J}\ll k_{B}T_{1,2}\ll\alpha_{\Sigma}E_{C}/\pi$. Within this model one can analytically derive the correction to the thermal conductance. Assuming that $\alpha_{\Sigma}\gtrsim 8$, we find

[TABLE]

The first term in this expression comes from the heat current at $E_{J}=0$ given by Eq. (38), and the second term provides a negative correction to it. This correction is a combined effect of the inelastic contribution (107), which gives positive correction to $\kappa$, and of the junction impedance $Z_{J}(\omega)$, which suppresses the transmission probability together with the thermal conductance.

VI Summary

We have studied the photonic heat transport across a Josephson junction coupled to the two linear electric circuits, acting as thermal baths, which are characterized by the impedances $Z_{1}(\omega),$ $Z_{2}(\omega)$ and by temperatures $T_{1}$, $T_{2}$. We have shown that linear approximation, in which the nonlinear junction is replaced by an inductor, provides a rather good estimate of the heat flux between the thermal baths for any relation between the parameters $E_{J},E_{C},k_{B}T_{1,2}$, for any coupling strength and for any frequency dependence of the impedances $Z_{1,2}(\omega)$. This approximation fully accounts for hybridization between the modes of the junction and of thermal baths.

However, simple harmonic approximation cannot capture subtle effects such as thermal rectification, for example. Therefore we have developed more elaborate approximations, which take into account the nonlinear nature of the Josephson junction. In the weak-coupling limit, the photon transmission probability of the system is determined by the Bloch band structure of the junction energy spectrum. At low frequencies it is predominantly determined by the transitions between the junction states within one energy band, while at high frequencies – by the transitions between different bands.

In the limit $E_{J}\gg E_{C}$ and at weak coupling the photon transmission probability is given by a series of narrow peaks with the positions corresponding to the splitting between the energy levels in the Josephson potential well. At stronger coupling between the junction and the environment, these peaks overlap forming a single broad peak. An additional peak associated with the intraband transitions is formed at low frequencies. Its shape is determined by the impedances $Z_{1,2}(\omega)$, see Eq. (88). The heat flux at strong coupling is given by the sum of harmonic contribution (2) and the contribution of phase slips (86), with the latter dominating at low temperatures.

In the opposite limit, $E_{J}\ll E_{C}$, the Bloch bands become wide and the energy gaps between them almost vanish. In this case, the dominating contribution to the heat flux comes from the Landauer formula (2), in which one should put $I_{C}=0$. In addition, there exists a small correction depending on the Josephson energy $E_{J}$. If $E_{J}\ll k_{B}T_{1,2}$ this correction can be split into two parts: the elastic one, coming from the junction impedance averaged over phase fluctuations (101), and the inelastic correction (107), which is associated with the absorption and re-emission of the photons by the junction.

Many of our predictions can be experimentally tested. Indeed, with reasonable values of the system parameters we have obtained the heat fluxes in the range 10 - 100 aW (0.03%-0.3% of the maximum value corresponding to a single thermal conductance quantum), see, e.g., Figs. 4 and 10b. Such heat fluxes can be reliably detectedAlberto . This opens up a possibility for quantum thermodynamics experiments with tunable Josephson junctions.

VII Acknowledgement

We acknowledge useful discussions with Keiji Saito, Erik Aurell and Bayan Karimi. This work was supported by the Academy of Finland Centre of Excellence program (project 312057) and the European Union’s Horizon 2020 research and innovation programme under the European Research Council (ERC) programme (grant agreement 742559).

Appendix A Heat transport and linearized Langevin equation

In this appendix, we will demonstrate how one can derive the Landauer formula for the heat current (2) by solving the Langevin equations, which exactly describe quantum dynamics of a linear system Schmid ; Dhar . Here, we follow Nyquist Nyquist and Pascal, Courtois, and Hekking Hekking .

Kirchhoff’s equations for the circuit of Fig. 14 read

[TABLE]

Here the admittances of the environments 1 and 2 in the time domain are defined as $Y_{j}(t)=\int d\omega\,e^{-i\omega t}/2\pi Z_{j}(\omega)$, the voltage drop $V$ across the junction is related to the phase by means of the Josephson relation $V=\hbar\dot{\varphi}/2e$, the noises $\xi_{1}$ and $\xi_{2}$ are the Gaussian stochastic processes fully characterized by their pair correlators

[TABLE]

[TABLE]

The currents $I_{1}(t)$ and $I_{2}(t)$ flow through the impedances $Z_{1}(\omega)$ and $Z_{2}(\omega)$ respectively. The heat flux in this framework is given by the Joule heating released in the circuit 1 and averaged over the noises $\xi_{j}$,

[TABLE]

Equations (120-122) can be readily solved by means of Fourier transformation. We find the Fourier component of the voltage, $V_{\omega}=\int dt\,e^{i\omega t}\,V(t)$ in the form

[TABLE]

The heat flux (125) can now be transformed as

[TABLE]

Combining this expression with Eq. (126), we obtain

[TABLE]

Using the expressions (124) for the spectral power of the noises, one can show that this expression is equivalent to the Landauer formula (2) with the transmission probability having the form (21).

Appendix B WKB approximation in the limit $E_{J}\gtrsim E_{C}$

Let us consider the junction with $E_{J}\gtrsim E_{C}$. It is well known that in this regime the eigen-energies of the $2\pi$-periodic junction Hamiltonian $\hat{H}_{J}$ form Bloch bands with the cosine dispersion (65). The approximation (65) applies if $E_{n}<2E_{J}$, i.e., for the wave functions localized inside the potential wells and having eigenenergies smaller than the height of the barrier separating neighboring wells. In this appendix, we summarize approximate analytical expressions for both energy levels $E_{n}$ and half-bandwidth $\delta_{n}$ of the corresponding Bloch bands, which can be derived by means of the Wentzel-Kramers-Brillouin (WKB) approximation.

The energies of the quasiclassical levels $E_{n}^{\rm qcl}$ are determined by the Bohr-Sommerfeld quantization rule,

[TABLE]

where $\varphi_{0}(E_{n}^{\rm qcl})=\arccos\left(1-{E_{n}^{\rm qcl}}/{E_{J}}\right)$ are the classical turning points in the Josephson potential. Evaluating the integral in Eq. (128), one arrives at the equation for $E_{n}^{\rm qcl}$ in the form

[TABLE]

where ${\rm K}(k)$ and ${\rm E}(k)$ are the complete elliptic integrals of the first and the second kind (here, we use the definitions of these functions given in the book GR , which differ from the definitions used in the popular Mathematica package), and $n$ is the non-negative integer number taking the values $0,1,2,\dots$. Equation (129) has solutions provided $E_{J}/E_{C}>\pi^{2}/32$, while for lower values of the ratio $E_{J}/E_{C}$ no discrete levels exists in the potential well. The total number of the energy levels in the well equals to $n_{\max}+1$, where $n_{\max}$ is given by

[TABLE]

Here the brackets $\lfloor\dots\rfloor$ imply the floor function. Having found the quasiclassical energy levels $E_{n}^{\rm qcl}$ from Eq. (129), we correct them as follows

[TABLE]

We have numerically verified that the accuracy of the approximation (131) is 4% for $E_{J}/E_{C}=3$, 0.3% for $E_{J}/E_{C}=50$, and even better than that for higher values of the ratio $E_{J}/E_{C}$. The maximum error occurs for the top most level $E_{n_{\max}}$, while the positions of the low lying levels are very accurate. Next, we expand Eq. (129) at small energies $E_{n}^{\rm qcl}\ll 2E_{J}$ and keep only the two leading terms in the small parameter $E_{C}/E_{J}$. In this way, we obtain an approximation for the low lying energy levels,

[TABLE]

It agrees with the well known result of perturbation theory in weak anharmonicity of the Josephson potential Koch .

WKB approximation also allows one to find the half-bandwidths $\delta_{n}$. For the lowest band with $n=0$, the result is well known Averin ; SZ ,

[TABLE]

In order to find $\delta_{n}$ for $1\leq n\leq n_{\max}$, we use the WKB formula for the level splitting in a double-well potential Grag ,

[TABLE]

Here, we have introduced the frequency of classical oscillations at the bottom the potential well for a particle with the energy $E_{n}$,

[TABLE]

Evaluating the integral in Eq. (134), we arrive at

[TABLE]

Comparison with the exact numerical simulation shows that the approximation (136) is quite accurate. It has the accuracy 10% or better for all bands with energies below the barrier top and for ratios $E_{J}/E_{C}>3$. The maximum relative error in $\delta_{n}$ again occurs for the highest Bloch band inside the potential well, namely for $\delta_{n_{\max}}$.

Now, we turn to the matrix elements of the phase operator (139). We use the well known quasiclassical approximation for the matrix elements LL ; Nikitin ,

[TABLE]

Here $\varepsilon_{mn}^{+}=(E_{n}^{\rm qcl}+E_{m}^{\rm qcl})/2$, the frequency $\omega(\varepsilon_{mn}^{+})$ is given by Eq. (135) with $E_{n}$ replaced by $\varepsilon_{mn}^{+}$, and

[TABLE]

is the solution of the classical equation of motion for a particle in the Josephson potential well having the energy $\varepsilon_{mn}^{+}$. It is expressed in terms of Jacobi amplitude function ${\rm am}(u,k)$. The integral (137) can be solved analytically using the properties of Jacobi functionsGR . We find $\varphi_{mn}=0$ if $|m-n|$ is an even number, and

[TABLE]

if $|m-n|$ is an odd number. The parameter $q_{mn}$ is known as the so called nome in the theory of Jacobi functions,

[TABLE]

For completeness, we also provide the quasiclassical matrix elements of the momentum (44) for the odd $|m-n|$ (for even $|m-n|\not=0$ they are also equal to zero),

[TABLE]

For the low lying levels with $E_{n}\ll 2E_{J}$ one can derive simpler expressions. In this limit, one can use an approximation $E_{J}(1-\cos\varphi)\approx E_{J}(\varphi^{2}/2-\varphi^{4}/24)$ and treat the term $\propto\varphi^{4}$ as a perturbation. Keeping the two lowest order terms, we find

[TABLE]

For comparison, we also expand the quasiclassical matrix element (139) in the same limit and find

[TABLE]

The two expressions agree quite well even for $n=1$, which confirms the accuracy of the approximation (139). Numerics shows that the approximate expression for the matrix element of momentum $p_{n,n-1}$, given by Eq. (141), is very accurate. Namely, for $E_{J}/E_{C}=3$, we find a tiny error of 0.3%, while for $E_{J}/E_{C}=50$ the maximum error, occurring at $n=n_{\max}=5$, turns out to be 2%.

The matrix element between the states $n$ and $n\pm 3$, derived by perturbation theory in the quartic term, reads

[TABLE]

Expanding the quasiclassical expression (139) in powers of the small parameter $E_{C}/E_{J}$, we find

[TABLE]

Again, we note that the two expressions (144) and (145) agree quite well. Numerically, we find that for $E_{J}/E_{C}=50$, in which case $n_{\max}=5$, the maximum relative error for the matrix element $p_{n,n-3}$, given by Eq. (141), equals to 7% and occurs at $n=3$.

Appendix C Solution of the master equation for a harmonic oscillator

In this Appendix, we find the stationary solution of the master equation (67) in harmonic approximation, and use it to derive the expression (40) for the heat flux. If one replaces the nonlinear Josephson potential by a harmonic one, $E_{J}(1-\cos\varphi)\to E_{J}\varphi^{2}/2$, the energy levels become equidistant, $E_{n}=\hbar\omega_{J}(n+1/2)$, and independent of $q$. In this approximation, one finds $\omega_{n+1,n}=\omega_{J}$ for all $n$. The phase matrix elements connecting neighboring levels are given by

[TABLE]

All other matrix elements vanish. Accordingly, the transition rates between the levels (54) take the form

[TABLE]

where we have defined

[TABLE]

In this approximation, the solution of Eq. (67) can be found analytically,

[TABLE]

Substituting this result in the general expression for the heat flux (68), one arrives at the formula (40).

If we allow anharmonicity, but consider only the transitions between the neighboring levels, the occupation probabilities for an N-level system are

[TABLE]

where $m=1,2,...,N-1$, and

[TABLE]

As $N\rightarrow\infty$ and with vanishing anharmonicity, (150) reduces to (149).

Appendix D Derivation of Eq. (86)

In this appendix, we provide the details of the derivation of the expression (86) for the intraband contribution to the heat flux. We ignore the interband transitions and put $\delta_{n}=0$. In this approximation, the Hamiltonian of the junction, $\hat{H}_{J}$, drops out from the full Hamiltonian (3). The operator of the Josephson phase can be expressed as $\hat{\varphi}=-2ie(\partial/\partial q)$, hence the interaction Hamiltonians may be written in the form

[TABLE]

Since after making these approximations the Hamiltonian has become quadratic in all operators, one can use Wick’s theorem and the Baker-Campbell-Hausdorff formula for the commutators. This leads to

[TABLE]

It is straightforward to show by solving equations of motion for the quantum operators, that for the Hamiltonian of the form $\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{1}^{\rm int}+\hat{H}_{2}^{\rm int}$, with the interaction terms given by Eq. (151), the charge correlators read

[TABLE]

where the functions ${\cal F}_{j}(t)$ and ${\cal K}_{j}(t)$ are given by the integrals (84) and (85). Now one can find the probabilities ${\cal P}_{j}(\omega)$ from Eq. (83), and ${\cal P}(\omega)$ from Eq. (82).

Next, we substitute the correlation function (80) in Eq. (14) and find the intraband contribution to the heat flux in the form

[TABLE]

We can further transform $J^{\rm b}$ by applying a useful property of the function ${\cal P}_{1}(\omega)$, which is derived by applying the Fourier transformation to both sides of the identity

[TABLE]

Since the functions ${\cal F}$ and ${\cal K}$ are given by the integrals (84) and (85), we arrive at the result

[TABLE]

With the aid of this identity one can easily transform the heat flux (154) to the form (86).

Appendix E Derivation of Eq. (107)

In this appendix, we derive the expression (100) for the phase-phase correlation function by means of perturbation theory in small Josephson energy $E_{J}$. Considering the term $-E_{J}\cos\varphi$ in the full Hamiltonian (3) as a perturbation and keeping the terms up to $E_{J}^{2}$ in the expansion, one can express the time dependent phase operator $\hat{\varphi}(t)=e^{i\hat{H}t/\hbar}\hat{\varphi}e^{-i\hat{H}t/\hbar}$ in the form

[TABLE]

Here the time evolution of the operator $\hat{\varphi}_{0}(t)=e^{i\hat{H}_{0}t/\hbar}\hat{\varphi}e^{-i\hat{H}_{0}t/\hbar}$ is determined by quadratic Hamiltonian $\hat{H}_{0}=\lim_{E_{J}\to 0}\hat{H}$. As in the previous appendix, this ensures the validity of the Wick’s theorem for various products involving phase operators. In particular, with the aid of the Wick’s theorem one can prove that

[TABLE]

[TABLE]

Here, we have also used the property typical for quadratic Hamiltonians, namely, we have used the fact that the commutator \big{[}\hat{\varphi}_{0}(t),\hat{\varphi}_{0}(t^{\prime})\big{]} is proportional to the identity operator $\hat{E}$ and commutes with all other operators,

[TABLE]

From the equation of motion for $\hat{\varphi}_{0}(t)$, or applying path integral techniques Grabert , one finds $K(t)=K_{1}(t)+K_{2}(t)$, with the functions $K_{1,2}(t)$ having the form (104).

Using the expansion (157), transforming the commutators in it according to the rules (158) and (159), and taking the long-time limit, we express the phase-phase correlation function in the form

[TABLE]

Here, we have omitted the terms containing $\sin[\hat{\varphi}_{0}(t^{\prime})+\hat{\varphi}_{0}(t^{\prime\prime})]$, which vanish upon averaging because at $E_{J}=0$ the phase fluctuations are unrestricted. Next, we apply Wick’s theorem once again and find

[TABLE]

where $F(t)=F_{1}(t)+F_{2}(t)$ and the functions $F_{1,2}(t)$ are given by Eqs. (104). Similarly,

[TABLE]

Substituting these expressions in Eq. (161), and taking the Fourier transformation over the time difference $t_{1}-t_{2}$, we find

[TABLE]

The junction impedance $Z_{J}(\omega)$, appearing here, is defined by Eq. (101), the function $K_{\omega}$ reads

[TABLE]

the Fourier transformed phase-phase correlation function evaluated at $E_{J}=0$ is given by

[TABLE]

and the correlation function of $\sin\hat{\varphi}_{0}$ is defined as

[TABLE]

One can straightforwardly show that $S_{\sin\varphi}(\omega)=\pi P(\omega)$, where $P(\omega)$ is the photon emission probability (102). Comparing Eqs. (164) and (100) one can verify that they coincide in the lowest nonvanishing order of the perturbation theory $\sim E_{J}^{2}$, but Eq. (100) has more compact and physically transparent form.

The inelastic contribution to the heat flux originates from the last term in the correlation function (164) containing $S_{\sin\varphi}(\omega)$. Substituting this term in the general expression (9), we obtain

[TABLE]

One can transform this expression to a more physically meaningful form (107) invoking the property analogous to (156),

[TABLE]