Thermally pumped on-chip maser
George Thomas Azat Gubaydullin Dmitry S. Golubev Jukka P. Pekola
cond-mat.supr-con quant-ph










^1Pico group, QTF Centre of Excellence, Department of Applied Physics, Aalto University, P.O. Box 15100, FI-00076 Aalto, Finland 
^2Moscow Institute of Physics and Technology, 141700 Dolgoprudny, Russia


We present a theoretical model of an on-chip three-level maser in a superconducting circuit based on a single artificial atom 
and pumped by a temperature gradient between thermal baths coupled to different interlevel transitions. 
We show that maser powers of the order of a few femtowatts, well exceeding the resolution of the sensitive bolometry,
can be achieved with typical circuit parameters. 
We also demonstrate that population inversion in the artificial atom can be detected without measuring 
coherent radiation output of the maser. For that purpose, the system should operate as a three-terminal heat transport device. 
The hallmark of population inversion is the influx of heat power into the weakly coupled output terminal even though 
its temperature exceeds the temperatures of the two other terminals.
The proposed method of on-chip conversion of heat into microwave radiation and  control of energy-level populations  
by heating  provide additional useful tools for circuit quantum electrodynamics experiments.


Thermally pumped  on-chip maser
    Jukka P. Pekola^1,2
    December 30, 2023
===============================





§ INTRODUCTION


In recent years, considerable progress has been achieved in  both 
circuit quantum electrodynamics (cQED)<cit.> and circuit quantum thermodynamics <cit.>.
The experiments with superconducting quantum systems have led to spectacular fundamental studies of dynamical Casimir effect <cit.>, 
Lamb shift <cit.>, etc.
Moreover, in cQED devices lasing and other quantum optical phenomena <cit.> can be more easily exploited 
in the strong coupling regime <cit.>. 
For example, single atom lasing, which requires strong coupling, and  
which has been  demonstrated in experiments with a cavity containing 
a single rubidium <cit.> or cesium <cit.>  atom, 
has also been achieved in cQED experiments with a microwave resonator coupled to a superconducting artificial atom <cit.>.
Lasing has also been observed in an alternative cQED setup with a voltage biased Josephson junction and a resonator <cit.>.
Microwave masers based on superconducting elements and pumped by microwave sources have been theoretically analyzed, e.g., 
in Refs. <cit.>, and those based on semiconducting quantum dots — in Refs. <cit.>.
Furthermore, Josephson cascade micromaser has also been proposed <cit.>. 



The implementation of a heat driven maser, which allows conversion
of heat into coherent radiation,  has not yet been discussed in the cQED context.
At the same time, significant progress has been recently achieved  in
integrating the tools of ultrasensitive nanoscale bolometry into cQED devices <cit.>, 
and a unique platform for  studying the heat transport in the quantum limit has been created.  
Motivated by this progress, here we develop a theory and propose an experimental realization of a thermally driven single atom maser based on 
superconducting components.
We follow the original proposal by  Scovil and  Schulz-DuBois <cit.>, who have considered
a three level atom (qutrit), in which each interlevel transition is coupled to a separate thermal reservoir with
tunable temperature. Creating proper temperature differences between the reservoirs, one can induce population
inversion between the two energy levels coupled to the output port and in this way achieve lasing.
This device is an example of a quantum heat engine with the efficiency bound by the Carnot limit <cit.>.
We propose to use a superconducting loop with three Josephson junctions as an artificial atom with its three lowest
levels labeled as 0, 1, and 2. 
The artificial atom is capacitively coupled to the three high quality factor resonators which are tuned in resonance with
01, 12, and 02 transitions, see Fig. <ref>a. The hot (02) and cold (12) resonators are terminated by
resistors, which can be heated by bias currents, and the 01 resonator is connected to
the output port. The temperatures of the resistors can be monitored by normal metal - superconductor tunnel junction
thermometers <cit.>. After proper calibration, the small variations of the resistor temperatures 
can be converted into power dissipated or emitted by them.
Our theory is based on the usual Lindblad equation formalism, which has been developed  
for a general laser setup by Mu and Savage <cit.>, and has been later applied to the specific case
of thermally driven masers by Scully et al <cit.> and by Li et al  <cit.>.
We derive an explicit expression for the maser power in terms of circuit parameters and
show that it can reach up to few femtowatts for typical experimental conditions. 
This power is  several orders of magnitude higher than the lowest power detectable in experiment <cit.>.
We further show that one can detect the population inversion
between the states 0 and 1 without measuring the coherent radiation output of the device. 
Instead, one can connect a resistor to the output resonator and monitor how  the power dissipated in it varies with its temperature.
In presence of the population inversion the  power keeps flowing into the output resistor even if its temperature exceeds
the temperatures of the other two resistors. One can  understand this feature by assigning a negative effective temperature 
to the artificial atom with population inversion. In the heat transport context,  
negative temperature is hotter than any positive temperature, that is why the heat always flows into the resistor connected to the output resonator.












