Thermometry Based on a Superconducting Qubit
Dmitrii S. Lvov Sergei A. Lemziakov Elias Ankerhold Joonas T. Peltonen Jukka P. Pekola


APS/123-QED




 dmitrii.lvov@aalto.fi








 
PICO group, QTF Centre of Excellence, Department of Applied Physics, Aalto University, P.O. Box 15100, FI-00076 Aalto, Finland



§ ABSTRACT 
We report temperature measurements using a transmon qubit by detecting the population of the first three levels of it, after employing a sequence of π-pulses and performing projective dispersive readout. We measure the effective temperature of the qubit and characterize its relaxation and coherence times τ_1,2 for three devices in the temperature range 20 – 300 mK. Signal-to-noise (SNR) ratio of the temperature measurement depends strongly on τ_1, which drops at higher temperatures due to quasiparticle excitations, adversely affecting the measurements and setting an upper bound of the dynamic temperature range of the thermometer. The measurement relies on coherent dynamics of the qubit during the π-pulses. The effective qubit temperature follows closely that of the cryostat in the range 100 - 250 mK. We present a numerical model of the qubit population distribution and compare it favorably with the experimental results.


Thermometry Based on a Superconducting Qubit
    J. P. Pekola
    September 9, 2024
============================================





§ INTRODUCTION


In recent years theoretical proposals for quantum thermometry have attracted a lot of attention due to prospects of monitoring temperature with minimal invasion, properties of quantum systems with high accuracy at extremely low temperatures and opportunities for studying open quantum systems and quantum thermodynamics <cit.>. However, the number of conceptually different experimental realizations of quantum thermometers based on coherent quantum objects remains quite limited. For example, quantum thermometry was demonstrated using a nuclear magnetic resonance (NMR) techniques for polarization measurements of spin ensembles in solutions at room temperature <cit.> and using coherent quantum N00N-states <cit.>. Generation and transfer of the latter was recently realized on the basis of a superconducting qubit <cit.>. Also there are noteworthy experiments on thermometry based on nitrogen-vacancy (NV) colour centres in diamond <cit.>. Another possible platform for quantum thermometers is semiconductor quantum dots, which were used for spectroscopy-type measurements in sub-kelvin temperature range in Ref. <cit.>, exhibiting good signal-to-noise ratio and fast operation but limited coherence.

Prospects of utilizing superconducting qubits as thermometers are of dual research interest: they can facilitate studies of open quantum systems and interaction between the qubits and their environment, which is vital in the context of building large multiqubit systems such as quantum computers and simulators. On the other hand qubit thermometers can be used for realization of quantum thermometry protocols. Quantum electrodynamics (QED) circuits based on coherent superconducting transmon qubits were primarily used to measure the residual qubit population and the corresponding residual temperature in a number of works, for example in <cit.>. The residual population of a transmon in a 3D cavity was measured by Rabi oscillations on the e-f transition in Ref. <cit.> and correlation measurements in Ref. <cit.>. In Ref. <cit.> another technique of applying sequences of π-pulses was used to measure the effective qubit temperature by observing the population distribution for the three lowest energy states of a transmon in a planar geometry. A different approach was demonstrated in Ref. <cit.>, where thermometry was done by reflection-type spectroscopy measurements of a transmon qubit strongly coupled to a waveguide.

Yet qubits are sensitive to various sources of dissipation and decoherence, and the diverse nature of their environment (see, for example, <cit.>). Therefore in practice, identification of a specific thermal bath, leading to thermalization of the qubit, can be quite challenging; i.e., it may be difficult to identify what temperature the qubit is actually measuring.

In this paper we present a comprehensive and concrete study of employing a superconducting qubit for quantum thermometry at sub-kelvin temperatures. The measurement technique used in our work relies on applying sequences of π-pulses for manipulating the qubit population distribution, requiring sufficient coherence of the qubit, and then performing consecutive qubit state readout, which was described in Ref. <cit.>. We also study temperature dependencies of qubit relaxation and coherence times. They determine the working range of the thermometer in terms of temperature, and they give conatraints on the pulse durations in the measurement. While showing a solid thermalization of the qubit to the cryostat stage, the technique allows us to measure the population distribution at temperatures up to 200-220 mK, when the second excited state of the qubit was already significantly populated. This exceeds the applicability range of methods based on consideration of only the lowest qubit levels. At the same time the technique is less dependent on the signal-to-noise ratio (SNR) of the measurement setup. We compare the experimental results with a numerical model of qubit population measurements, being in a good qualitative agreement. Within the model we investigate the influence of the measured parameters on the population and effective qubit temperature obtained as the result. Finally, we address the factors limiting the dynamical range of a superconducting qubit thermometer. The paper is organized as follows: in Section <ref> we provide a brief description of a realistic model of a qubit being in thermal equilibrium with heat baths. Section <ref> explains the principle of the population distribution measurements. Section <ref> provides a brief description of the samples used in the experiment, and discussion on the main experimental results is given in Section <ref>, and finally, limits of the thermometer and ways of mitigation are discussed in Section <ref>. The paper is followed by appendices providing additional experimental data, details of the numerical simulations, experimental setup and sample fabrication.





§ QUBIT POPULATION AND EFFECTIVE TEMPERATURE


Let us consider a qubit in thermal equilibrium with a heat bath having temperature T. In this case the qubit is in a thermal state, which means that its density matrix is diagonal and the probabilities p_g and p_e, respectively, to be in the ground |g⟩ or excited state |e⟩ (qubit population) are Boltzmann distributed as

    p_e/p_g = e^-ħω_ge/k_B T,

where ħω_ge is the separation between these two states, and k_B is the Boltzmann constant. Thus, by measuring the population ratio of the qubit one obtains the temperature of the thermal bath. 

In a more realistic scenario when a qubit is in steady state with respect to several uncorrelated heat baths (with temperatures T^(i), where i enumerates baths), the state of the qubit is still a thermal-like state. But the temperature characterizing this state is not necessarily the temperature T^(i) of any of the thermal baths. Then Eq. (<ref>) could be interpreted as giving an effective qubit temperature T_eff. In this sense, Eq. (<ref>) can be rewritten as a definition of T_eff as

    T_eff = ħω_ge/k_B(- lnp_e/p_g)^-1,

which coincides with the true temperature in equilibrium.

For the case of ohmic heat baths, one can find the energy relaxation rate γ_1 for the qubit towards T_eff using Fermi golden rule approach. The detailed description is given in Appendix <ref>. The qubit energy relaxation rate writes

    γ_1 = ∑_i(Γ_↓^(i)+Γ_↑^(i))=∑_iγ_1^(i)[2n(T^(i), ω_ge)+1]
    
        =∑_iγ_1^(i)ħω_ge/2 k_B T^(i),

where n = [exp(ħω_ge/k_B T) - 1]^-1 is the Bose-Einstein distribution.

The population ratio for the qubit then writes

    p_e/p_g=Σ_i Γ_↑^(i)/Σ_i Γ_↓^(i)=Σ_i γ_1^(i) n(T^(i), ω_ge)/Σ_i γ_1^(i)[n(T^(i), ω_ge)+1]




§ POPULATION MEASUREMENTS


In this work we measured population distribution of several transmon qubits to obtain temperature using the technique described in <cit.>, which is based on applying different sequences of π-pulses and then performing a dispersive readout of the qubit state. Here we do not assume that the transmon is a two-level qubit, but include the higher levels, especially the third one |f⟩ in the analysis and protocols. We presume, that the population of a qubit thermalized at temperature T follows the Maxwell-Boltzmann distribution p_i = exp(-E_i/k_B T) /Z, with the partition function Z = ∑_i=1^3exp(-E_i/k_B T) and the i-th level energy E_i.
A transmon qubit typically has a weak anharmonicity α ( |α| ≪ħω_ge) <cit.>, so that the energy separations ħω_ge and ħω_ef are related by ω_ef = ω_ge + α/ħ≅ω_ge. For example, if a transmon with ω_ge/2π=6.65GHz and α/h = -230MHz is thermalized at temperature 300mK we can estimate the population ratios p_e/p_g ≈ 0.35 and p_f/p_g ≈ 0.12. Thus, it is important to take into account the finite population of the f-level. Limiting to the three lowest energy states of a transmon, the density matrix of a thermalized transmon reads 

    ρ̂=p_g|g⟩⟨g|+p_e|e⟩⟨e|+p_f|f⟩⟨f|.


Let us denote the readout response of a projective measurement of the qubit being in state |i⟩ as φ_i, which we will call “pure state response” further in the text. Then the result φ of the qubit state measurement is a linear combination of the form φ = p_g φ_g+p_e φ_e+p_f φ_f. The π_ge(π_ef)-pulses effectively swap pairs of the population probabilities p_g and p_e (p_e and p_f), which results in getting different readout outcomes, which we denote x_i and y_j, i,j = 0,1,2. An example of such a measurement is shown in Fig. <ref> a, c. We can then write a system of linear equations as

    x_0 = p_g φ_g + p_e φ_e + p_f φ_f,    (no pulse)
       x_1 = p_e φ_g + p_g φ_e + p_f φ_f,      (π_ge)
       x_2 = p_e φ_g + p_f φ_e + p_g φ_f,     (π_geπ_ef)
       y_0 = p_g φ_g + p_f φ_e + p_e φ_f,      (π_ef)
       y_1 = p_f φ_g + p_g φ_e + p_e φ_f,     (π_efπ_ge)
       y_2 = p_f φ_g + p_e φ_e + p_g φ_f.    (π_efπ_geπ_ef)




We can find functional dependencies between p_i, i∈{g,e,f} and x_j, y_k with j,k = 0,1,2 which are valid for arbitrary pure state responses φ_i. These quantities are temperature dependent given by A=p_g-p_e/p_g-p_f, B=p_e-p_f/p_g-p_e and C=p_e-p_f/p_g-p_f, which can also be expressed through the readout results x_j, y_k (see Table <ref>).
To understand the temperature behaviour of functions A, B and C, we can take the low temperature limit ħω_ge≪ k_B T, also assuming p_f → 0. Then we 
can rewrite A ≈ 1 - p_e/p_g, B ≈ C ≈ p_e/p_g or 

    A ≈ 1 -    exp(-ħω_ge/k_B T), 
    
        B ≈ C ≈   exp(-ħω_ge/k_B T).




§ SAMPLE DESCRIPTION


We used transmon qubits on three different chips I, II and III (which is reflected in the name, e.g. R2-I, R4-I, R3-II, Q2-III) made in different fabrication cycles. All transmon qubits had a fixed transition frequency ω_ge/2π≈6.5GHz or flux-tunable spectrum (ω_ge/2π = 5.5-6.8GHz) and charging energy E_c/h=210-230MHz (see Fig. <ref>) and they were weakly coupled (g/2π = 35-45 MHz) to a coplanar waveguide (CPWG) λ/4 readout resonator with the internal quality factor Q_int∼ 10^5, external quality factor Q_ext≈ 8· 10^4 and fundamental frequency ω_r/2π=4.5-6.0GHz for a dispersive state readout (|Δ_r| = |ω_ge-ω_r| ≫ g).
Experimentally found parameters of the devices are shown in Table <ref>.
The fabrication details are given in Appendix <ref>.





§ RESULTS AND DISCUSSION


In the experiment we measured qubit relaxation and coherence times τ_1 and τ_2, and qubit effective temperature T_eff, while controlling the temperature of the mixing chamber (MXC) stage of the dilution refrigerator. After a certain MXC stage temperature T_MXC was set, we waited for 10 min for thermalization of the sample stage and then performed a series of measurements, consisting of five repetitions of the following measurement sequence: (i) measurement of the MXC stage temperature T_MXC, (ii) population distribution measurement x_i, y_j, (iii) measurement of qubit relaxation time τ_1, (iv) Ramsey oscillations measurement to get τ_2^R and ω_ge and (v) Hahn echo measurement with the decay time τ_2^E. Each of these iterations took approximately 5 min. Qubit excitation pulses were applied via an AC drive antenna, and the readout pulses were sent through a waveguide. The measurement setup is shown in Fig. <ref>.

First we address the temperature dependence of the qubit decay time τ_1, shown in Fig. <ref> a. It stays almost constant in temperature range 20-150 mK. Higher than T_MXC≈170mK a noticeable effect of quasiparticle induced relaxation appears and already at T_MXC≈250-270mK decay time drops to below 0.5 μs, which is shorter than the readout pulse duration Δ t_RO=1.5-2.0 μs and becomes comparable to the duration of π-pulse sequences used for the population measurements. To compare τ_1 with theoretical prediction we use the quasiparticle relaxation model in Ref. <cit.> that gives the relaxation rate for a transmon as 

    γ_1^qp = 1/π   ω_p^2 /ω_ge{ x_qp√( 2Δ/ħω_ge). 
        + . 4e^-Δ/k_B Tcosh( ħω_ge/2 k_B T)   K_0 ( ħω_ge/2 k_B T) }


with the plasma frequency ω_p = √(8 E_J E_C)/ħ, equilibrium quasiparticle density x_qp = √(2π k_B T/ Δ) exp(-Δ/ k_B T), modified Bessel function of the second kind K_0 and the superconducting gap Δ. Expression (<ref>) provides a good correspondence with the experimentally observed drop of τ_1 at temperatures above 220 mK for the gap value Δ_Al/e = 180 μV (see Fig. <ref> b).  Thus, the sharp exponential drop of quasiparticle-induced suppression of τ_1 restricts the measurement range for Al-based qubits using long pulse sequences. This is clearly demonstrated in the experimental data discussed below.