§ MODEL

 
In our model, the artificial atom is realized as a superconducting loop 
with three Josephson junctions, similar to the conventional persistent-current
flux qubit <cit.>, as shown in Fig. <ref> (b). 
The loop is formed by the two bridges connecting a superconducting island, which is
restricted by the capacitors C_c,C_h,C_d and the junctions 1 and 3 and shown by dotted line in Fig. <ref> (b),
to the superconducting ground electrode. The left bridge contains only one Josephson junction
(junction 3), while the right bridge – two other junctions (junctions 1 and 2).
The capacitors C_c,C_h,C_d couple the island
to the cold, hot, and output (or detector) resonators respectively. The latter can be realized as
usual coplanar waveguide λ/2-resonators widely used in cQED circuits <cit.>.
The fundamental modes of the resonators have the frequencies ω_r, where the index r can be h,c or d.
We denote the characteristic impedances of the resonators as Z_0r.
On the other side, the cold and the hot resonators are coupled to the resistors R_c
and R_h via the capacitors C'_c and C'_h.  
Finally, the capacitor C'_d couples  
the output resonator to a resistor R_d or to a transmission line with the impedance Z_0. 
The latter should guide the output radiation to the spectrum analyser. 
The temperatures T_r of all resistors can be varied by applying DC bias currents to them without influencing the qutrit,
and can be monitored with normal metal - superconductor tunnel junction thermometers <cit.>. 
The maser is pumped by heating up the hot resistor R_h.
The first three levels of the artificial atom have the energies E_0, E_1 and E_2, and the corresponding transition frequencies are
ω_ij=(E_i-E_j)/ħ (i,j= 0,1,2).
These frequencies can be controlled by the magnetic flux Φ applied to the loop, and may be slightly detuned from the
frequencies of the corresponding resonators, so that ω_h=ω_20+Δ_h, ω_c=ω_21+Δ_c,
and ω_d=ω_10+Δ_d. 
The asymmetry of the loop with three junctions makes the transitions between the states 0 and 2 possible.
Such transitions are required for the pumping of the maser.
The quality factors of the fundamental modes of the resonators <cit.>, Q_r=π/(2ω_r^2Z_0rR_r C_r^' 2), 
should be sufficiently high to avoid their overlap with unintended transition frequencies,  Q_r≫ω_r/|ω_10-ω_21|.
However, Q_r should not be too high to allow some tolerance in the detunings Δ_r. 

 
One can understand the operation principle of the maser first assuming  that the qutrit is 
uncoupled  from the output resonator and coupled only to the hot and the cold ones. The hot resistor induces
the transitions from the state 0 to the state 2 with the rate Γ_20, and the inverse transitions 2→ 0 with the rate Γ_02.  
Likewise, the cold resistor is responsible for the transitions 1→ 2 with the rate Γ_21 and 2→ 1 with the rate Γ_12.
These rates satisfy the detailed balance relations Γ_02/Γ_20=exp(ħω_20/k_BT_h) 
and Γ_12/Γ_21=exp(ħω_21/k_BT_c) . 
One can easily verify that in this case the population inversion between the  
states 0 and 1 is achieved if <cit.> 

    ħω_21/k_BT_c≥ħω_20/k_BT_h.
 
Once the population inversion is there, one can weakly couple the output resonator to the qutrit.
The latter will be driven to a coherent state and will emit coherent radiation to the output transmission line. 



 


We now turn to the rigorous theoretical framework,
which allows one to evaluate maser power for given circuit parameters.  
We assume that all three junctions in the loop have the  same critical current I_C.
In this case, the Hamiltonian of the loop reads

    H_J    =    -4E_C∂_φ^2 
    +E_J[3cos(φ_x/3)
    
          
    - cosφ- 2cos((φ-φ_x)/2)],

where φ is the Josephson phase difference between the central superconducting island of the device
and the ground, E_C is the charging energy of the island, φ_x=2πΦ/Φ_0  
is the phase shift induced by the magnetic flux, Φ_0 is the magnetic flux quantum and E_J=ħ I_C/2e is the Josephson energy
of a single junction.
Here we follow Ref. <cit.> 
and assume that the phases of the junctions 1 and 2 are equal, φ_1=φ_2=φ/2.
The dependence of the potential on the phase difference θ=φ_1-φ_2 may be
ignored if E_J≫ E_C, the self-capacitance of the island between the junctions is small, and one considers only
few low lying energy levels in the potential well. In this case, motion in θ-direction should be frozen
because the resonators are coupled only to the phase φ.
Expanding the potential energy near the minimum φ_x/3, 
and keeping only the third and the fourth order terms in the phase difference φ̅=φ-φ_x/3,
which is justified provided 3E_Jcos(φ_x/3)/2≫ E_C, one obtains
the energy level spacings as 
ħω_10=√(12 E_J E_Ccos(φ_x/3))-3E_C/4
and ħω_21=ħω_10-3 E_C/4 (see Appendix for details). 
The potential energy of the loop, numerically evaluated energy levels of the Hamiltonian (<ref>) and the square of the corresponding wave functions  
are shown in Fig. <ref>(c). 
  
Considering only the three lowest levels of the qutrit and high quality factor resonators, one can express the Hamiltonian of the system in the form

    H = ∑_n=0^2 E_n|n⟩⟨ n|+ ∑_r=h,c,dħω_r(a^†_r a_r +1/2)
    
          
    + ħ g_d[(a^†_d +a_d)(σ^+_01  + σ^-_01)]
    + ħ g_h[(a^†_h +a_h)(σ^+_02  + σ^-_02)]
    
          
    + ħ g_c[(a^†_c +a_c)(σ^+_12  + σ^-_12)].

Here g_r are the coupling constants between the qutrit and the resonators  given by 

    g_mn=ω_r^2C_r/2e√(ħ Z_0r/π)φ_mn,

with g_h=g_20, g_c=g_21 and g_d=g_10 (for explicit approximate expressions see Eqs. (<ref>)),   
φ_ij=⟨ i| φ̂|j⟩ 
are the matrix elements of the Josephson phase between the states of the qutrit,  
a_r and a^†_r are the ladder operators of the resonators, and
σ_ij^+=|i ⟩⟨ j|
and σ^-_ij=|j⟩⟨ i| are the transition operators between the levels of the qutrit.
The damping rates of the fundamental modes of the resonators, κ_r=(2/π)ω_r^3 Z_0rR_rC_r^' 2, 
are induced by their coupling to the resistors R_r, or, in case of the output resonator – to the outer transmission line with the impedance Z_0,
in that case, one should replace R_d by Z_0.
Below we assume κ_h/g_h≳ 1 and κ_c/g_c≳ 1 and ignore 
coherent coupling between the qutrit and the hot and cold resonators, 
treating the latter as parts of the environment inducing the transitions between the qutrit levels. 
However, at this stage we keep the ratio κ_d/g_d arbitrary. In this case, 
the  Lindblad equation governing the time evolution of the density matrix ρ of the "qutrit plus output resonator" system has the form

    d ρ/dt   =   -i/ħ[H_q-d,ρ]+κ_d (N_d^ eq+1)(aρ a^†-1/2{a^†a,ρ})
    
          
     + κ_d N_d^ eq(a^†ρ a-1/2{a a^†,ρ})
    
          
    + ∑_i=1^2∑_j=0^i-1Γ_ij(σ_ij^-ρσ_ij^+
    -1/2{σ_ij^+σ_ij^-,ρ})
    
          
    + ∑_j=1^2∑_i=0^j-1Γ_ij(σ_ij^+ρσ_ij^--1/2{σ_ij^-σ_ij^+,ρ}).

Here  H_q-d is the Hamiltonian describing "qutrit plus output resonator" system having the form (<ref>) with g_h=g_c=0, and r=d,
N_d^ eq=1/(e^ħω_d/k_BT_d-1) is the equilibrium photon number in the output resonator, 

    Γ_ij   =    |φ_ji|^2 ħω_ji/2 e^2  Re[1+n_h(ω_ji)/Z_h(ω_ji)+1+n_c(ω_ji)/Z_c(ω_ji)],  (i<j)
    
    Γ_ij   =    |φ_ij|^2 ħω_ij/2 e^2  Re[n_h(ω_ij)/Z_h(ω_ij)+n_c(ω_ij)/Z_c(ω_ij)],  (i>j)

are the transition rates from the state j into the state i,  
Z_h(ω),Z_c(ω) are the impedances of the hot and cold resonators (<ref>), 
and n_h(ω) and n_c(ω) are the effective photon distributions 
in these resonators. According to Eqs. (<ref>) both hot and cold resonators
contribute to all transitions, thus the expressions (<ref>) are valid
even for low quality resonators with overlapping spectrum lines.
In the limit Γ_02≪κ_h and Γ_12≪κ_c the resonators
are almost decoupled from the qutrit and n_h,c(ω)=1/(e^ħω/k_BT_h,c-1).
For Γ_02≳κ_h, Γ_12≳κ_c the distributions n_h,c(ω) should be found self consistently 
from the equation of motion for the average values of the ladder operators ⟨ a_h,c⟩.  
For high quality factor resonators, with non-overlapping spectral lines, and in the semiclassical approximation, 
the equations for ⟨ a_h,c⟩ can be re-written
in terms of the average numbers of photons in the resonators N_h = ⟨ a^†_h a_h⟩≈ n_h(ω_20) and 
N_c=⟨ a^†_c a_c⟩≈ n_c(ω_21), see  Eqs. (<ref>,<ref>).
In the experimentally relevant limit 3E_Jcos(φ_x/3)/2≫ E_C 
the phase matrix elements can be approximately evaluated by means of the perturbation theory in φ̅. In this way, one finds
|φ_10|^2 = √(4E_C/3E_Jcos(φ_x/3)), |φ_21|^2 = 2|φ_10|^2 
and |φ_20|^2=(E_C/54E_J)tan^2(φ_x/3)/cos(φ_x/3)
(see Sec. <ref> of the Appendix for details). 
We note that the matrix element φ_20, needed for pumping the maser, 
differs from zero only in presence of magnetic field.
In the numerical simulations we compute the matrix elements exactly instead of using these approximate expressions.
We use the parameters such that the third energy level is close to the barrier top of the double-well potential profile, as shown in Fig. <ref> (c). This
increases the anharmonicity ω_10-ω_21 above the value 3E_C/4ħ, predicted by the perturbation theory in φ̅,
and helps to increase the maser power.

Instead of numerically solving the Lindblad equation (<ref>),
for most practical purposes one can use a simpler semiclassical model of the device 
valid for high quality factor resonators without overlap of the spectral lines 
and for high numbers of photons in the output resonator N_d=⟨ a^†_d a_d⟩≳ 1.
Similar models are often used in the literature,
see e.g. Refs. <cit.>.
In this approximation, the system dynamics is governed by master equation for the vector of
occupation probabilities of the qutrit states p=(p_0,p_1,p_2)^T,

    ṗ=Γ p

with the transition rates matrix

    Γ=[ [ -Γ_10-Γ_20       Γ_01       Γ_02;       Γ_10 -Γ_01-Γ_21       Γ_12;       Γ_20       Γ_21 -Γ_02-Γ_12;            ]],

and by the equations for the numbers of photons,

    Ṅ_d    =   Γ_01p_1 - Γ_10p_0 - κ_d (N_d - N_d^ eq),
    
    Ṅ_h    =   Γ_02p_2 - Γ_20p_0 - κ_h (N_h - N_h^ eq),
    
    Ṅ_c    =   Γ_12p_2 - Γ_21p_1 - κ_c (N_c - N_c^ eq).

Here N_r^ eq=1/(e^ħω_r/k_BT_r-1) are the equilibrium photon numbers.
Since the resonator lines do not overlap, the transition rates (<ref>) take the simple form

    Γ_01 = γ_d (1+N_d),  Γ_10 =γ_d N_d,
    
    Γ_02 = γ_h (1+N_h),  Γ_20 =γ_h N_h,
    
    Γ_12 = γ_c (1+N_c),  Γ_21 =γ_c N_c,
 
where the rates of downward transitions in the absence of photons in the resonators are

    γ_h =  ħω_20/2 e^2  Re[|φ_20|^2/Z_h(ω_20)], γ_c =  ħω_21/2 e^2  Re[|φ_21|^2/Z_c(ω_21)],
    
          γ_d =  ħω_10/2 e^2  Re[|φ_10|^2/Z_d(ω_10)].

These rates can be simplified for small detunings Δ_r, 

    γ_r = g_r^2κ_r/Δ_r^2 + κ_r^2/4,   r=h,c,  or  d.

Eqs. (<ref>) and (<ref>-<ref>) can be easily solved numerically in the stationary case.
Once the photon numbers N_r are known, all other parameters can  be evaluated. 
For example, the output power is given by the  expression (see Appendix <ref> for details):  

    P_out/ħω_d = γ_hγ_cγ_d[(1+N_d)N_h(1+N_c)-N_d(1+N_h)N_c]/B,

where

    B = γ_dγ_h[(1+2N_h)(1+N_d)+(1+N_h)N_d] 
    
          
    + γ_dγ_c[(1+N_c)(1+N_d)+(1+2N_c)N_d] 
    
          
    + γ_hγ_c[N_h (1+N_c)+(1+N_h)N_c+N_hN_c].






§ RESULTS AND DISCUSSION


We numerically find the stationary solution of Eq. (<ref>) in the standard way <cit.>.
In Figs.  <ref>a- <ref>f 
we plot the resulting photon number distributions and Wigner functions of the output resonator  
for several values of the temperature bias and assuming typical circuit parameters given in the caption.
These parameters result in the following values of
the coupling constants and the damping rates for hot and cold resonators:
g_h/2 π≈45.34  MHz < k_h/2 π≈ 161.17 MHz  and
g_c/2 π≈ 22.68   MHz< k_c/2 π≈ 60.90 MHz. For the  output resonator,  we obtain 
g_d/2 π≈ 29.88   MHz≫ k_d/2 π≈ 58.21 KHz. We also obtain γ_h=52.25 MHz, γ_c=33.94 MHz,
and γ_d=52.88 MHz.
Although for these parameters, the ratios γ_h/κ_h and γ_c/κ_c are small, 
they are not small enough to assume equilibrium photon numbers in the hot and cold resonators. That is why  
we have used Eqs. (<ref>,<ref>) in oder to find the stationary values of N_h and N_c.
We observe that the photon number distribution evolves from the Boltzmann form e^-ħω_dn/k_BT_d for T_d=T_h=T_c (Fig. <ref>b) to
the Poissonian form P^ Pois_n=e^-⟨ n⟩⟨ n⟩^n/n! typical for a lasing state (Fig. <ref>f)
as the temperature bias becomes stronger. Simultaneously, the Wigner function evolves from a gaussian thermal distribution (Fig. <ref>a)
to a ring-shaped form indicating high amplitude harmonic oscillations in the coherent state (Fig. <ref>e). 

Following the recepie of Refs. <cit.>, 
we derive an approximate expression 
for the output power of the maser in the lasing regime N_d≫ 1, N_h≳ N_c and γ_d≫κ_d(1+3N_c).
Simplifying the exact expression (<ref>) in this limit we obtain

    P_out=γ_cγ_h(N_h-N_c)/γ_c(2+3 N_c)+γ_h(2+3 N_h)ħω_10.

Eq. (<ref>) rather  well agrees with the exact numerics for the chosen parameter values, see Fig. <ref> (g).
It provides a good approximation in the high temperature bias regime T_c/T_h≪ω_10/ω_20, in which lasing
is expected. In Fig. <ref> (g) we also plot the exact power values [Eq. (<ref>)] with red line, which in high 
temperature bias regime exactly matches with the power obtained by solving the Lindblad equation (<ref>), while
the green line represents high power regime  [Eq. (<ref>)] for a certain choice of parameters.
In  the plot for high power (green line), we use  g_h/2 π=191.73 MHz, g_c/2 π=258.16 MHz, g_d/2 π=74.71 MHz,
κ_h/2 π=577.86 MHz, κ_c/2 π=311.46 MHz, and κ_d/2 π=58.22 KHz.

For typical system parameters ω_10/2π≈ 5 GHz and E_C/2πħ≈500 MHz, T_h≈ 0.35 K, and
T_c≲ 0.05 K,
the maser operates in the regime N_h∼ 1 and N_c≪ 1.
For sufficiently strong coupling to the hot resonator, γ_h≫γ_c,
the power (<ref>) approaches the limiting value P̃_out=ħω_10γ_c/3. Increasing the rate γ_c
by design, one can push this power up. However, in our model γ_c cannot exceed κ_c,
otherwise the cold resonator cannot be efficiently cooled. The damping rate κ_c, in turn, is limited by the condition
κ_c≲ (ω_10-ω_21)/2 = 3E_C/8. Indeed, for higher κ_c  
the output resonator becomes coupled to the 1↔ 2 transition and gets cooled. 
Thus, at most, one can achieve γ_c∼ 3E_C/8 
and the maximum output power of the device can be estimated as P_max≈ω_10E_C/8.
For typical device parameters given above one finds P_max=1.3 fW. 
This value is comparable with the maser power of 0.7 fW reported in Ref. <cit.>, 
but it is significantly lower than 255 fW observed in a laser pumped by a voltage biased Josephson junction <cit.>. 
The linewidth of the maser can be estimated 
according to Schawlow-Townes formula <cit.>, δω=κ/ 2N_d. For maximum maser power one finds
N_d=P_max/ħω_dκ_d≈ E_C/8ħκ_d and 
the linewidth reduces to δω≈4ħκ_d^2/ E_C= 4ħω_d^2/E_C Q_d^2.
For the parameters quoted above and for Q_d∼ 10^3 we estimate δω/2π∼ 0.2 MHz.
The efficiency of the maser is defined as the ratio of the output maser power to the input heat absorbed
from the hot resistor, η=P_out/P_in. One can show <cit.> that in the lasing regime  
η≈ω_10/ ω_20. This value, however, cannot exceed the Carnot efficiency
1-T_c/T_h. Indeed, for ω_10/ ω_20 > 1-T_c/T_h, no population inversion in the qutrit
can be created (see Eq. (<ref>)) and, hence, no lasing can occur. Instead, in this regime, the trivial heat transfer from the hot resistor
to the output port is going on. The efficiency of the maser is plotted in Fig. <ref> (h).
Finally, the lasing threshold is determined by the condition (<ref>),
which can be satisfied only if κ_d≪γ_d.

 

Ideally, the output radiation of the maser should be detected with a spectrum analyser.
We will now demonstrate that one can replace the latter by the simplest broadband detector – a resistor R_d –
and still detect the population inversion in the
qutrit. Let us  assume that the coupling between all three resonators and the resistors is 
stronger than the coupling between them and the qutrit,  C_r'≫ C_r. No lasing can occur in this regime, and the photons in all
resonators are equilibrated with the corresponding resistors, N_r=N_r^ eq. 
The output power of the device (<ref>)
is proportional to the combination P_out∝ (1+N_d^ eq)N_h^ eq(1+N_c^ eq) - N_d^ eq(1+N_h^ eq)N_c^ eq. It
is the net contribution of the forward process, in which one photon comes from the hot resonator
and two photons are subsequently emitted by the qutrit into cold and output resonators, and of the reverse process.  
Substituting Bose functions in the above expression, one can check that  P_out>0  if and only if

    ħω_d/k_BT_d + ħω_c/k_BT_c > ħω_h/k_BT_h.

Equation<ref> ensures that the entropy production for the net forward process is positive.
As we have discussed above, for g_d≪ g_c,g_h the population inversion
in the qutrit appears if the condition (<ref>) is satisfied. Obviously, in this case
the inequality (<ref>) is satisfied, and hence P_out>0, for any temperature of the output resonator T_d
including high temperatures T_d>T_h. This unusual behaviour discussed by Geusic et al <cit.> is illustrated in the context 
of superconducting circuits in Fig. <ref> for a certain choice of the system parameters.
In Figs. <ref>a,c we plot the power dissipated in the resistor R_d and 
the ratio of the occupation probabilities  p_1/p_0 for a fixed 
temperature T_d=0.03 K as functions of the magnetic flux. 
In Figs. <ref>c and <ref>d, we take the value of the magnetic flux Φ/Φ_0=0.32 close to the maximum power, 
and  plot this power together with the ratio p_1/p_0 versus the temperature T_d.
P_out indeed remains positive even for T_d>T_h=0.38 K.
Thus, we have shown that in presence of the population inversion between the levels 0 and 1 heat can flow from the colder resistor 
R_h to the hotter resistor R_d. One can also interpret this observation as follows: population inversion
implies negative temperature of the qutrit, and the latter is formally hotter than any positive temperature.
Clearly, this property also holds in the lasing regime where the photon number distribution in the output resonator strongly 
deviates from equilibrium. 






§ CONCLUSION



In conclusion, we have proposed a model of thermally driven on-chip three level maser
in  a superconducting circuit,  which converts heat into coherent radiation.
We have also proposed an experimentally feasible method of detecting
the population inversion in the artificial atom by thermometry.
We have derived a simple analytical expression for the output power of the maser in terms of the circuit parameters. 
We have shown that powers up to few femtowatts can be achieved in a typical cQED setup.
We believe that the proposed way of converting heat produced by DC-biased resistors into coherent radiation on-chip 
can be very useful for quantum circuit applications. The same technique can be used to control and detect the populations of
energy levels of quantum superconducting devices.

Our model can be generalized to a heat accelerator and a heat switch setups if one allows 
manipulating the coupling between the qutrit and one of the heat baths. 
One can also include more than one qutrit into it
and consider more subtle effects like, for example, 
superradiance of several artificial atoms <cit.>.
Finally, our analysis can be generalized to a four level laser model,
which is, however, more challenging to realize in the experiment.




§ ACKNOWLEDGEMENT

 
This work was supported by the Academy of Finland Centre of Excellence
program (project No. 312057), European Union's Horizon 2020 research 
and innovation programme under the European Research Council (ERC) 
programme (grant agreement No. 742559) and Marie Sklodowska-Curie grant agreement No. 843706. 
We acknowledge the computational resources provided by the Aalto Science-IT project.








§ APPROXIMATE THEORY OF THE QUTRIT




We consider the setup shown in Fig. <ref>b and assume that all Josepshon junctions in the superconducting loop have the same critical current  I_C.
In this case, the phase difference on the junction 3 equals to φ and on the junctions 1 and 2 it is φ'/2.
We can assume that the phase drops across junction 1 and 2 are equal when 
 E_J≫ E_C 
and capacitance to the ground (C_0) by the island between the junctions 1 and 2 is C_0≪ C_J/4, where C_J is the capcitance of a junction <cit.>.
The total current through the loop is I= I_Csin(φ'/2)+I_Csinφ.
The flux quantization condition implies φ_3=φ '+2πΦ/Φ_0. This leads to

    I= I_C[sin(φ-φ_x)/2+sinφ],

where φ_x=2πΦ/Φ_0. The potential energy of the loop reads

    U=ħ/2e∫_φ_x/3^φ I_C(sin(φ_1-φ_x)/2+sinφ_1) dφ_1
    
     =E_J[3cos(φ_x/3)-cosφ-2cos((φ-φ_x)/2)].

Combining this expression with the charging energy of the island -4E_C∂^2_φ, we arrive at the
Hamiltonian of the artificial atom (<ref>).

The potential energy (<ref>) achieves the minimum value U=0
at φ=φ_x/3 provided -π≤φ_x≤π. 
Expanding it near the minimum 
up to the quartic order in φ̅=φ-φ_x/3, 
we approximate the Hamiltonian (<ref>) as 

    H=-4E_C∂^2/∂φ̅^2+Ẽ_J[φ̅^2/2
    -tanφ_x/3φ̅^3/12-φ̅^4/32],

where  Ẽ_J=3E_Jcos(φ_x/3)/2.
In the lowest order perturbation theory in ∝φ̅^3,φ̅^4 we obtain
the usual expression for the low lying energy levels of a quarctic oscillator
  
    E_n=√(8 Ẽ_J E_C)(n+1/2)-E_C/16(6n^2+6n+3),

which is valid provided E_C≪Ẽ_J. The correction ∝φ̅^3 does not affect
the energy levels due to symmetry.
 The transition frequencies between the three lowest levels are

    ħω_10   =    E_1-E_0=(8 Ẽ_J E_C)^1/2-3E_C/4,
    
    ħω_21   =    E_1-E_0 = (8 Ẽ_J E_C)^1/2-3E_C/2,
    
    ħω_20   =    E_2-E_0 =2(8 Ẽ_J E_C)^1/2-9E_C/4.




The corrected nth eigenstate is given by
|n⟩=|n⟩_2+|n⟩_3 +|n⟩_4, where |n⟩_2 is
the eigenstate of the unperturbed quadratic Hamiltonian,
the correction coming from ∝φ̅^3 term is denoted as |n⟩_3 and the
corrections due to ∝φ̅^4 term — as |n⟩_4.
These corrections read
 
    |n⟩_3   =   -1/24(4E_C/3 E_J)^1/4tan(φ_x/3)/(cos(φ_x/3))^1/4
          × [(√(n(n-1)(n-2))/3)|n-3⟩_2
          -(√((n+1)(n+2)(n+3))/3) |n+3⟩_2
          +3n√(n)|n-1⟩_2 
     -3(n+1)√((n+1))|n+1⟩_2],


    |n⟩_4    =    -1/64√(2)√(E_C/E_Jcosφ_x/3)[(4n-2)√(n(n-1))|n-2⟩_2 
          - (4n+6)√((n+1)(n+2)) |n+2⟩_2
          -(1/2)√((n+1)(n+2)(n+3)(n+4))|n+4⟩_2 
           +(1/2)√(n(n-1)(n-2)(n-3))|n-4⟩_2].

Having found the corrections to the eigenfunctions, we can evaluate
the phase matrix elements:

    ⟨ n-1| φ̂|n⟩
    = (4 E_C/3 E_Jcosφ_x/3)^1/4√(n)[ 1+n√(E_C/32E_J)],
    
          ⟨ n-2| φ̂|n⟩
    =√(E_C/3 E_Jcosφ_x/3)tanφ_x/3/6√( n(n-1)),
    
          ⟨ n-3| φ̂|n⟩
    = (4 E_C/3 E_Jcosφ_x/3)^3/4√(n(n-1)(n-2))/48.


With the explicit expressions for the matrix elements (<ref>), we find the approximate expressions for the
coupling constants (<ref>) between the qutrit and the resonators,

    g_h    =   ω_h^2C_h/2e√(ħ Z_0h/π)√(2E_C/3 E_Jcosφ_x/3)tanφ_x/3/6,
    
     
    g_c    =   √(2) ω_c^2C_c/2e√(ħ Z_0c/π)(4E_C/3E_Jcosφ_x/3)^1/4,
    
     
    g_d    =   ω_d^2C_d/2e√(ħ Z_0d/π)(4E_C/3E_Jcosφ_x/3)^1/4.




§ IMPEDANCES OF THE RESONATORS


The transition rates between the states of the qutrit (<ref>) are expressed in terms
of the impedances of the resonators. These impedances include the resonators themselves as well as
the coupling capacitors and the resistors, which terminate the resonators and serve as thermal baths
in our setup. The impedance of the resonator r (where r=h,c or d) reads

    Z_r(ω)   =   Z_0r(R_r + 1/-iω C'_r)cosω t_r -i Z_0rsinω t_r/Z_0rcosω t_r -i (R_r + 1/-iω C'_r)sinω t_r
          
    + 1/-i ω C_r.

Here t_r are the travel times of the microwaves along the resonators.









§ MASER OUTPUT POWER AND LASING THRESHOLD




Solving the system (<ref>), we find the steady state occupation probabilities of the qutrit states in the form

    p_0=Γ_01 +  Γ_02Γ_21/Γ_02+Γ_12/A,  
    p_1 = Γ_10 + Γ_12Γ_20/Γ_02+Γ_12/A,
    
          
     A=(1+Γ_20/Γ_02+Γ_12)(Γ_01 +  Γ_02Γ_21/Γ_02+Γ_12)
    
          
    + (1+Γ_21/Γ_02+Γ_12)(Γ_10 + Γ_12Γ_20/Γ_02+Γ_12),
    
          
    p_2 = 1-p_0-p_1.

The steady state numbers of photons in the resonators follow from Eqs. (<ref>),

    N_d    =    N_d^ eq+(Γ_01p_1-Γ_10p_0)/κ_d,
    
    
    N_h    =    N_h^ eq+(Γ_02p_2-Γ_20p_0)/κ_h,
    
    
    N_c    =    N_c^ eq+(Γ_12p_2-Γ_21p_1)/κ_c.

Since the rates Γ_ij (<ref>-<ref>) depend on the photon numbers N_r, 
these equations should be numerically solved for N_r. 

The output power of the maser reads

    P_out =  ħω_d(Γ_01p_1-Γ_10p_0).

With the aid of the explicit solutions (<ref>) and the expressions for the rates (<ref>-<ref>), 
this expression reduces to the form (<ref>) given in the main text.

Lasing threshold is defined as the value of the temperature bias at which N_d, which satisfies Eq. (<ref>), 
begins to exponentially grow in time provided it was small at the beginning. This  condition  can be formulated as

    d/dN_d[ Γ_01p_1 - Γ_10p_0 - κ_d (N_d - N_d^ eq) ]|_N_d→ 0>0.

This condition can be simplified in the limit γ_d≪γ_h,γ_c, which is  favorable for lasing.
In this limit we find

    N_h-N_c/N_h(1+N_c)+(1+N_h)N_c+N_hN_c > κ_d/γ_d.

In practical terms, since the combination in the left side of the equation does not exceed 1, one
should meet the conditions κ_d≪γ_d≪γ_h,γ_c while designing the sample.







§ WIGNER FUNCTION


The Wigner function of the output resonator
is evaluated as

    W(x',p')=∑_k,l f_kl(x',p')⟨ k|ρ^*_ res|l⟩,

where the sum runs over the states of the harmonic oscillator, 
the density matrix elements ⟨ k|ρ^*_ res|l⟩ are found by solving
the Lindblad equation (<ref>) and taking trace over all variables except the states of the output resonator,
and the function f_kl(x',p') is defined as <cit.>

    f_kl(x',p')=∫_-∞^∞dy/2πψ_k^*(x'-y/2)ψ_l(x'+y/2) e^-i y p'.

In Eqs. (<ref>,<ref>) both the coordinate and the momentum are normalized by their zero point fluctuation values.
Next, ψ_l(x) is the wave function of the harmonic oscillator corresponding to its energy level with the number l.
The integral (<ref>) can be solved exactly <cit.>

    f_kl(x',p')   =   (-1)^k/π√(k!/l!)(4H)^(l-k)e^-2He^i(l-k)arctanp'/x'
          ×
    L_l^(l-k)(4H),

where H=(x'^2+p'^2)/2 and L_l^k(x') are the generalized Laguerre polynomials. 





99

Wallraff  A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar, S. M. Girvin and R. J. Schoelkopf, 

https://doi.org/10.1038/nature02851 Nature 431, 162-167 (2004). 

Devoretap M. Devoret, S. Girvin, and R. Schoelkopf, 

https://doi.org/10.1002/andp.200710261 Ann. Phys. (Leipzig) 16, 767 (2007).  


pekola2015 J. P. Pekola, 

https://doi.org/10.1038/nphys3169 Nat. Phys. 11, 118 (2015).   

Hoffer2016P. P. Hofer, M. Perarnau-Llobet, J. B. Brask, R. Silva, M. Huber, and N. Brunner, 

https://link.aps.org/doi/10.1103/PhysRevB.94.235420
Phys. Rev. B 94, 235420 (2016).

vinjanampathy_quantum_2016 S. Vinjanampathy, J. Anders, 

https://doi.org/10.1080/00107514.2016.1201896 Contemp. Phys. 57, 545 (2016). 

Wilson C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori and P. Delsing,  

https://doi.org/10.1038/nature10561 Nature 479, 376 (2011).  

Fragner A. Fragner,M. Göppl, J. M. Fink, M. Baur, R. Bianchetti, P. J. Leek, A. Blais, A. Wallraff, 

https://doi.org/10.1126/science.1164482 Science 322, 1357 (2008).  

Silveri M. Silveri, S. Masuda, V. Sevriuk, K. Y. Tan, M. Jenei, E. Hyyppä, 
F. Hassler, M. Partanen, J. Goetz, R. E. Lake, L. Grönberg and M. Möttönen, 

https://doi.org/10.1038/s41567-019-0449-0 Nat. Phys. 15, 533 (2019). 


Haroche J. M. Raymond, M. Brune, and S. Haroche, 

https://doi.org/10.1103/RevModPhys.73.565 Rev. Mod. Phys. 73, 565 (2001) ; 
S. Haroche and J. M. Raymond, Exploring the Quantum: Atoms, Cavities and Photons (Oxford University Press, New York, 2006).

Meschede1985 D. Meschede, H. Walther,  and G.  Müller, 

https://link.aps.org/doi/10.1103/PhysRevLett.54.551Phys. Rev. Lett. 54, 551 (1985).

Kimble2003 J. McKeever, A. Boca, A. D. Boozer, J. R. Buck and H. J. Kimble, 

https://doi.org/10.1038/nature01974 Nature 425, 268 (2003).

Astafiev2007 O. Astafiev, K. Inomata, A. O. Niskanen, T. Yamamoto, Yu. A. Pashkin, Y. Nakamura and  J. S. Tsai, 

https://doi.org/10.1038/nature06141Nature 449, 588-590 (2007).

Astafiev2015 Z. H. Peng, Yu-xi Liu, J. T. Peltonen, T. Yamamoto, J. S. Tsai, and O. Astafiev, 

https://link.aps.org/doi/10.1103/PhysRevLett.115.223603Phys. Rev. Lett. 115, 223603 (2015).

Cassidy M. C. Cassidy, A. Bruno, S. Rubbert, M. Irfan, J. Kammhuber, R. N. Schouten, A. R. Akhmerov, L. P. Kouwenhoven, 

https://doi.org/10.1126/science.aah6640 Science 355, 6328 (2017). 

Hekking2010 N. Didier, Ya. M. Blanter, and F. W. J. Hekking, 

https://doi.org/10.1103/PhysRevB.82.214507 Phys. Rev. B 82, 214507 (2010). 

Nori2007 J. Q. You, Yu-xi Liu, C. P. Sun, and F. Nori,

https://link.aps.org/doi/10.1103/PhysRevB.75.104516 Phys. Rev. B 75, 104516 (2007).

Nori2009 S. Ashhab, J.R. Johansson, A.M. Zagoskin, F. Nori, 

https://doi.org/10.1088/1367-2630/11/2/023030 New J. Phys. 11, 023030 (2009).

Lukin2004 L. Childress, A. S. Sørensen, and M. D. Lukin, 

https://doi.org/10.1103/PhysRevA.69.042302 Phys. Rev. A 69, 042302 (2004).

Marthaler2011 P.-Q. Jin, M. Marthaler, J. H. Cole, A. Shnirman, and G. Schön, 

https://doi.org/10.1103/PhysRevB.84.035322 Phys. Rev. B 84, 035322 (2011).

Segal2019 B. K. Agarwalla, M. Kulkarni, and D. Segal, 

https://link.aps.org/doi/10.1103/PhysRevB.100.035412Phys. Rev. B 100, 035412 (2019).

Rastelli2019 G. Rastelli and M. Governale, 

https://doi.org/10.1103/PhysRevB.100.085435 Phys. Rev.  B 100, 085435 (2019).

Kurihara1996 N. Hatakenaka and S. Kurihara, 

https://link.aps.org/doi/10.1103/PhysRevA.54.1729Phys. Rev. A 54, 1729 (1996).

Schwab K. Schwab, E. A. Henriksen, J. M. Worlock, and M. L. Roukes, 

https://doi.org/10.1038/35010065 Nature 404, 974 (2000).	

Meschke M. Meschke, W. Guichard, and J. P. Pekola, 

https://doi.org/10.1038/nature05276 Nature 444, 187 (2006). 

Timofeev A. V. Timofeev, M. Helle, M. Meschke, M. Möttönen, and J. P. Pekola, 

https://doi.org/10.1103/PhysRevLett.102.200801 Phys Rev. Lett. 102, 200801 (2009).  


Partanen M. Partanen, K. Yen Tan, J. Govenius, R. E. Lake, M. K. Mäkelä, T. Tanttu, and M. Möttönen, 

https://doi.org/10.1038/nphys3642 Nat. Phys. 12, 460 (2016).

QHV A. Ronzani, B. Karimi, J. Senior, Y.-C. Chang, J.T. Peltonen, C. Chen and J.P. Pekola, 

https://doi.org/10.1038/s41567-018-0199-4 Nat. Phys. 14, 991 (2018).  

QHR J. Senior, A. Gubaydullin, B. Karimi, J. T. Peltonen, J. Ankerhold, and J. P. Pekola, 

https://doi.org/10.1038/s42005-020-0307-5 Commun Phys 3, 40 (2020).

TPG2019 G. Thomas, J. P. Pekola, and D. S. Golubev, 

https://link.aps.org/doi/10.1103/PhysRevB.100.094508
Phys. Rev. B 100, 094508 (2019).

Scovil1959 H. E. D. Scovil and E. O. Schulz-DuBois, 

https://link.aps.org/doi/10.1103/PhysRevLett.2.262Phys. Rev. Lett. 2, 262 (1959).

Scovil1967 J. E. Geusic, E. O. Schulz-DuBios, and H. E. D. Scovil, 

https://link.aps.org/doi/10.1103/PhysRev.156.343Phys. Rev. 156, 343 (1967).

Giazotto2009 F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola,

https://doi.org/10.1103/RevModPhys.78.217 Rev. Mod. Phys. 78, 217 (2006). 

Savage1992 Y. Mu and C. M. Savage, 

https://link.aps.org/doi/10.1103/PhysRevA.46.5944Phys. Rev. A 46, 5944 (1992).

Scully2011M. O. Scully, K. R. Chapin, K. E. Dorfman, M. B. Kim, and A. Svidzinsky, 

https://doi.org/10.1073/pnas.1110234108Proceedings of the National Academy of Sciences, 108 15097 (2011).

Scully2017 Sheng-Wen Li, M. B. Kim, G. S. Agarwal, and M. O. Scully, 

https://link.aps.org/doi/10.1103/PhysRevA.96.063806Phys. Rev. A 96, 063806 (2017).

Kokkoniemi2019 R. Kokkoniemi, J. Govenius,  V. Vesterinen  R. E. Lake, A. M. Gunyhó, K. Y. Tan, S.
Simbierowicz, L. Grönberg, J. Lehtinen, M. Prunnila, J. Hassel, A. Lamminen, O-P. Saira and M. Möttönen,

https://doi.org/10.1038/s42005-019-0225-6Commun Phys 2, 124 (2019).

Giazotto2020F. Paolucci, N. Ligato, V. Buccheri, G. Germanese, P. Virtanen, and F. Giazotto,

https://arxiv.org/abs/2003.05966arXiv:2003.05966 [cond-mat.mes-hall] (2020).


Orlando T. P. Orlando, J. E. Mooij, Lin Tian, Caspar H. van der Wal, L. S. Levitov, S. Lloyd, and J. J. Mazo, 

https://doi.org/10.1103/PhysRevB.60.15398 Phys. Rev. B 60, 15398 (1999).  

Devoret2017 N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M. Sliwa, M. H. Devoret, 

 https://doi.org/10.1063/1.4984142Appl. Phys. Lett. 110, 222603 (2017)..
 
Esteve2008  A. Palacios-Laloy, F. Nguyen, F. Mallet, P. Bertet, D. Vion and D. Esteve, 

 https://doi.org/10.1007/s10909-008-9774-xJ. Low Temp. Phys.  151 1034 (2008).


SiegmanBook A. E. Siegman  Lasers (University Science Books, Mill Valley, 1986).

Ke2019 Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, and A. N. Poddubny, 

https://link.aps.org/doi/10.1103/PhysRevLett.123.253601Phys. Rev. Lett. 123, 253601 (2019).

Wigner1932 E. Wigner, 

https://link.aps.org/doi/10.1103/PhysRev.40.749Phys. Rev. 40, 749 (1932).

Groenewold1946  H. Groenewold, 

https://doi.org/10.1016/S0031-8914(46)80059-4 Physica,  12, 405 (1946).




 

 