As to the dephasing during the π-pulses, typical length Δ t_π of single π-pulses varied between 20 and 220 ns, with the lower limit Δ t_min determined either by the qubit anharmonicity (Δ t_min> h/|α| ≃ 1/200MHz = 5ns) or by the dynamical range of instruments and total line attenuation of the setup (see Fig. <ref>). We present experimental data of the temperature dependence of the dephasing time τ_φ in Fig. <ref> b. The times τ_φ^R and τ_φ^E are extracted from the measurements of Ramsey oscillations and Hahn echo measurements, respectively, as τ_φ^R,E = 1/γ_φ^R,E = 1/(γ_2^R,E - γ_1/2), where the decoherence rates γ_2^R,E describe Ramsey oscillations and Hahn echo experiments. While in the low-temperature limit τ_φ^R,E⪆τ_1, with τ_φ^R exhibiting a clearly visible saturation below 50-70 mK, the dephasing time starts decreasing at lower temperatures than τ_1 and reaches 1-2 μs at 150 mK, becoming difficult to measure at 220-250mK due to the finite ring-up time of the readout resonator ∼ Q_ext/ω_r≃ 1 μs. Similarly to the analysis of the relaxation time, we can estimate influence of quasiparticles on the dephasing, using the approach presented in Ref. <cit.>. The pure dephasing rate
due to tunneling of equilibrium quasiparticles is given by expression

    γ_φ = E_c /πħk_B T/Δexp(-Δ/k_B T) ,

and appears to be negligible in the experimental temperature range. However, we can estimate the dephasing mechanism caused by the Josephson frequency fluctuations due to quasiparticles occupying the Andreev bound states in the tunnel junction of the transmon. The latter writes <cit.>

    γ_φ∼ 4 πω_p^2/ω_ge√(x_qp^A/N_e),

where x_qp^A = exp(-Δ/k_B T) is the equilibrium quasiparticle occupation of the Andreev states and N_e = ζ^-1 g_T/2 g_K is the effective number of channels in the junction with the factor ζ^-1∼ 10^3 - 10^5 describing transparency of the junction in the subgap regime <cit.>, g_T = 1/R_n being the junction conductance and the conductance quantum g_K = e^2/h. In our case R_n ≈ 4.4-5.5kΩ, giving g_T/2 g_K ≈ 2.3-2.9, and the dependence (<ref>) for the realistic parameters is shown as dashed lines in Fig. <ref> b.

Population ratios A, B and C can be extracted from the measurement data x_i, y_j by three different methods each, which we denote as A_1, A_2, A_3, B_1, …, C_3 (see Tab. <ref>). Qualitatively speaking, the difference between the methods is the order of the π-pulses in the control sequence. Because of the qubit relaxation and relatively low population of the f-level, the resulting signals can have different signal-to-noise (SNR) ratio, giving us some guidance on how to choose the best one (see Fig. <ref> and discussion in Sec. <ref> and Appendix <ref>).

An example of population measurements in terms of datasets A_1,A_2,… C_3 (see exact expressions in Appendix in Table <ref>) for device R4-I can be seen in Fig. <ref> c. While according to expressions (<ref>), (<ref>), in the low-temperature limit A→1 and B,C → 0, there is a clear saturation at p_e/p_g ≈ 0.024, which agrees with T_eff≈ 85mK, calculated by Eq. (<ref>). Above 170 mK the population readings have increased deviation from the analytical expressions caused by the enhanced relaxation. 



Measurements of the qubit effective temperature for devices R2-I, R4-I and Q2-III extracted from A_2 are shown in Fig. <ref> d and demonstrate a linear behaviour in the temperature range 120-230 mK close to one-to-one correspondence with the temperature of the mixing chamber stage T_MXC of the cryostat. At lower temperatures of the cryostat T_eff exhibits saturation, reaching T_eff^0≈ 65-85mK, which coincides with the saturation temperature estimated from the population in Fig. <ref> c. This elevated temperature compared to T_MXC is most probably caused by imperfect shielding against thermal photons and insufficient filtering of the measurement lines. For comparison, T_eff extracted using all nine methods A_1, A_2, A_3, B_1, …, C_3 are shown in Fig. <ref> a. A thorough explanation of the deviations between the experimental data and the expected behaviour above 170 mK in Figs. <ref> c, d requires consideration of the time evolution of the qubit population due to its finite lifetime, which we discuss in the following subsection.



 §.§  Comparison of the Experimental Data with the Numerical Model of the Qubit Population Evolution


In order to understand the deviation of the population functions A,B and C from the exact formulas in Tab. <ref> and the qubit effective temperature T_eff from the cryostat temperature T_MXC at higher temperatures, we implement a numerical model which accounts for time evolution of the population distribution of a qubit. First, let us discuss important assumptions underlying the model, which are based on the experimental data. We assume, that the control sequences of π-pulses are short compared to the qubit dephasing and relaxation times (3Δ t_π < τ_φ, τ_1 ). Then we are interested in tracking the time evolution of the population distribution, described by the diagonal terms of the density matrix. The duration of the qubit state readout Δ t_RO is larger than or of the order of the relaxation time (Δ t_RO≳τ_1).
The model considers population of the three lowest states with finite decay times τ_1^ge and τ_1^ef of the ge- and ef-transitions, respectively. 

We begin with a qubit, which is thermalized at a temperature T, setting initial populations p_g,e,f^(0)(T) and then imitate application of π_ge- and π_ef-pulses by instantaneous swapping of the corresponding level populations with finite efficiencies, designated as δ_ge,ef∈ [0,1] (see Fig. <ref> e). For example, a π_ge-pulse with δ_ge=0.9 means that a qubit state with p_g=1.0, p_e=0.0, p_f=0.0 would turn result into p_g=0.1, p_e=0.9, p_f=0.0. Then we calculate the time evolution of the population during the readout time, and multiply it by arbitrary pure state responses φ_i. Next we construct the population ratios A, B and C and calculate the effective temperature (see Fig. <ref> b). A detailed description of the model is provided in Appendix <ref>.

The numerical model can qualitatively reproduce the deviation of the qubit population functions A,B and C from the analytical expressions (see Fig. <ref> in App. <ref>), as well as the resulting deviations of the temperatures T_eff extracted with different methods in the higher temperature range, which is presented in Fig. <ref>. This divergence is a consequence of the temperature-dependent qubit decay times τ_1^ge and τ_1^ef, caused by quasiparticles. For the simulation results presented in Fig. <ref> b, we use frequencies ω_ge and ω_ef and the temperature dependent relaxation time τ_1^ge(T)=2π/(γ_1^qp(T)+γ_1^0), which correspond to the sample R4-I shown as the red line in Fig. <ref> b, and τ_1^ef≈τ_1^ge/2, valid for a transmon qubit <cit.>. For ideal π-pulses, the effective temperatures T_A1, T_A2, …, T_C3 extracted from different ratios A_1,A_2,…,C_3, respectively, demonstrate almost one-to-one behaviour for the temperature range 30-200 mK, when the qubit lifetimes are relatively long. But above 200 mK, when τ_1 is significantly suppressed, there is an increasing deviation between these dependencies. Importantly, this does not indicate lack of thermalization of the qubit, rather than showing the limits of applicability of the given thermometry protocol. Moreover, in this temperature limit the values T_A1, T_A2, …, T_C3 start to depend on the pure state responses φ_i, which eliminates the main advantage of the protocol. In fact, Eqs. (<ref>) are no longer valid, as the qubit population changes significantly during the readout pulse.

Another error in determining the effective temperature of the qubit is the finite efficiency of the pulses. As seen in Fig. <ref> b, the divergence of T_A1, T_A2, …, T_C3 grows monotonously with the temperature and is noticeable already at 100-150 mK, which is a characteristic of the protocol. The plot, showing the relative error between the extracted temperature and the qubit thermalization temperature averaged over different methods A_1,A_2,…,C_3  is shown in Fig. <ref> in Appendix <ref>. 








 §.§  Qubit Thermalization


In order to make sure, that the qubits used in our measurements thermalize to a single thermal bath, which is the substrate, we analyze the relaxation rate γ_1 and photon number n(ω_ge, T) as functions of the qubit effective temperature T_eff. For example, if a qubit is in a steady state with several ohmic baths, then from (<ref>) and (<ref>) we can write down the following temperature dependencies for γ_1 and n(ω_ge, T_eff):

    γ_1(T   )= γ_1^0[ 2 n(ω_ge, T_eff) +1 ]=γ_1^0 ħω_ge/2k_B T_eff ,
       n(ω_ge, T_eff) = γ_1^0/∑_i γ_1^(i) n(ω_ge, T_MXC) + n_0,


where we denote the relaxation rate and the residual qubit population at the base temperature as γ_1^0 and n_0, respectively.
The corresponding experimental data is presented in Fig. <ref>. The photon number dependence of γ_1 (see Fig. <ref> a) exhibits two regimes: the linear growth described by Eq. (<ref>), where the derivative ∂γ_1/∂ n is twice bigger than the offset n_0, and the sharp rise due to the quasiparticle relaxation.

The dependence between the photon number n(ω_ge, T_eff), which corresponds to the qubit effective temperature, and n(ω_ge, T_MXC), defined by the cryostat temperature, characterizes the general process of the qubit thermalization (see Fig. <ref> b). In case of a single thermal bath coupled to the qubit, the dependence would be exactly one-to-one, whereas the presented data shows finite photon population n_0 in the low-temperature limit and the slope slightly less than unity, witnessing coupling to other thermal baths. One of those could be a residual radiation field, as well as a bath of TLSes or nonequilibrium quasiparticles with the effective temperature 65-85 mK, which sets the lower limit for the population of the qubit at the lower temperatures of the cryostat.





§ THERMOMETER LIMITATIONS AND WAYS OF IMPROVEMENT




 §.§ Temperature Range


The Al-based transmon thermometers studied in this work showed a relatively narrow working temperature range 60-220 mK. In this subsection we discuss possibilities of widening this regime.

A higher residual temperature of the qubits T_eff = 65-85 mK compared to the base cryostat temperature T_base≈ 20 mK is most likely caused by both improper shielding of the sample stage and filtering of the measurement lines <cit.>, likely due to nonequlibrium quasiparticles <cit.> or TLSes <cit.>. Recent works demonstrated the possibility of achieving effective qubit temperatures of 30-45 mK <cit.>.

When considering the higher temperature limit of the presented method the following aspects should be taken into account: i) consideration of only the three lowest states is valid until the population of the next d-level is negligible, otherwise, at higher temperatures, the system (<ref>) should be supplemented with appropriate terms having φ_d and p_d; ii) suppression of the qubit relaxation time τ_1 due to growing quasiparticle density above 200 mK; iii) sufficiently long coherence time for the qubit state preparation.
In agreement with the model, quasiparticle suppression of τ_1 leads to a significant influence on the population measurements, which in turn affects adversely on determining the effective temperature. This is a hard limit due to the exponential growth of the quasiparticle density with temperature. 

Another possibility is employment of materials with a larger superconducting gap to shift the quasiparticle suppression limit to higher temperatures. This might be done by using Nb-based tunnel junctions <cit.> with the critical temperature T_c^Nb = 7.6-8.7 K or by replacing the tunnel junction with a nanowire made of granular Al having critical temperatures T_c^grAl up to 3.15 K <cit.>, which is 2.5 times larger than critical temperature of bulk aluminium, T_c^Al, bulk = 1.2 K.

The qubit state preparation is also of big importance in the population distribution measurement presented here. In this sense, the qubit coherence time, or more strictly speaking, the dephasing time should be longer than the control pulse sequences (τ_φ≳ 3 Δ t_π). Since the dephasing rate is growing with the temperature increase, it limits duration of the pulse sequences at higher temperatures. 



 §.§ Signal-to-Noise Ratio and Accuracy


The parameters of the experimental setup can be fine-tuned to improve both the temperature measurement accuracy and the range of the cryostat temperatures, so that it reaches the quasiparticle limit. First of all, improvement of the π-pulse efficiency affects positively on the measurement accuracy. Secondly, as mentioned above, the short lifetime of the qubit is the main restriction of the method. We can divide the qubit state decay into two parts: i) decay during the control sequence and ii) decay during the readout. The first process leads to change of the qubit population distribution p_i and can be associated to imperfect control pulses, leading to increased deviation between T_A1, T_A2,…,T_C3. Decay during the readout pulse is defined mainly by the ring-up time of the resonator, and diminishes distinguishability between the qubit states (i.e. difference between φ_i on the I-Q plane), thus reducing the SNR, and also distorting the extracted effective temperature. The main method of improving the SNR is optimizing the ratio χ/ κ for the readout resonator, where χ is the dispersive shift and κ is the linewidth, as well as an appropriate readout frequency detuning for a given χ/ κ, which directly affects the pure state responses φ_i <cit.>.

Another factor directly affecting the efficiency of the π-pulses is temperature dependent properties (specifically S21) of the measurement lines of the cryostat and CPWGs on the chip. In particular, it is known, that superconducting CPWG resonators have a temperature dependent Q-factors and fundamental frequencies <cit.>, which can lead to changes in amplitude of the qubit drive signal. Moreover, the qubit frequency even for the single-junction transmons can slightly vary with temperature (less than 50 kHz in our case). All these factors can lead to an increased deviation of the π-pulse parameters if a calibration procedure is not applied at every cryostat temperature point. 



The minimal temperature measurement error is defined by the Cramer-Rao bound from the quantum Fisher information (QFI) <cit.>, where we can use an expression for the relative error |Δ T/⟨ T⟩|_sm of a single measurement for a two-level thermometer with  energy separation ħω_ge and N-fold degenerate excited level as 

    (Δ T/⟨ T⟩)^2_sm = (N - 1 + e^x_ge)^2/(N - 1) x_ge^2 e^x_ge,

x_ge=ħω_ge/k_B T.
An analogous expression can be derived for the three-level case considered in our work,

    (Δ T/⟨ T⟩)^2_sm = 
        (e^x_ge+x_gf + e^x_ge + e^x_gf)^2/[x_ge^2 e^x_gf + x_gf^2 e^x_ge + (x_ge + x_gf)^2 ] e^x_ge+x_gf,

here x_gf=ħω_gf/k_B T and ω_gf = 2ω_ge - |α|.
Taking into account averaging over N = 2^17 repetitions, one gets the lower bound for the relative error (Δ T/⟨ T⟩)^2 = (Δ T/⟨ T⟩)^2_sm/N ≈ 5.0· 10^-5 or |Δ T/⟨ T⟩| ≈ 0.7 % for ⟨ T⟩ = 65 mK and qubit frequency ω_ge = 7.04 GHz. The noise equivalent temperature NET=√(|Δ T|^2 Δ t_meas) with the measurement time Δ t_meas = 29 s for our setup in the QFI-limited case would be NET_QFI≈ 2.5 mK/√(Hz). These estimates are an order of magnitude lower than the experimental errors. The data shown in Fig. <ref> a has the relative errors |Δ T_A/⟨ T_A ⟩| = 14-42 %, |Δ T_B/⟨ T_B ⟩| = 8.6-9.4 % and |Δ T_C/⟨ T_C ⟩| = 8.7-9.4 % at the base cryostat temperature T_MXC = 22 mK, which corresponds to NET≈ 28 mK/√(Hz).




Notably, if the thermometer has a constant level separation, the relative error grows exponentially at T ≪ħω_ge/k_B 
(see Fig. <ref>), 
which is, however, not observed experimentally in our case, as the qubit effective temperature saturates at T_eff^0≈ 65-85 mK.
Our system is a three-level quantum thermometer with a constant anharmonic spectrum, which refers to sub-optimal quantum thermometers <cit.>. Knowing the dependence of QFI on the spectrum of the system, it is possible to optimize it for a better performance in the experimentally realistic temperature range by choosing an appropriate level separation and anharmonicity.



 §.§ Operation Speed

First of all, the presented method refers to the equilibrium type of quantum thermometry, meaning that the population measurement is performed after the qubit thermalization to the bath, which happens on time scales not shorter than the relaxation time τ_1. In the presented experiment the time delay between the consecutive measurements, which is necessary for the qubit relaxation after the readout and thermalization of the qubit, was ∼ 10 τ_1, which takes up most of the measurement time, and defines the speed of the response of the system. Pulse application and the qubit state readout can be limited to sub-microsecond range. Finally, the total measurement time is defined by the averaging number t_total∼ N_av· 10 τ_1 = 10^6 τ_1, with N_av≃ 10^5, meaning a time scale of seconds or minutes per point, depending on τ_1.



§ CONCLUSION


In this work we experimentally realized equilibrium quantum thermometry of the cryostat temperature, based on the population distribution measurements of the three lowest energy levels of a transmon qubit. The presented technique allowed for temperature measurements in the range of 60-220mK, in which the qubit population follows the Maxwell-Boltzmann distribution. Analysis of temperature dependencies of the qubit relaxation rate and the photon number showed that the qubit was thermalized to the cryostat. The working temperature range was defined by saturation of the effective qubit temperature in the low temperature limit and the quasiparticle suppression of the qubit relaxation time above 200 mK. We compared the experimental data to a numerical model, describing time evolution of the qubit population, showing good agreement. We believe, a similar approach could be exploited for quantum thermometry of various on-chip structures and for studies of coupling between quantum objects to different thermal baths. Finally, this experimental platform allows for exploitation of quantum thermometry algorithms including nonequilibrium techniques for faster measurements and dynamical control of the thermometers for enhanced precision.


The authors are grateful to Gershon Kurizki for fruitful discussions and Yu-Cheng Chang for assistance with measurements and fabrication.
We acknowledge the provision of facilities and technical support by Aalto University at OtaNano - Micronova Nanofabrication Centre and OtaNano - Low Temperature Laboratory. 







§ RELAXATION AND EXCITATION RATES FOR COUPLING TO A THERMAL BATH


Let us consider a quantum oscillator based on a resonator with frequency ω_0 = 1/√(LC) and characteristic impedance Z_0=√(L/C)=50 Ω, where L and C are the inductance and capacitance of the resonator, respectively, to which we connect an ohmic bath with resistance R and temperature T. Then the excitation and relaxation rates Γ_↑, ↓ caused by the bath can be written as <cit.>

    Γ_↑,↓ = Z_0/R∓ω_0/1-exp(±ħω_0β),


where β = 1/k_B T is the inverse temperature and k_B is the Boltzmann constant. For a more general case the resistance R should be replaced by the real part of the bath impedance Z. Thus, we can introduce the quality factor Q=Re  Z/Z_0 and the resonator linewidth γ = ω_0/Q, giving us

    Γ_↑,↓ = γ∓ 1/1-exp(±ħω_0β).

This can be further shortened for the sake of notation brevity to

    Γ_↓ = γ[ n(ω, T)+1 ],
       Γ_↑ = γ  n (ω, T),

where n=1/[1-exp(-ħω_0β)] is the Bose-Einstein distribution.
The total decay rate in this case would be 

    Γ_1 = Γ_↓ + Γ_↑ = γ[ 2n(ω, T)+1 ].


Analogous expressions can be readily applied to the qubit-bath interaction (see for example <cit.>).



 §.§ Effective Temperature of a Two-Level System

Let us consider a qubit with the ground state |g⟩ and excited state |e⟩, thermalized with bath with temperature T. Then the population ratio of the qubit levels should follow the detailed balance principle p_e/p_g=e^-ħωβ = Γ_↑/Γ_↓, and <cit.>

    p_g = Γ_↓/Γ_↓+Γ_↑,  
        p_e = Γ_↑/Γ_↓+Γ_↑.

Also, we can write it in terms of the qubit polarization:

    ⟨σ_z ⟩ = -p_g + p_e = -tanh(ħωβ/2).


Energy relaxation rate is just a sum of the two rates:

    γ_1 = 2 π/τ_1 = Γ_↓+Γ_↑.





 §.§ Transmon Qubit Coupled to Several Baths


Here we will extend the results from the previous section to the case when a quantum object interacts with a number of uncorrelated baths with different temperatures T^(i) and coupling rates Γ_↑,↓^(i). As in previous section, the only requirement for the baths is that the corresponding rates obey the detailed balance principle. In this case we can just sum up the rates when calculating the level population:

    p_g = ∑_iΓ_↓^(i)/∑_i(Γ_↓^(i)+Γ_↑^(i)),
       p_e = ∑_iΓ_↑^(i)/ ∑_i(Γ_↓^(i)+Γ_↑^(i)),
       p_e/p_g=Σ_iΓ_↑^(i)/Σ_iΓ_↓^(i)=Σ_iγ_1^(i) n(T^(i), ω)/Σ_iγ_1^(i)[n(T^(i), ω)+1]


The qubit energy relaxation rate writes

    γ_1 = ∑_i(Γ_↓^(i)+Γ_↑^(i))=∑_iγ_1^(i)[2 n(ω, T^(i))+1].


One can use the same considerations for a resonator (harmonic oscillator). Then again the detailed balance principle should be applied to get p_n+1/p_n = ∑_iΓ_↑^(i)/ ∑_iΓ_↓^(i), which is positive, less then unity and does not depend on n. Thus we can find such a number β_eff = 1/k_B T_eff, to have p_n+1/p_n = exp(-ħωβ_eff), leading to

    n = ∑_iΓ_↑^(i)/∑ _i(Γ_↓^(i)-Γ_↑^(i)) = ∑_iγ_i n(T^(i), ω)/∑_iγ_i= n(T_eff, ω).


Effective resonator temperature is then dependent on the number of photons n:

    T_𝑒𝑓𝑓^r = ħω_r/k_B(lnn+1/n)^-1.


All considerations above can be generalized to the case of an arbitrary number of baths.









§ NUMERICAL MODEL

To describe the dynamics of the qubit population distribution, firstly we neglect effects of the readout resonator, and secondly, restrict ourselves by considering only the three lowest energy states of the transmon, |g⟩, |e⟩, |f⟩.
We assume that the initial qubit state is thermal, which means that its density matrix is diagonal. Then we track the dynamics of only these diagonal terms, characterizing the population distribution, which significantly simplifies the calculations.
We include only sequential decay and excitation processes with rates Γ^ge_↓,↑ and Γ^ef_↓,↑ for g-e and e-f transitions, respectively. Then the free time evolution of the qubit populations can be described by a system of differential equations:

    ṗ_f    = p_eΓ^ef_↑-p_fΓ^ef_↓, 
    ṗ_e    = p_gΓ^ge_↑ + p_fΓ^ef_↓-p_e(Γ^ge_↓ + Γ^ef_↑), 
    ṗ_g    = p_eΓ^ge_↓ -p_gΓ^ge_↑ ,

where dotted functions denote time derivatives. The system has an analytical solution:

    p_f(t)    = ζ_f e^α_0 t + η_f e^α_1 t + ξ_f, 
    
        p_e(t)    = ζ_e e^α_0 t + η_e e^α_1 t + ξ_e, 
    
        p_g(t)    = ζ_g e^α_0 t + η_g e^α_1 t + ξ_g.


We find the following coefficients:

    α_0,1 = 1/2   {-(Γ^ef_↑ + Γ^ef_↓ + Γ^ge_↑ + Γ^ge_↓). 
       ±[(Γ^ef_↑ + Γ^ef_↓ + Γ^ge_↑ + Γ^ge_↓)^2 . 
       - .. 4(Γ^ge_↓Γ^ef_↓ + Γ^ge_↑Γ^ef_↑ + Γ^ge_↑Γ^ef_↓)]^1/2},


    ζ_f    = Γ^ef_↑(p_e^(0) - ξ_e) - (Γ^ef_↓ + α_1)(p_f^(0)-ξ_f)/α_0-α_1, 
    ζ_e    = α_0+Γ^ef_↓/Γ^ef_↑ζ_f, ζ_g = Γ^ge_↓/Γ^ef_↑α_0+Γ^ef_↓/α_0+Γ^ge_↑ζ_f, 
    η_f    = (Γ^ef_↓ + α_0)(p_f^(0)-ξ_f) - Γ^ef_↑(p_e^(0) - ξ_e)/α_0-α_1, 
    η_e    = α_1+Γ^ef_↓/Γ^ef_↑η_f, η_g = Γ^ge_↓/Γ^ef_↑α_1+Γ^ef_↓/α_1+Γ^ge_↑η_f,


    ξ_f = 1/𝒵Γ^ef_↑Γ^ge_↑/Γ^ef_↓Γ^ge_↓, ξ_e = 1/𝒵Γ^ge_↑/Γ^ge_↓, ξ_g = 1/𝒵,

where the canonical partition function 𝒵=1+Γ^ge_↑/Γ^ge_↓ + Γ^ef_↑Γ^ge_↑/Γ^ef_↓Γ^ge_↓ imposes the normalization condition p_g + p_e + p_f = 1.

The schematic representation of the model is presented in Fig. <ref> c. We prepare the qubit with initial population p_g,e,f^(0) and, for a fully thermalized qubit in the steady state (t→∞),  observe residual populations p_i^(∞) = ξ_i. To mimic state manipulation with microwave pulses, we define 
pulse matrices acting on a population vector p⃗(t) = (p_g(t), p_e(t), p_f(t))^T:

    M_ge   = [ 1-δ_ge   δ_ge      0;   δ_ge 1-δ_ge      0;      0      0      1 ],
    
        M_ef   = [      1      0      0;      0 1-δ_ef   δ_ef;      0   δ_ef 1-δ_ef ],

where we introduce a phenomenological efficiency parameter 0≤δ≤ 1 to simulate effects such as variation of pulse length and amplitude, detuning and finite fidelity. 
Unlike the real control pulses, these matrices act on the qubit population instantaneously. It is a realistic assumption if the pulse lengths are much smaller than the relaxation times. In this case the π-pulse duration corresponds to the time delay between application of corresponding instantaneous population transformations.



We studied influence of the pulse efficiency on the error of the qubit temperature measurement (see Fig. <ref>). Our model predicts the relative error less than 10% at 200 mK, which is comparable with the experimental data.





§ ERROR ESTIMATION



 §.§ Mutual Dependencies

As described in the main text, for thermometry we use temperature dependent functions A, B and C and each of them could be found in three different ways (see Table <ref>), which gives nine different values for temperature. Though the calculation of these parameters based on six independent x_j, y_k, only six of these functions A_i, B_i, C_i are mutually independent. That could be easily shown from relation A_i× B_i=C_i, i ∈{1,2,3}. That also leads to the mutual dependencies for the corresponding measurement errors. Indeed, one can get the same relation between the relative errors of A, B and C:

    (Δ C/C)^2 = (Δ A/A)^2 + (Δ B/B)^2.

From this equation one can get the relation between the absolute errors:

    (Δ C)^2 = (C/A)^2 (Δ A)^2 + (C/B)^2 (Δ B)^2.


In the low temperature limit we assume p_f = 0, p_e≪ p_g ≈ 1. Using Eqs. (<ref>), (<ref>) we can get C/B ≈ 1 and C/A ≈ p_e, and as a result (Δ C)^2 ≈ (Δ B)^2 + p_e^2 (Δ A)^2. If absolute errors for A, B and C are of the same order of magnitude (it is valid for most of the experimental results), absolute and relative errors for B and C are approximately equal: (Δ C)^2 ≈ (Δ B)^2 and |Δ C/C| ≈ |Δ B/B|.

In order to find the relation between the errors of A, B and C and the measured temperature T one can use logarithmic derivatives of the functions A(T), B(T) and C(T). In order to simplify the final expressions, in the low temperature limit expressions we can write

    |Δ T/T|_A = e^x/x|Δ A| = e^x-1/x|Δ A/A|, 
       |Δ T/T|_B = (e^x-1)^2/x e^x|Δ B| = e^x-1/x e^x|Δ B/B|, 
       |Δ T/T|_C = e^x/x|Δ C| = 1/x|Δ C/C|,

where x = ħω_ge / k_B T. At low temperatures x →∞, Eqs. (<ref>) and (<ref>) give approximately equal errors, meaning that  the relative errors of the temperatures T_B and T_C are similar in the relevant temperature range. That agrees very well with the data presented in this paper (see Sec. <ref> and <ref>) and in Ref. <cit.>. 

Note, that although the coefficient between the relative error of T and the relative error of A in Eq. (<ref>) is e^ħω_ge / k_B T times higher than the same coefficients in Eqs. (<ref>) and (<ref>), the A itself is e^ħω_ge / k_B T times higher than B and C, so it means that relative error for T_A could be the same as for T_B and T_C. The exact relation between errors for T_A and for T_B and T_C depends on particular pure state responses and voltage noise, as will be shown in the next subsection.



 §.§ Further Error Analysis


One can notice that columns 1, 2 and 3 in Table <ref> have the same structure: A, B and C from one column could be written in a form A = b/c, B = a/b and C = a/c, where a, b and c are the measured voltage differences. For example, in the first column a = x_2 - y_2, b = x_0 - x_1 and c = y_0 - y_1. Also, one can write these voltage differences in the following form: a = (p_e - p_f)Δφ_i, b = (p_g - p_e)Δφ_i, c = (p_g - p_f)Δφ_i. Here Δφ_i = φ_j - φ_k for i,j,k ∈{g,e,f}, i≠ j ≠ k. 

One can then make an additional assumption that the variance of a single measurement x_i, y_i is a sum of the voltage noise variance (defined by the measurement setup) and the variance of protective measurements with operator M̂ of the qubit in a thermal state with density matrix ρ̂_th: σ_x,y^2 = Tr(ρ̂_thM̂^2) - Tr(ρ̂_thM̂)^2. Based on this, one can find absolute errors for measurement differences a, b, c:

    (Δ a)^2     = 2 p_e p_f Δφ_i^2
        + (p_e+p_f)p_g(Δφ_j^2 + Δφ_k^2) + σ_a^2, 
    
        (Δ b)^2     = 2 p_e p_g Δφ_i^2
        + (p_e+p_g)p_f(Δφ_j^2 + Δφ_k^2) + σ_b^2, 
    
        (Δ c)^2     = (p_g+p_e) p_f Δφ_i^2
       + (p_g+p_f)p_e(Δφ_j^2 + Δφ_k^2) + σ_b^2,

here σ_a^2, σ_b^2, σ_c^2 are the measurement voltage variances, which in general are not the same for different parameters a, b, c and could be found from the measurements. Relative errors for the voltage differences a, b and c could be written as:

    (Δ a/a)^2 = 2 p_e p_f/(p_e-p_f)^2    + (p_e+p_f)p_g/(p_e-p_f)^2 F(Δφ_i)
    
        +    σ_a^2/Δφ_i^2 (p_e-p_f)^2, 
    (Δ b/b)^2 = 2 p_e p_g/(p_g-p_e)^2    + (p_g+p_e)p_f/(p_g-p_e)^2 F(Δφ_i)
    
        +    σ_b^2/Δφ_i^2 (p_g-p_e)^2, 
    (Δ c/c)^2 = (p_e+p_g) p_f/(p_g-p_f)^2    + (p_g+p_f)p_e/(p_g-p_f)^2 F(Δφ_i)
    
        +    σ_c^2/Δφ_i^2 (p_g-p_f)^2,

where F(Δφ_i) = (Δφ_j^2 + Δφ_k^2)/Δφ_i^2, i≠ j ≠ k. Note, that function F(Δφ_i) reaches the minimum value of 0.5 when Δφ_i = 2Δφ_j = 2Δφ_k, i≠ j≠ k. This gives a condition for the optimal difference of the pure state responses. Importantly, this condition can be satisfied only for the one out of three differences Δφ_i at the same time. 

In order to find relative errors for the A, B, C, one can use the following expressions:

    (Δ A/A)^2 = (Δ b/b)^2 + (Δ c/c)^2, 
    (Δ B/B)^2 = (Δ a/a)^2 + (Δ b/b)^2, 
    (Δ C/C)^2 = (Δ a/a)^2 + (Δ c/c)^2.

Using Eqs. (<ref>, <ref>, <ref>) for relative errors for a, b and c, one can substitute it to the equations above and get the final (but extremely bulky) result.

As in a previous subsection, one can find absolute errors for A, B, C in the low-temperature limit (p_f = 0, p_e≪ p_g ≈ 1):

    (Δ A)^2 ≈ 2 p_e + F(Δφ_i) + σ_b^2/Δφ_i^2 + σ_c^2/Δφ_i^2, 
       (Δ B)^2 ≈ p_e F(Δφ_i) + σ_a^2/Δφ_i^2, 
       (Δ C)^2 ≈ p_e F(Δφ_i) + σ_a^2/Δφ_i^2.

These equations illustrate several important features of errors of the functions A, B and C. Firstly, the relative errors |Δ A/⟨ A ⟩|, |Δ B/⟨ B ⟩| and |Δ C/⟨ C ⟩| grow as 1/√(p_e)∝exp(ħω_ge/2k_B T) when T → 0. Secondly, precision of the method could be limited by voltage noise, as it is for the measurements presented in this paper. Thirdly, in general, measurement precision highly depends on the differences Δφ_i of the pure state responses. On one hand, it gives some room for the measurement protocol optimization in terms of maximization of the Δφ_i^2 and minimization of F(Δφ_i). But, unfortunately, this optimization should be done for each qubit separately. Finally, one can obtain the minimum possible values for the absolute errors in the low temperature limit: (Δ A)^2 ≈ 2 p_e +0.5,  (Δ B)^2 ≈ 0.5 p_e, (Δ C)^2 ≈ 0.5 p_e. 

To sum up, the precision of errors for the temperature estimations from A, B and C measurements are limited by voltage noise and strongly depend on the differences Δφ_i of the pure state responses. So, in general, out of nine possible ways to extract temperature, there is no unique or universal answer as to which one is the optimal one. Instead, for every particular sample and measurement setup one needs to compare them and choose the best one.



§ SAMPLE FABRICATION


The fabrication was done by a multi-step electron-beam lithography (EBL) in a 100 kV EBL system Raith EBPG 5200.
First, pristine high-resistive (ρ > 10^4 Ohm·cm) undoped Si ⟨100⟩ 200 mm wafers were RCA-cleaned and magnetron sputtered with 10 nm of Al and 100 nm of Nb without a vacuum break. Then the wafers were spincoated with 400 nm thick layer of AR-P 6200.13 e-beam resist by spinning at 4000 rpm for 60 s and then baked at 160^∘C for 9 min. After that the samples were cleaved into ∼35×35 mm^2 pieces and loaded into the EBL system. Each wafer piece was used for fabrication of nine 7× 7 mm^2 chips. The exposure of the ground plane was done with the dose 350 μC/cm^2, 200 nA beam current and a 50 nm step size. 
After that the e-beam resist was developed in AR-P 546-600 for 3 min 30 s and rinsed with IPA for 3 min. The etching of Nb was done in a RIE tool Oxford Instruments Plasmalab 80 Plus. Prior to the etching the empty chamber with a dummy quartz wafer was precleaned with a CF_4 +O_2 process using 100 sccm and 15 sccm flows of CF_4 and O_2 respectively at the pressure 600 mTorr, bias voltage 200 V and RF power 200 W during 10 min. Then the samples were loaded into the chamber and pumped to 1·10^-5 mbar. The RIE etching was done with a SF_6 +Ar process at 20 sccm flow of SF_6 and 10 sccm flow of Ar at the pressure of 15 mTorr, V_bias=360 V and RF power 100 W during 90 s. After that residues of the resist were removed with an O_2 descum process at the 40 sccm flow, pressure 250 mTorr, V_bias=320 V and RF power 150 W.
Then the sample was put into acetone for 4 min, rinsed with IPA and dried with N_2 gun. Removal of the Al stopping layer was done using AZ351B developer during 2 min. Then the sample was thoroughly rinsed with deionised water (DIW) for 10 min and dried by putting in IPA for 5 min with consecutive N_2 drying.



The next step was forming of Josephson junctions (JJs) (see Fig. <ref> c). For that the sample was spincoated by a MMA/PMMA e-beam resist bilayer. A 800 nm thick layer of MMA EL-11 was applied by two steps of spincoating at 4000 rpm and baking at 160^∘C for 2 min. Then a 200 nm layer of PMMA A4 was applied by spinning at 4000 rpm and baking at 160^∘C for 10 min. 
Finally, an EBL step was performed to create a suspended mask for further JJ shadow deposition. The main and undercut doses for exposure of submicron structures was 950 μC/cm^2 and 150 μC/cm^2, respectively, with 1 nA beam current and 4 nm stepsize.
The development of the mask was done in a two-step process: (i)
development in MIBK:IPA 1:3 mixture for 25 s and (ii) in methylglycol:methanol 2:1 mixture for 15 s with a following rinsing in IPA for 30 s.
After that the samples were loaded into an e-beam evaporator Plassys MEB700S2-III UHV for JJ deposition. First, the sample was pumped down to 5·10^-7 mbar in the loadlock and an Ar milling was done to remove the resist residues at Ar pressure of 4·10^-4 mbar with ion-beam gun (IBG) parameters of 60 mA beam current and 200 V beam voltage at angles ± 12^∘ for 30 s at each angle. Then the substrate was moved into the oxidation chamber and pumped to 1·10^-7 mbar. After that the chambers were connected and Ti gettering was performed at a rate of 0.5 Å/s for 5 min. Then the first layer of 25 nm thick Al was deposited at angle -12^∘ at 1 Å/s rate at pressure below 5·10^-8 mbar in the oxidation chamber, after which the evaporation chamber was closed and an oxidation was done at a pressure of 18 mbar during 7 min. The second layer of 45 nm thick Al was deposited at +12^∘ at 1 Å/s rate at a pressure below 1·10^-7 mbar in the oxidation chamber. Finally, another oxidation was performed at pressure 20 mbar during 10 min.
The lift-off procedure was done in hot acetone at 52^∘C
for 3 h, then the sample was rinsed with IPA and dried with a N_2 gun.

Now the third step of EBL was done to create patches (see Fig. <ref> c) for a galvanic connection between JJs and the Nb ground plane. The MMA/PMMA e-beam resist bilayer was applied as described above and after that the sample was loaded into the EBL tool. The patch exposure was done with the dose of 950 μC/cm^2 with 30 nA beam current and 20 nm stepsize. 
The resist was developed by putting into MIBK:IPA 1:3 mixture for 30 s and then into methylglycol:methanol 2:1 for 30 s followed by rinsing in IPA for 30 s.
For the patch deposition we again used the Plassys tool. Analogously to the previous step, the sample in the loadlock was pumped to 5·10^-7 mbar and a stronger Ar milling at zero angle was done for 2 min to remove the surface Nb and Al oxides. The Ar pressure was 4·10^-4 mbar and the IBG parameters were the following: 120 mA beam current, 400 V beam voltage.
Then the sample was moved to the oxidation chamber and pumped to 1·10^-7 mbar, after which a deposition of 150 nm of Al was performed at the rate of 2 Å/s and pressure below 1·10^-7 mbar in the oxidation chamber with consecutive oxidation at 20 mbar during 10 min.
Another lift-off step was done similarly to the one described above.

For the low-temperature measurements the wafers were scribed with a diamond tip and cleaved into separate chips. The chips were mounted on a low-temperature sample holder using BF-6 glue and then wire-bonded using FS Bondtec 5330 with Al wire.

*


