Singlet-Doublet Transitions of a Quantum Dot Josephson Junction Detected in a Transmon Circuit

Singlet-Doublet Transitions of a Quantum Dot Josephson Junction Detected in a Transmon Circuit


INTRODUCTION
Superconducting pairing and charging energy are two fundamental interactions that determine the behavior of mesoscopic devices.Notably, when a quantum dot (QD) is coupled to a superconductor, they compete to determine its ground state.A large charging energy favors single-electron doublet occupancy of the dot and thus a spin-1/2 ground state, while a strong coupling to the superconducting leads favors double occupancy in a singlet configuration with zero spin.A quantum phase transition between the singlet and doublet ground state can occur as system parameters such as the dot energy level and the coupling strength are varied.The latter also controls the nature of the singlet ground state, which can be either of the Bardeen-Cooper-Schrieffer (BCS) type or of the Kondo type.The rich phase diagram of the system, as well as its transport properties, are theoretically well captured by an Anderson model with superconducting leads [1][2][3][4][5][6][7][8][9][10].
Recent experiments [32][33][34] on Andreev pair and spin qubits [35][36][37][38] have renewed the interest in quantum dot junctions due to the possibility of tuning the ground state of the system to be in addressable spin states.Knowledge of the phase diagram of the quantum dot junction is also beneficial for realizing proposals for a quantum-dot-based readout of topological qubits [39][40][41].
These developments have highlighted the need for a better fundamental understanding of the quantum dot junction and its dynamics, requiring tools which are not limited by the long integration times of lowfrequency measurements nor by the invasiveness of transport probes.To address this need, we have embedded a fully controllable quantum dot in a microwave superconducting circuit.This experimental choice is motivated by the success of circuit quantum electrodynamics (QED) techniques in the investigation of mesoscopic effects in Josephson junctions [32][33][34][42][43][44][45][46][47][48][49], which stems from its enhanced energy and time resolution compared to lowfrequency transport techniques.In this context, the microwave response of a quantum dot junction has attracted recent theoretical [50,51] and experimental [52] attention.
The core of our experiment is a transmon circuit formed by an island with charging energy E c , coupled to ground via a superconducting quantum interference device (SQUID) formed by a junction with a known Josephson energy E J and a quantum dot junction (Fig. 1(a-c)).The energy-phase relation of the quantum dot junction depends on whether it is in a singlet or doublet state, with a characteristic π-phase shift between the two relations (Fig. 1(d)) [53].The Josephson energies of the two junctions either add or subtract, depending on whether the quantum dot junction is in the singlet or doublet state, and on the value of the applied flux bias (Fig. 1(f)).The two branches of the spectrum give rise to two distinct transition frequencies of the transmon circuit, which can be detected and distinguished via standard circuit QED techniques [54].Therefore, a transition from a singlet to a doublet state will appear as a discontinuous jump in a measurement of the transmon frequency spectrum.Using this method, we have detected the singletdoublet transition and reconstructed the phase diagram of a quantum dot junction as a function of all experimentally controlled parameters in a single device: the energy level of the dot, the tunnel couplings to the superconducting leads, the superconducting phase difference across the quantum dot junction, and also an external Zeeman field.The measured phase boundaries are in agreement with the single-impurity Anderson model with superconducting leads as calculated via the numerical renormalization group (NRG) [2,[55][56][57] methods, and include parameter regimes that have experimentally not been explored before.Finally, we have investigated the rates at which the quantum dot switches between doublet and singlet occupation via real-time monitoring of the transmon circuit, allowing us to determine the switching time-scales of the quantum dot junction parity across the phase transition.

II. DEVICE OVERVIEW
The quantum dot junction under investigation is formed in an epitaxial superconductor-semiconductor InAs/Al nanowire [58].It is defined in an uncovered section of the nanowire where the Al has been etched away, and controlled by three electrostatic bottom gates (Fig. 2(d)).As shown in the circuit of Fig. 2(a), this quantum dot junction is placed in parallel to a second Josephson junction, hereafter referred to as the "reference junction", to form a SQUID.The reference junction consists of a second uncovered segment of InAs on the same nanowire as the quantum dot junction.Its Josephson energy E J can be tuned with a single electrostatic gate via the field effect.
The SQUID connects a superconducting island to ground, resulting in a transmon circuit [59] governed by the Hamiltonian where E c = e 2 /2C Σ , C Σ is the total capacitance of the island to ground.The Josephson potential V (φ) is determined by the phase-dependent energies of the reference junction, V J (δ) = E J (1 − cos δ), and of the quantum dot junction, V s,d (φ): Here, the phase drops across the quantum dot junction (φ) and across the reference junction (δ) are connected according to φ − δ = φ ext , where φ ext = (2e/ )Φ ext is the phase difference resulting from the externally applied magnetic flux through the SQUID loop, Φ ext .The presence of the reference junction serves several purposes.First, it allows us to tune the phase difference at the quantum dot junction by changing Φ ext with the B y component of the magnetic field [60].We operate the device in a regime where the reference junction has a Josephson energy that is substantially larger than that of the quantum dot; for most of the parameter regime explored in our measurements it is larger by more than an order of magnitude [60].This ensures that the phase drop across the reference junction is close to zero, while the phase difference across the quantum dot junction is close to φ ext [61].Second, the ability to tune E J independently of the quantum dot junction configuration ensures that the transition frequencies of the transmon circuit remain inside the measurement bandwidth for all parameter regimes of the quantum dot junction.Finally, the Josephson energy of the reference junction is such that E J /E c > 25, suppressing unwanted sensitivity to the offset charge of the superconducting island, justifying its absence in the Hamiltonian of Eq. ( 1) [59].
In order to perform microwave spectroscopy measurements, the transmon is capacitively coupled to a readout resonator which is in turn coupled to a transmission line.This allows us to measure the circuit's complex microwave transmission S 21 through the transmission line's input (1) and output (2) ports.
We implement the circuit as shown in Fig. 2(b-d).The device differs from conventional circuit QED geometries in several ways [62], in order to allow the application of magnetic fields in excess of 100 mT.Apart from the Josephson junctions, all circuit elements are made out of field compatible 20 nm-thick NbTiN films [63].We additionally incorporate vortex pinning sites in the ground plane, the transmission line, the resonator island and the transmon island [64].We use a lumped element readout resonator, which has previously been successfully utilized in flux-sensitive devices up to 1 T [65].Its capacitance is formed by an interdigitated capacitor to ground, while its inductance is formed by a 200 nm wide NbTiN nanowire, which has a kinetic inductance of 15 pH/ .This design localizes the regions of high current density at the narrow inductor where vortices are less likely to nucleate due to its reduced width [66].For the transmon circuit the SQUID loop area is chosen to be small, ∼ 5 µm 2 , in order to suppress flux noise from misalignment in large parallel magnetic fields.Finally, InAs/Al nanowires, in which both junctions are defined, have been shown to support sizeable Josephson energies in fields in excess of 1 T [48,65].Further details about device fabrication as well as the cryogenic and room temperature measurement setup can be found in the Supplementary Information (Section II) [60].

III. ANDERSON MODEL FOR A QUANTUM DOT JUNCTION
As we will show, the quantum dot junction can be described by a single Anderson impurity tunnel-coupled to two superconducting leads (Fig. 1(b)).We review its most important properties to facilitate the discussion of the experimental results that follow.
The Hamiltonian of the model takes the form The first term describes a single-level quantum dot, Here, ↑,↓ = ± E Z /2 gives the single-particle energies: is the dot energy level measured with respect to the Fermi level in the leads, which can be controlled via the central electrostatic gate, and E Z = gµ B B is the Zeeman energy.In the latter, g is the effective g-factor of the level, µ B is the Bohr magneton, and B is the magnetic field strength.In the experiment we choose the B-field direction to be parallel to the nanowire in order to maximize the magnetic field compatibility.Finally, U > 0 is the repulsive Coulomb interaction between the electrons, which disfavors the double occupancy of the impurity, while n σ = d † σ d σ are number operators for the dot level, with d σ (d † σ ) the electron annihilation (creation) operators.The resulting energy diagram described by Eq. ( 4) is shown in Fig. 1(c).
The many-particle energy levels of Eq. ( 4) are divided in two sectors, corresponding to their fermion parity, or equivalently, to their total spin S. The singlet sector includes the states of even parity, which have S = 0: the empty state |0 and the pair state The second term in Eq. ( 3) describes two superconducting reservoirs, where i = L, R labels the left and right leads, k labels spin-degenerate single-particle states, and ∆e −iφi is the s-wave pairing potential in each reservoir.The gaugeinvariant phase difference between them is It is made experimentally tunable with magnetic flux as discussed in section II.We assume the reservoirs to have identical gap ∆ and density of states ρ; this assumption should be reasonable since in the experiment the two leads are made out of a single hybrid nanowire.We further take the g-factor of the reservoirs to be zero, capturing the magnetic field dependence of the combined system in the effective quantum dot g-factor of Eq. ( 4).Finally, the third term is the tunneling Hamiltonian coupling the dot and the reservoirs, where t i are the dot-reservoir tunnel coupling strengths, which, for simplicity, we choose to be independent of k and spin.The tunneling rate across each barrier is given by Γ i = πρ |t i | 2 .In practice these rates can also be tuned with electrostatic gates.The tunneling terms in H T break the conservation of the parity and spin in the quantum dot.Nevertheless, the notion of singlet and doublet sectors introduced for the dot Hamiltonian of Eq. ( 4) is inherited by the total Hamiltonian of Eq. ( 3), provided that the spin S is now regarded as the total spin of the system, including that of quasi-particles in the reservoirs.
Over the years, the model of Eq. (3) (or immediate extensions of it) has become paradigmatic to describe quantum dots coupled to superconducting leads.It has been studied in different limits and using a variety of numerical methods, often requiring advanced many-body methods such as NRG and quantum Monte Carlo for full quantitative descriptions [10].
A salient feature of the model is that a quantum phase transition between doublet and singlet ground state can occur upon changing several experimentally-tunable parameters.The dot energy level ξ, the coupling strengths Γ L,R , as well as the superconducting phase difference φ and the magnetic field B all act to shift the relative positions of V s and V d and cause an energy crossing between the ground states of the two sectors.In the measurements reported in Sec.IV and V, we vary all these parameters and compare the extracted phase boundaries to theory.
For the theoretical comparison we use the NRG method [55,67,68] to compute the lowest-lying eigenvalues in the singlet and doublet spin sectors for the Hamiltonian in Eq. ( 3) as a function of the phase difference φ.This results in the Josephson potentials V s (φ) and V d (φ), which are then used as input to the model of Eq. ( 2) to calculate the transmon transition frequencies [60].The projection onto the lowest-energy state of the Josephson junction in each sector is enough to capture the salient features of our experiment, although the inclusion of excited Andreev states of the quantum dot junction in the circuit model is theoretically possible [35,50,[69][70][71].
Experimentally, the observation of the phase transition is facilitated by the presence of a π-phase shift between V s (φ) and V d (φ).The phase shift originates from the required permutation of a spin-up and a spin-down electron when a Cooper pair sequentially tunnels through a dot that initially is occupied by one quasiparticle [53].Thus, while V s (φ) has a minimum at φ = 0, as encountered for conventional Josephson junctions, V d (φ) has a minimum at φ = π (Fig. 1(d)).Therefore, a quantum dot junction in a doublet state is often denominated as a π-junction, and the singlet-doublet transition is also referred to as the 0-π transition.In the following sections we will use the presence or absence of such a π-phase shift to identify regions with a singlet or a doublet ground state [72].

IV. TRANSMON SPECTROSCOPY OF THE QUANTUM DOT
To perform spectroscopy of the resonator, we monitor the microwave transmission S 21 across the transmission line while varying the frequency of a single continuous microwave tone, f r .This results in a dip with Lorentzian lineshape around the resonance frequency of the lumped-element resonator.Two-tone spectroscopy is subsequently performed by fixing the frequency of this first tone, f r , at the minimum of the transmission amplitude, |S 21 |, while varying the frequency of a second tone, f t , also sent through the transmission line.When the second tone matches the frequency of the ground to first excited transmon transition, f t = f 01 , a peak in |S 21 | is observed due to the transmon-state-dependent dispersive shift of the resonator [54].This gives us access to the transmon transition frequency.
We are interested in the behavior of the device when a single level of the quantum dot provides the dominant contribution to the Josephson effect [72].To find such a regime, we search for an isolated resonance in the gate dependence of the frequency spectrum.Isolated resonances often occur when the gate voltages controlling the quantum dot are set close to their pinch-off val-ues [46,47], here operationally defined as the voltage values below which the quantum dot junction does not contribute appreciably to the transmon's transition frequency.In order to identify the right gate configuration, we perform the following sequence of calibration measurements.First, we characterize the reference junction with the quantum dot pinched-off; second, we explore the sizeable parameter space governed by the three quantum dot gates; third, we identify the relation between B y and φ ext through the transmon frequency's SQUID oscillations; and finally we define appropriate gate coordinates to account for cross-couplings.These calibration measurements are detailed in the Supplementary Material (Section III) [60].As a result of this procedure, the gate voltage of the reference junction V j is fixed such that the transmon frequency when the quantum dot junction is pinched-off is f 0 01 ≈ 4.4 GHz.Furthermore, we fix V l = 470 mV and introduce virtual plunger (V p ) and right tunnel (V t ) gates as a linear combination of V c and V r , such that in what follows ξ is mostly independent of V t .
We then move on to study the quantum dot junction.We first monitor the resonator frequency for φ ext = 0 while the plunger gate voltage V p is varied (Fig. 3(a)).This reveals a resonant shape which is discontinuously interrupted near its peak at V p = 395 mV, followed by other discontinuous jumps in the resonator frequency.A zoom into the resonance is shown in Fig. 3(b) and the corresponding transmon transition frequency, exhibiting the same discontinuity as the resonator, is shown in Fig. 3(d).We identify regions in V p where the transmon frequency f 01 is larger and smaller than the reference frequency f 0 01 .This hierarchy is reversed upon changing the applied flux to φ ext = π, as shown in Figs.3(c, e).
These observed discontinuities in frequency are a signature of a singlet-doublet transition.The change of the ground state of the quantum dot junction determines a sudden switch in the branch of the Josephson potential of Eq. ( 2) (from V s to V d or vice-versa) and, thus, a sudden change in the transmon frequency.This is illustrated theoretically in Figs.3(f-g), which show the expected evolution of the transmon frequencies as a function of ξ.
Here, the transition occurs as the single-particle energy level is tuned toward the electron-hole symmetry point ξ = 0, where the doublet ground state is energetically favorable.The resemblance of the dispersion of the transition frequencies in Figs.3(f-g) to the experimental data in Figs.3(d-e) confirms that V p primarily tunes ξ, as was intended with the virtual gate voltage definition.
The occurrence of the singlet-doublet transition requires a change of the fermion parity of the quantum dot junction.In the S-QD-S setup, this is possible in the presence of a population of excited quasiparticles in the superconducting leads, providing the required fermion parity reservoir.The presence of these quasiparticles should further result in a finite occupation of both the singlet and doublet states when their energy difference is small compared to the effective temperature of the quasiparticle bath, namely in the vicinity of the transition.Indeed, 01 , the transmon frequency set by the reference junction when the quantum dot is pinched off.(e) Same as panel (d) but for φext = π.For panels (a-e) Vt = 182 mV and V l = 470 mV.(f) Theoretical estimates of the singlet (orange), doublet (purple) and reference junction-only (dotted, grey) transmon frequencies as ξ is varied for φext =0.Solid (dashed) lines indicate which quantum dot occupation corresponds to the ground (excited) state.(g) Same as panel (f) but for φext = π.For panels (f-g) ∆/h = 46 GHz, U/∆ = 12.2, ΓL/∆ = 1.05 and ΓR/∆ = 1.12.ft and f01 denote, respectively, the frequency of the second tone in two-tone spectroscopy and the first transmon transition frequency (see text).
upon closer inspection of the data of Figs.3(b-c), both branches of the spectrum are visible in a small frequency window surrounding each discontinuous jump.This is because these transition spectra are obtained by averaging over many subsequent frequency sweeps, thus reflecting the occupation statistics of the junction.This feature is further discussed in the next Section.
In Figs.3(d-e), the fact that the frequency shift of the transmon has the opposite sign for the singlet and doublet sectors is a consequence of the π-phase shift in the Josephson potential between the two sectors.For the case φ ext = 0, the singlet potential interferes constructively with the reference junction potential, while the doublet potential interferes destructively, resulting in f 01 > f 0 01 for the singlet and f 01 < f 0 01 for the doublet.This behaviour is reversed when φ ext = π, and thus serves as a method for identifying the quantum dot junction state.

V. SINGLET-DOUBLET TRANSITION BOUNDARIES
Having established a method for identifying singlet and doublet states by transmon spectroscopy, we now experimentally investigate the phase diagrams of the quantum dot junction.We focus on the behaviour around V p = 395 mV and monitor singlet-doublet transitions versus multiple different control parameters.

A. Plunger gate and Flux
We first study the singlet-doublet phase map in V p and φ ext space.Fig. 4(a) shows the transmon frequency offset with respect to the frequency set by the reference junction, ∆f 01 = f 01 − f 0 01 , as a function of V p and φ ext .As discussed in the previous Section, positive values of ∆f 01 result from constructive interference between the two junctions, while negative ∆f 01 values result from destructive interference.Going from left to right, three distinct plunger regions can be observed, with a sudden flux offset of exactly π between them (Fig. 4(b,d)).Based on the preceding discussion of Fig. 3, we identify the outer two regions as phases with a singlet ground state and the inner region as a doublet ground state.We note that the change in contrast between the two singlet regions suggests that V p also weakly tunes Γ L,R in addition to ξ.
For values of V p close to the singlet-doublet transition we also observe a sinusoidal dependence of the transition boundary on the external flux, resulting in an enhanced region of doublet occupation around φ ext = π with respect to φ ext = 0.This comes about from interference between tunneling processes involving the two superconducting leads of the quantum dot junction [4,29], as further discussed in Sec.V B. At a value of V p fixed near this boundary one thus also observes a singlet-doublet transition versus the external flux (Fig. 4(c)).
In Fig. 4 and subsequent figures, the transition boundary between the singlet and doublet phase appears to be sharp and not affected by the thermal broadening typical of transport experiments [20][21][22][28][29][30].The sharpness is a result of a selective spectroscopy technique.As detailed in the Supplementary Information (Section III.E) [60], in the vicinity of a transition two resonant dips appear in single-tone spectroscopy, one for the singlet and one for the doublet.In this circumstance, the center frequency of either dip can be chosen as the readout frequency for the subsequent two-tone spectroscopy measurement.This binary choice selects the transmon transition frequency belonging to the corresponding quantum dot junction state.It is reasonable to assume that the most prominent dip corresponds to the state of the quantum dot junction which is more prominently occupied, and thus lower in energy.If this is the case, the extracted phase boundaries are a close approximation of the zero-temperature phase diagram of the quantum dot junction.When the occupations of singlet and doublet states are almost equally probable, the selective spectroscopy method is affected by selection errors, which leads to the pixelation effects visible in Fig. 4 near the phase boundaries.In Sec.VI, we will explicitly measure the lifetimes of the quantum dot in the singlet and doublet states, substantiating the latter statements.

B. Tunnel gate
Next, we explore the singlet-doublet transition in plunger and tunnel gate space, where the tunnel gate is expected to control Γ L,R .Fig. 5(a) shows ∆f 01 versus plunger and tunnel gates at φ ext = 0. We can identify the region where ∆f 01 > 0 as the singlet phase and the region where ∆f 01 < 0 as the doublet phase.The region of doublet occupancy takes the shape of a dome, similar to the one coarsely seen in flux-insensitive tunneling spectroscopy experiments [14,17].This shape is in accordance with theoretical expectations for the boundary in the ξ − Γ plane.Its physical origin depends on the parameter regime [9].For U ∆ it arises due to an increase in induced superconductivity on the dot with increasing values of Γ, favoring BCS-like singlet occupation.For U ∆ it instead comes about from increased anti-ferromagnetic Kondo exchange interactions between the spin on the dot and the quasiparticles in the leads, favoring a Yu-Shiba-Rusinov (YSR)-like singlet occupation.In both regimes the singlets compete with doublets, ultimately determining the transition to a singlet ground state at large enough Γ = Γ L + Γ R .
We investigate the same plunger and tunnel gate dependence at an external flux φ ext = π, see Fig. 5(b).We find that the doublet phase is enhanced considerably compared to φ ext = 0, due to the previously mentioned interference between tunneling processes to the superconducting leads.Notably, rather than a dome-like shape, the phase boundary takes a characteristic "chimney" shape that was theoretically predicted [4] but, to our knowledge, not yet confirmed experimentally before these measurements.Unlike the dome, the chimney does not close for any V t .In an extended gate range, it is seen to connect to another doublet region of the parameter space which was disconnected from the dome of Fig. 5(a) at φ ext = 0 (Supplementary Information Section III.D [60]).
The chimney at φ ext = π is much less thoroughly researched than the dome at φ ext = 0.The open questions include that of the exact nature of the doublet states as a function of the U/∆ and Γ/U ratios, and the role of the flux bias [73][74][75][76].In particular, when U ∆, the doublet state for small Γ is a decoupled doublet state with a single local moment in the quantum dot.On the other hand, in the same limit but at large Γ (i.e. in the neck of the chimney), the strong exchange interaction with both superconductors is expected to lead to some mixing with the doublet states that involve one Bogoliubov quasiparticle from each lead [77], causing an overscreening of the local moment in the quantum dot.The role of the exchange interaction is more pronounced at φ ext = π also because the anomalous component of the hybridisation (describing the proximity effect) is suppressed due to the cancellation of contributions from the left and right leads [75], where the cancellation is exact when Γ L = Γ R .This further stabilizes the spin-doublet states.The experimental observation of the chimney calls for more thorough theoretical studies of this parameter regime of the model.
We compare the results at both values of external flux to the expected transition frequencies obtained from NRG calculations using Eq. ( 3).We assume that ξ = 0 at V p = 395 mV since this is the symmetry point of the experimental data.At this point, by requiring simultaneous agreement between experiment and theory for both values of external flux (Fig. 5(c)), we are able to extract several of the model parameters.We find that ∆/h = 46 GHz (190 µeV), close to the bulk value of Al.We furthermore extract U/∆ = 12.2, corresponding to a sizeable charging energy of 2.3 meV.It places the nature of the singlets near ξ = 0 in the strongly correlated regime, with a YSR-like character rather than a BCS-like one.By matching values of Γ L,R to V t we then find that Γ/U varies between 0.05 and 0.4, while Γ R /Γ L ≈ 0.75 − 1 in the range of gates explored (Fig. 5(d)).The details of the numerical procedure as well as error estimation can be found in the Supplementary Information (Section I.C), including estimates based on an alternative potential shape for the reference junction [60].
The extracted set of parameters is consistent with the observed dome shape at φ ext = 0, as shown in Fig. 5(e).Additionally, as a result of the ratio Γ R /Γ L remaining close to 1, the extracted parameters also match the observed diverging behaviour at φ ext = π (Fig. 5(f)), which was not enforced in the parameter extraction.We note that in these panels we did not map V p to ξ beyond identifying V p = 395 mV with ξ = 0.A unique mapping could not be constructed due to the unintended dependence of Γ on V p .We speculate that this causes the remain- ing discrepancies between the measured and calculated boundaries in the horizontal direction.

C. Parallel magnetic field
Finally, we investigate the effect of a parallel magnetic field on the phase transition boundaries.Here, we expect a magnetic-field induced singlet-doublet transition to occur [14,18,19].As B z increases, the doublet sector separates into spin species that are aligned and antialigned with respect to the magnetic field, dispersing in opposite energy directions.The singlet ground state energy, on the other hand, is approximately independent of magnetic field.Given an appropriate zero-field energy level configuration, for some B z value the energy of one of the two doublet states will thus become lower than that of the singlet, and become the ground state instead (see Fig. 1(c)).
Such a transition will only occur for specific configurations of V p and V t in the experimentally accessible range of magnetic fields.We therefore start by applying B z = 200 mT parallel to the nanowire axis, a sizeable magnetic field, yet one for which the E J of the reference junction is not yet substantially suppressed.At this field we investigate the effect on the V p and V t phase map.The result, shown in Figs.6(a-b), reveals an expansion of the doublet region for both φ ext = 0 and φ ext = π.We can classify different regions in the parameter space by comparing the phase boundaries at B z = 10 mT and B z = 200 mT.There are regions in which a singlet ground state remains a singlet ground state, independent of the flux and the magnetic field, as well as regions where a singlet-doublet transition occurs depending on the value of the flux.There is also a region that starts off as a singlet ground state and ends up as a doublet ground state at high field, for all values of the flux.Thus, fixing V p and V t in this region, we expect to observe a transition with B z for any value of φ ext .A measurement of ∆f 01 versus φ ext and B z (Fig. 6(c)) indeed reveals such a transition, occurring at a different magnetic field depending on the flux value.For details about the data analysis and identification of the flux axis we refer to the Supplementary Information (Section IV) [60].

VI. DYNAMICS OF THE SINGLET-DOUBLET TRANSITION
In the preceding sections we made use of selective spectroscopy to reconstruct the phase transition boundaries.We now turn to time-resolved spectroscopy techniques to study the parity dynamics of the quantum dot junction close to the transition, aiming to characterize the lifetimes of singlet (even parity) and doublet (odd parity) states.These methods have previously been used to study quasiparticle dynamics in superconducting qubits [48,78], and recently also applied to a nanowire junction to study the poisoning of Andreev bound states [33,45,79].
To resolve individual switching events we use a second device (device B) with a larger signal-to-noise ratio (SNR) than the device used for the preceding sections (device A), enabling the use of short acquisition times.Device B is nearly identical to device A, except for two features meant to increase the SNR: (1) a stronger coupling between the resonator and the transmission line; (2) an additional capacitor at its input port, which increases the directionality of the outgoing signal [80].On device B we perform measurements on microsecond timescales by directly monitoring changes in the outgoing signal at a fixed readout frequency.A continuous measurement of the outgoing microwave field then reveals a random telegraph signal between two different levels, a consequence of the switches in the quantum dot junction parity (Fig. 7(a-b)).Owing to the increased temporal resolution of the detection method, even short-lived excited state occupation of a few µs can now be detected.The characteristic time scales of the telegraph signal reflect the underlying lifetimes of the singlet and doublet states, T s and T d , or equivalently their decay rates, Γ s = 1/T s and Γ d = 1/T d .These quantities can be extracted via a spectral analysis of the time traces, as described in detail in the Supplementary Information (Section V) [60].
To investigate the switching dynamics we tune device B to a regime similar to that of Sec.V A studied in device A. By measuring S 21 with single-tone spectroscopy we once-more find ground state transitions between singlet and doublet as a function of V p (Fig. 7(c)).The discontinuous resonant shape, akin to that of Fig. 3(b), is symmetric around V p = 546 mV, which we identify with ξ = 0.The singlet and doublet resonant frequencies are simultaneously visible close to the discontinuity at the transition.The time-resolved measurements over the same gate voltage range reveal a smooth but strong evolution of the parity lifetimes with V p (Fig. 7(d)).The hierarchy of lifetimes inverts as V p is tuned across the phase transition, reflecting the change in the ground state parity.Away from the transition, in either the singlet or doublet phase, we observe the lifetime in the ground state sector to be on the order of several milliseconds, exceeding that of the excited state by more than an order of magnitude.These numbers are very favorable for the implementation of Andreev pair qubits [32] as well as Andreev spin qubits [33,34,45], whose control has so far been limited by microsecond parity lifetimes .
We further explore the evolution of the relative lifetimes versus V p and φ ext .Fig. 7(e) shows a twodimensional map of log 10 (T d /T s ), which is a measure of the lifetime asymmetry.We find behaviour similar to that previously seen in Fig. 4, with a sinusoidal boundary of equal rates, indicative of the singlet-doublet transition.Furthermore, we observe a strong polarization of the junction parity inside the doublet phase (T d T s ), where the signal-to-noise ratio (SNR) eventually becomes limited by our ability to resolve the rare and short-lived switches out of the ground state.Additionally, we find a modulation of T d with flux, with longer lifetimes at φ ext = π [60].This flux dependence likely originates from the oscillation of the singlet-doublet energy gap with flux, but might also be indicative of a coherent suppression of the tunneling rates, as previously observed in Al tunnel junctions [81].The polarization of the junction parity also occurs inside the singlet phase, where T s T d for V p values away from the transition (Fig. 7d).
Strong parity polarization may not be surprising for a system in thermal equilibrium at temperatures below 100 mK, typical of these experiments, corresponding to a thermal energy small compared to the singlet-doublet energy difference away from the transition.However, parity lifetimes in Josephson junctions are seldom determined by thermal fluctuations, but rather by highly energetic non-equilibrium quasiparticles often observed in superconducting circuits [82].While non-equilibrium quasiparticles are most likely also present in our device, we believe that their influence is suppressed by the large charging energy of the quantum dot junction.
Finally, we observe a non-monotonic variation of the rate asymmetry inside both the singlet and doublet phase, forming apparent contours of fixed lifetimes.We hypothesize two possible reasons, distinct in origin but possibly co-existing, behind this structure in the data: it could be caused by parity pumping mechanisms where the readout tone is resonant with the energy difference between singlet and doublet [79], as well as by the spectral density of the non-equilibrium quasiparticles present in the environment [83].Further investigation of the tunnel gate, power, and temperature dependence of the rate asymmetry can be found in the Supplementary Information (Section VI) [60]; we leave a more detailed study for future work.

VII. CONCLUSIONS
We have demonstrated the use of a transmon circuit to sensitively detect the ground state parity of a quantum dot Josephson junction.The transition frequency of the transmon exhibits a discontinuity if the ground state of the device changes from a singlet to a doublet, due to the presence of a π-phase shift in the Josephson potential of the junction.This allowed us to accurately reconstruct the occurrence of the singlet-doublet transition as a function of all control parameters available in a single device, matching them to those expected from NRG calculations of an Anderson impurity model.In particular, we have observed the flux-induced enhancement of the doublet phase, in the form of the striking transformation of a dome-shaped phase boundary at φ ext = 0 into a chimney-shaped phase boundary at φ ext = π (Fig. 5).
In future research, this singlet-doublet tuning capability could become beneficial for several applications.First, it can be used to define and control Andreev pair and spin qubits, and to couple them to conventional superconducting qubits.Second, tuning the dot to the doublet phase is a robust way to induce a π phase shift, which could be exploited to define a hybrid 0−π qubit that does not rely on the fine-tuning of the applied flux [84].Third, it can facilitate the bottom-up realization of a topological superconductor from a chain of proximitized quantum dots [85][86][87].Finally, fast gate or flux-based switching between the 0 and π shift of the dot can also be of interest for applications in Josephson magnetic random access memory (JMRAM) technologies [88].
We have subsequently used continuous time-domain monitoring of the transmon resonant frequency to determine the lifetimes of singlet and doublet states.We find that the time between switching events is strongly enhanced when the quantum dot is tuned away from the phase transition.Since our estimates indicate that U ∆ in our devices, we attribute this effect to the large energy difference associated with charging the quantum dot.These findings are encouraging for Andreev qubits, which benefit from long parity lifetimes, and suggest that large-U quantum dots could be effective as filters for high-energy quasiparticles.However, further work is required to understand the full dependence of parity lifetimes on U .
In this work we have focused on the study of a singlelevel quantum dot by tuning our junction very close to pinch-off.Looking forward, there is much left to explore in the parameter space of such a device.To begin with, it would be interesting to understand whether the crossover from the BCS-like to the YSR-like singlet has any signature in the microwave response of the system.Second, opening the junction further brings the quantum dot into a multi-level regime, not captured by the single impurity Anderson model, and still largely unexplored.Finally, while we have primarily studied the ground state properties of the quantum dot junction, microwave spectroscopy should allow to study its excitations, as e.g. recently demonstrated in Refs.[49,52], particularly at φ ext = π.Further work will also aim at elucidating the role of spin-orbit coupling in the quantum dot junction.It is well known that, when time-reversal invariance is broken, spin-orbit coupling can induce a spin-splitting of energy levels in the doublet sector [36,37,44], essential for Andreev spin qubits.While this effect could have been expected to occur in the measurements presented here, it was not detected; we speculate that the level spacing in the dot was too large to result in a significant splitting [37].
Important extensions of our work could arise if the hybrid nanowire in our microwave circuit was driven into the Majorana topological phase [69][70][71]89], which is currently challenging because of a large parameter space [90] and because of demanding disorder requirements [91].Including a quantum dot in a Josephson junction between two topological superconductors could be beneficial for the detection of the 4π Josephson effect: as we have seen, it mitigates quasiparticle poisoning, although it would not resolve [92] the problem of distinguishing Majorana zero modes from trivial zero-energy Andreev bound states [93].Finally, the manipulation of quantum dots coupled to superconducting leads is an essential ingredient of scalable proposals for topological quantum computation [40]. .

I. NUMERICAL MODELING
As discussed in the main text, we model the quantum dot junction as a single Anderson impurity coupled to two superconducting leads.Using numerical renormalization group (NRG) methods we extract the energies of the singlet and doublet states for any combination of the model parameters, and subsequently incorporate them into a DC SQUID transmon Hamiltonian.This Hamiltonian is then used to match the experimental data and extract estimates of the model parameters.These procedures are detailed in this section.

A. NRG calculation
The NRG method is an iterative procedure for solving quantum impurity problems involving a localized few-level system coupled to a continuum of itinerant electrons (fermionic bath, normal-state or mean-field BCS superconductor).It consists of several steps: 1) discretization of the continuum parts of the Hamiltonian using a geometric-progression mesh with an accumulation point at the Fermi level (the so-called logarithmic discretization), 2) unitary transformation of the resulting discretized Hamiltonian from the star-geometry (impurity coupling to each representative mesh point) to a linear tight-binding chain representation (the so-called Wilson chain), 3) iterative diagonalization in which the Wilson chain sites are taken into account consecutively [1][2][3][4][5].The discretization is controlled by the discretization parameter Λ > 1 which controls the coarseness of the grid.When the discretization is coarse, the results can be improved by twist averaging, which consists of performing the same calculation for several different discretization grids and averaging the results [5,6].The growth of the Hilbert space is controlled by the truncation parameters which control the number of states retained after each step of the iteration.
The calculations in this work have been performed with the NRG Ljubljana code [7].Since the main quantities of interest are the ground state energies in each spin sector, very high quality results can be obtained even with coarse discretization (Λ = 8) and keeping no more than 3000 states (spin multiplets) in the truncation.We have verified that the twist averaging is not required.The BCS gap was chosen to be ∆ = 0.1D, where D is the half-bandwidth.The calculations were performed for a problem with symmetric hybridisations, Γ L = Γ R .This is sufficient, because the results for an arbitrary coupling asymmetry can be obtained from the following mapping [8]: where a = Γ L /Γ R is the asymmetry, φ is the BCS phase difference in the asymmetric problem, and φ S is the effective BCS phase difference in the effective symmetric problem.Such calculations were performed for a set of values of the interaction strength U (from very low values U = 0.1∆ that correspond to ABS-like subgap states, up to U = 30∆ that correspond to YSR-like subgap states).In every value of U , a grid of ξ and Γ parameters was set up, and a sweep of φ between 0 and π (50 points) has been performed for each (ξ, Γ) pair.The ground state energies are obtained as the sum of all energy shifts [2] performed during the NRG evolution, which has been shown to produce extremely accurate results [6].Some calculations have also been performed in the presence of a small Zeeman splitting.The results have been collected, documented, and made available on a public repository [9].The full set of input files and scripts is provided for running the calculations for different parameters or for different Hamiltonians.
Having developed the NRG calculation, we can gain insight into the expected boundaries between singlet and doublet occupation.In Fig. S1(a), we show the phase diagram for the symmetric configuration Γ L = Γ R at fixed φ = 0 and U/∆ = 5.In the (ξ, Γ) plane, the phase diagram takes a dome-like shape with the transition value of Γ being the highest at the electron-hole symmetry point ξ = 0.At this point, the transition value of Γ diverges if the phase difference between the reservoirs is changed to φ = π, because in this case a destructive interference between tunneling events to the left or right occurs.This causes the "dome" in the (ξ, Γ) plane to turn into the "chimney" shown in Fig. S1(b).
As mentioned above, at Γ = 0 the ground state is in the doublet sector for |ξ/U | < 1/2.Upon increasing Γ, the Kondo coupling favours the binding of a Bogoliubov quasiparticle in the superconductor to the impurity local moment ("Yu-Shiba-Rusinov" screening), ultimately determining the transition to a singlet ground state at a value Γ c .The value of Γ c depends on ξ, φ, U and ∆, as well as on the asymmetry between Γ L and Γ R .This implies that the singlet-doublet transition can be observed varying any of these parameters individually.Since in the experiment the values of U and ∆ are fixed, being determined by the materials and the geometry of the physical device, we focus here on variations in ξ, φ, Γ L and Γ R .
In Fig. S1(c), we show the singlet-doublet transition boundary in the ξ−φ plane.The interference effect is modulated continuously by the value of the phase difference φ, resulting in periodic oscillations of the boundary.The average position of the oscillating boundary is determined by Γ.In Fig. S1(d), we show the effect of a Zeeman energy E Z in the case when the ground state is singlet at B = 0.As mentioned in the main text, a singlet-doublet transition is induced at finite E Z due to the spin-splitting of energy levels in the doublet sector.

B. Transmon diagonalization
Having established how to calculate singlet and doublet potentials using the NRG method, we now turn to their inclusion in the Hamiltonian of the transmon circuit (main text Eq. ( 1)).To numerically solve the Hamiltonian for an arbitrary potential term V (φ) we make use of the Fourier decomposition (note that the potential can include an external flux φ ext ): with the components where we assume the potential to be a real-valued 2π-periodic function.We can then express the full Hamiltonian in the charge basis as with E J,n = E cos J,n − iE sin J,n , N the charge operator and N+ the charge raising operator defined via N+ |N = |N + 1 .Upon substituting the potential of main text Eq. ( 2) into Eq.(S4) and diagonalizing the Hamiltonian, we find the eigenvalues and obtain the energy levels of the combined reference junction and quantum dot junction system.Their difference then results in the transmon's transition frequencies.To numerically compute the eigenvalues we truncate the number of charge states and Fourier coefficients to N = 35 for all calculations [10].We verify that this leads to good convergence for the eigenvalues.We further note that while the presence of the potential offset E J,0 does not affect the transmon transition frequencies, its inclusion is crucial: it plays a large role in determining whether the ground state of the combined system corresponds to singlet or doublet occupation for a given set of quantum dot junction parameters.

C. Parameter matching routine
To match the numerical model to the experimental data we have to overcome several complications.First, the mapping between experimental control parameters and those present in the model is not always trivial.As discussed  in the main text, V p appears to not only tune ξ but also Γ L,R .In turn V t is constructed in such a way that (to first approximation) it does not tune ξ, but it does act on both tunnel rates simultaneously with different, unknown lever arms.For mapping the magnetic field axis to the Zeeman energy the challenge lies in determination of the effective g-factor of the quantum dot, known to be a strongly gate and angle-dependent quantity [11].Only the flux axis allows for a simpler identification, in particular if one assumes that in the singlet configuration the combined DC SQUID Josephson potential takes its minimal (maximal) values at 0 (π), which should hold for even modest SQUID asymmetry.A separate challenge comes from the large number of parameters of the model: ∆, U , ξ, Γ L , Γ R , and φ ext .With 6 potentially correlated parameters to match one has to carefully assess whether the fit is under-determined.
Given these considerations, we identify a specific gate point in the experimental data that could result in a wellconstrained situation: the top of the dome shape of Fig. 5(a) in the main text.Here we have access to three measured quantities at a known flux φ ext : the singlet and doublet qubit frequencies f s 01 (0) and f d 01 (0) measured at the boundary of the transition, and also the doublet qubit frequency f d 01 (π).We furthermore know that here E s ≈ E d for φ ext = 0, since the data lies on the boundary of a singlet-doublet transition versus tunnel gate.Finally, based on the symmetry of the dome shape we identify that this V p should correspond to ξ ≈ 0. We can therefore eliminate two of the model parameters (ξ and φ ext ) and are left to determine ∆, U , Γ L and Γ R .
In a first step we tackle the condition of a singlet-doublet transition occurring versus tunnel gate.For each value of ∆, U and Γ L we numerically diagonalize the Hamiltonian of Eq. (S4) to determine the lowest energy level of the total circuit for both the singlet and doublet states and find the value of Γ R for which these energies are equal.For this we use a reference junction potential V J = E J (1 − cos φ) with E J = 12.8 GHz and E c /h = 210 MHz as determined in Sec.III.Shown in Fig. S2(a), this results in a U -dependent range of Γ L for which there is indeed a value of Γ R that leads to a singlet-doublet transition.Outside of this range Γ L is so large that the ground state is always a singlet.
Having determined these possible values of Γ R we calculate the three relevant transmon frequencies f s 01 (0), f d 01 (0), and f d 01 (π).These are then compared to the measured values, and an optimal solution is sought that minimizes the sum of the absolute difference between calculation and measurement of all three quantities (S5) In Figs.S2(b-d) we plot a sample of this three-dimensional optimization, while Fig. S3 shows how each panel is constructed from the individual singlet and doublet qubit frequencies.Other than the trivial symmetry between Γ L,R , it appears that there is indeed a single region of parameters matching our data.At its global minimum we find ∆/h = 46 GHz (190 µeV), U = 12.2∆, Γ L = 1.19∆ and Γ R = 1.47∆, which results in a precise match to the measured qubit frequencies.
Having determined ∆, U , Γ L , and Γ R at this single point in gate space, we attempt to match the model to the V t axis of the data.To do so we fix ∆ and U to the determined values and for each value of V t find the best set of Γ L,R to match the data.To determine these two parameters we have two measured quantities: up to the transition we have f d 01 (0) and f d 01 (π), and after the transition we have f s 01 (0) and f d 01 (π).This procedure results in good correspondence to the experimental results, as shown in main text Fig. 5(c,d).We note that by construction this captures all the granularity and measurement uncertainty of the experimental data, even though the underlying quantities might have been more smooth.A subsequent procedure that attempts to match V p to ξ did not turn out to be unique, as V p appears to also act on Γ L,R .We therefore leave this mapping undetermined.
The uncertainty in the extracted quantities is affected by several factors.The first is the measurement accuracy; we measure the qubit frequency with MHz-scale precision.Based on numerical evaluation of the model, this precision in qubit frequency should limit the extracted parameter accuracy to several GHz.A more substantial uncertainty comes from the determination of the transmon island charging energy E c , which is typically determined from the transmon transition anharmonicity α = f 12 − f 01 .While the anharmonicity can be measured to high precision, a complication arises from the usage of a nanowire based Josephson junction as the reference junction.Up to now we have assumed its potential to take the form V (δ) = E J (1 − cos δ); that of a conventional superconductor-insulator-superconductor (SIS) tunnel junction governed by many weakly transparent channels.In this case we find that E c /h = 210 MHz, resulting in the parameter estimates given above.However, previous work has found that nanowire-based Josephson junctions are better described by several or even a single transport channel, such that V (δ) = − n ∆ 1 − T n sin 2 δ/2.This change in potential shape can lead to a strong reduction in the anharmonicity, and thus an underestimation of E c when using the SIS potential [12].We therefore also match our reference junction dependence to a single transport channel, which is the most extreme case for a reduction in the anharmonicity, finding good agreement with a single transport channel of T = 0.58.This in turn leads to an extracted E c /h = 306 MHz (see Sec. III), resulting in a different set of extracted quantum dot parameters.In particular, we now find ∆ = 30.5GHz and U = 17.3∆.This value of the induced gap in the InAs-Al nanowire is on the low end of what is typically found in DC transport experiments, which might hint at a reduced proximity effect in the ungated leads [13,14].
Capacitance simulations of the full circuit do not provide an unambiguous answer for which of the two limits is more appropriate, as the circuit was designed to target E c /h = 250 MHz which falls in the center of the estimated range.As it stands we therefore do not have to uniquely determine the experimentally realized E c and thereby resolve the uncertainty in the extracted quantities.However, future works could make use of additional circuit QED compatible quantum dot probes such as direct DC access [15] or dispersive gate-sensing techniques [16] to independently characterize several model parameters and further constrain the matching.

D. Calculated 2D maps
Having established how to match the model parameters to the data, we now turn to the reconstruction of the full 2D dependencies measured in the experiment (Fig. S4).For the plunger versus tunnel gate dependence, we calculate both the singlet and doublet qubit frequencies for all values of Γ L,R encoded by V t for a range of ξ at both φ ext = 0 and π.We subsequently mask the data according to the ground state of the combined transmon Hamiltonian, and obtain a result that closely approximates the measured data (main text Fig. 5).Using the same set of quantum dot junction parameters, we also perform a similar procedure for the 2D map of plunger gate and external flux, resulting good correspondence with main text Fig. 4.

E. State population
We now turn to the singlet and doublet lifetimes determined in device B. For this device we could not identify a measurement point where a unique set of parameters matched the measured data, and can therefore not make a quantitative comparison to the numerics.Instead, we attempt to gain some intuition about the obtained results based on the parameters of device A.
In main text Fig. 7 we extract log 10 (T d /T s ), the ratio of the lifetimes of singlet and doublet occupation.If the system was in thermal equilibrium with a bath of temperature T , one would naively expect that the relative lifetimes should follow the state populations P s,d as described by a Maxwell-Boltzmann distribution: where g i is the degeneracy of the state, E s,d are the singlet and doublet energies, and k B is the Boltzmann constant.
We take ), where we neglect potential other many-body states which should be unoccupied at the experimentally relevant temperatures.In Fig. S5(a) we then plot log 10 (P d /P s ), choosing a bath temperature of 400 mK.Qualitatively this follows the same trend as observed experimentally, with a sharp boundary at the phase transition and a saturated population imbalance away from that.We stress once-more that this is not a quantitative comparison.However, the need for a temperature far in excess of the refrigerator's base temperature of 20 mK could hint at a non-thermal origin such as non-equilibrium quasiparticles [17].
In the main text we also speculate that non-thermal effects lie at the origin of the experimentally observed contours of fixed lifetime ratio's.We corroborate this in Fig. S5(b), where we plot the energy difference between singlet and doublet occupation of the quantum dot junction.This quantity exhibits distinct contours of equal energy difference that qualitatively match those found in the experiment.If the environment has spectral components resonant with these specific energies, one could expect these to modify the dynamics.

II. DEVICE AND EXPERIMENTAL SETUP A. Nanofabrication details
The device fabrication occurs in several steps using standard nanofabrication techniques, and it is identical for device A and B. The substrate consists of 525 µm-thick high-resistivity silicon, covered in 100 nm of low pressure chemical vapor deposited Si 3 N 4 .On top of this a 20 nm thick NbTiN film is sputtered, into which the gate electrodes and circuit elements are patterned using an electron-beam lithography mask and SF 6 /O 2 reactive ion etching.Subsequently, 30 nm of Si 3 N 4 dielectric is deposited on top of the gate electrodes using plasma enhanced chemical vapor deposition and then etched with buffered oxide etchant.The nanowire is then deterministically placed on top of the dielectric using a nanomanipulator and an optical microscope.For this we use an approximately 10 um-long vapour-liquid-solid (VLS) hexagonal InAs nanowire with a diameter of 100 nm and a 6 nm-thick epitaxial Al shell covering two facets [18].After placement, two sections of the aluminium shell are removed by wet etching with MF-321 developer.These sections form the quantum dot junction and the reference junction, with lengths 200 nm and 110 nm respectively.A zoom-in of the the quantum dot junction is shown in Fig. 2(d) of the main text.The reference junction is controlled by a single 110 nm-wide electrostatic gate, set at a DC voltage V j .The quantum dot junction is defined by three 40 nm-wide gates separated from each other by 40 nm, set at DC voltages V l , V c and V r .Note that in Fig. 2(d) the gates appear wider (and the gaps between gates appear smaller) than stated due to distortion by the Si 3 N 4 layer; the given dimensions are therefore determined from a scanning electron microscopy image taken before the deposition of the dielectric.After the junction etch the nanowire is contacted to the transmon island and to ground by an argon milling step followed by the deposition of 150 nm-thick sputtered NbTiN.Finally, the chip is diced into 2 by 7 millimeters, glued onto a solid copper block with silver epoxy, and connected to a custom-made printed circuit board using aluminium wirebonds.

B. General chip overview
Optical microscope images of the chips containing devices A and B are shown in Figs.S6(a) and (b), respectively.Each chip, 7 mm long and 2 mm wide, consists of four devices coupled to the same transmission line.For the chip containing device A, only one device was functional.Out of the other three, one did not have a nanowire, another contained three nanowires stuck together, and for the third device a gate electrode showed no response.The chip of device B includes an on-chip capacitor on the input port of the transmission line to increase the signal-to-noise ratio.For this chip only two of the devices were bonded: device B, which was functional, and another device that did not show any response to the electrostatic gates.The two unbonded devices were dismissed based on prior optical inspection, containing two and no nanowires respectively.

C. Flux control with in-plane magnetic field
In all measurements we control the external flux φ ext with the in-plane component of the magnetic field perpendicular to the nanowire, B y , as illustrated in Fig. S7 [19].This is done since flux tuning with the out-of-plane magnetic field B x led to strong hysteric behaviour in the resonator as well as flux jumps in the SQUID loop.We attribute these effects to Abrikosov vortex generation and the presence of superconducting loops on the chip, causing screening currents.

D. Flux jumps in device A when |B| < 9 mT
For all measurements of device A, the value of the applied magnetic field is kept above 10 mT to prevent flux jumps observed when |B| < 9 mT.In particular, for Figs.3-5 in the main text, B z = 10 mT.The reason for this is purely technical.Device A contains various on-chip aluminium wire-bonds connecting separate sections of the ground plane together.Below the critical magnetic field of aluminium (∼10 mT [20]) these wire bonds create superconducting loops close to the device region, and have a significant cross-section perpendicular to the chip plane.In this regime, the application of an in-plane magnetic field B y generates unwanted currents across these superconducting loops, which in turn result in multiple jumps observed in the flux through the SQUID loop (Fig. S8), making it impossible to reliably control φ ext .Applying a field |B| > 9 mT turns the aluminium wire bonds normal and prevents the unwanted flux jumps, as shown in Fig. S8(a).As this magnetic field is small compared to other energy scales involved, it should not affect the physics under study.We further note that the absence of superconducting loops containing wire-bonds in device B made it possible to measure this device at B z = 0 mT without suffering from similar flux jumps.A three-axis vector magnet (x-axis not shown) is thermally anchored to the 4 K temperature stage, with the device under study mounted at its center.The Bz component of the magnetic field is controlled with a MercuryiPS current source while the Bx and By axes are controlled with Yokogawa GS200 and GS610 current sources respectively.At room temperature a vector network analyzer (VNA) is connected to the input and output RF lines for spectroscopy at frequency fr.On the input line, this signal is then combined with the qubit drive tone at frequency ft for two-tone spectroscopy.A separate tone at fr only used for time-domain measurements is also combined onto this line.For time-domain measurements the output signal is additionally split off into a separate branch and down-converted to 25 MHz to be measured with a Zurich Instruments ultra-high frequency lock-in amplifier.

III. BASIC CHARACTERIZATION AND TUNE UP OF DEVICE A A. Reference junction characterization
In this section we investigate the basic behaviour of the reference junction versus junction gate voltage V j and magnetic field B z when the quantum dot junction is completely closed.This information is used to choose a V j setpoint, V j = 640 mV, which maintains a good SQUID asymmetry in all regimes of interest.Figs.S10(a) and (b) show the V j dependencies of the resonator and transmon frequencies, respectively.As V j is varied, different junction channels open sequentially [21,22], with transparencies that increase non-monotonically due to mesoscopic fluctuations at the junction.This in turn affects the transmon's E J and results in the observed fluctuations of its frequency.
The B z dependencies of f 01 and f 02 /2 at V j = 640 mV are shown in Fig. S10(e).From this we estimate both the transmon island charging energy E c (not to be confused with U , the charging energy of the quantum dot junction) and the parameters of reference junction potential used in Sec.I C to match the measurements to the numerical calculations.Illustrated in this figure is a fit of the data with a Josephson potential governed by a single Andreev level is the field dependent superconducting gap [20], ∆ is the superconducting gap at zero field, B c is the critical magnetic field and T is the transparency of the junction.As the fit is not constrained well enough to provide a unique solution, we fix ∆/h = 60 GHz based on recent experiments on the same nanowires [23].We obtain E c /h = 306 MHz, T = 0.58, and B c = 413 mT, resulting in an effective E J ∼ ∆T /4 = 8.7 GHz.A similar procedure is then performed for We can use these parameters to estimate the experimentally-realized SQUID asymmetry α S = E J /E J,QD where E J,QD denotes the effective quantum dot junction Josephson energy.To do so we estimate E J,QD from the calculated qubit frequencies of the singlet and doublet obtained in Sec.I C through the relation ω 01 ≈ 8E J,QD E c − E c [10].We find that α S > 10 for almost all of the parameter range, exceeding 30 for low values of V t .The asymmetry is at its smallest for the upper values of V t in the vicinity of ξ = 0, where we find a minimum asymmetry α S = 4.We note that the effects of these variations in asymmetry are fully captured by the numerical model; its effects are predominantly on the modulation of the qubit transition frequency with flux and not on the position of the singlet-doublet transition boundaries.

B. Quantum dot junction characterization
In this section we show the basic behaviour of the quantum dot gates when the reference junction is closed.Fig. S11 shows effective pinch-off curves for all three quantum dot gates ramped together (a) and for each of them separately, when the other two are kept at 1250 mV (b-d).This shows that each of the three quantum dot gates can independently pinch off the quantum dot junction even if the other gates are in the open regime, signifying strong lever arms and good gate alignment.We note that these are not pinch-off curves as encountered in conventional tunnel spectroscopy.They reflect the voltages at which there is no longer a measurable transmon transition frequency mediated by the quantum dot junction, which could either be due to low tunneling rates or a full depletion of the quantum dot.

C. Device tune up
This section describes the process of tuning up the quantum dot gates to the setpoint used for the main text figures.We start by closing the reference junction (V j = −200 mV) and going to a point in quantum dot gate voltages near pinchoff (V c = 100 mV, V l = 250 mV and V r = 400 mV, see Fig. S11).Monitoring the frequency of the resonator while varying one of the gates reveals small shifts away from its bare frequency which resemble the shape expected for quantum dot resonances (Fig. S12(a)).Fixing the readout frequency f r at the bare frequency of the resonator, one can map out the regions where these shifts happen on a two-dimensional map versus the left and right gates (Fig. S12(b)).In such maps, a pixel with a dark color indicates the resonator is not shifted from its bare frequency while a bright pixel indicates a shift of the resonator frequency, which we can use to identify potential regions of interest.
After identifying such a region in V l -V r space, we open the reference junction to its set-point V j = 640 mV, which lifts the reference transmon frequency to f 0 01 = 4.4 GHz, closer to the bare resonator frequency.This magnifies the dispersive shift of the resonator and, furthermore, brings the external flux into the picture.As shown in Fig. S12(e), the asymmetric SQUID behaves as expected for different quantum dot gate setpoints.The reference junction sets the reference value for the transmon frequency, f 0 01 , and the quantum dot contributes with small variations above or below this setpoint due to constructive or destructive interference, respectively.Fixing φ ext = 0 and repeating the initial measurement versus V r with the reference junction open reveals much stronger deviations of the resonant frequency than before (Fig. S12(c)).Importantly, the observed resonant frequency is now discontinuous, which, as detailed in the main text, is a signature of a singlet-doublet transition of the quantum dot junction.We tentatively identify the regions for which the resonator frequency is shifted to lower values as doublet regions and perform single frequency readout versus V r and V l , now with f r fixed at the resonator frequency corresponding to doublet regions (Fig. S12(d)).The resulting two-dimensional map reveals regions for which the transmission amplitude signal is low (dark regions in Fig. S12(d)) which we identify as potential regions with a doublet ground state.The next step for tuning up is identifying an isolated region where the quantum dot is in a doublet ground state and exploring the behaviour versus the central quantum dot gate.This is shown in Fig. S13.As V c is varied at φ ext = 0 (Fig. S13(c)), the resonator first shows a displacement towards higher frequencies to then abruptly drop to a lower frequency, to then finally go back to the higher frequencies once-more.As detailed in the main text, we identify this behaviour with a singlet-doublet transition as the relative level of the quantum dot ξ is being varied.Figs.S13(a) and  (b) show how this central doublet ground state region varies with each of the two lateral quantum dot gates.In both cases we observe a dome shape, resembling the behaviour we would expect when varying the tunnel coupling between quantum dot and leads.However, these dome shapes are rotated in V c -V r and V c -V l space.This is understood as the result of cross coupling between the different quantum dot gates.
After identifying the cross coupling effect between different quantum dot gates, we define a new set of virtual gates in an attempt to tune the model parameters independently.We fix V l =470 mV (set-point kept for all results shown in the main text) and focus on V r -V c space.Fig. S14(b) shows the dome shape previously identified in V r -V c space.We identify a line along the dome (indicated with a dashed line) for which the quantum dot level appears to be fixed and define new plunger virtual gate (V p , perpendicular to this line) and right tunnel virtual gate (V t , along this line) (see Fig. S14(d)).This rotated gate frame is the one used for the main text.Note that this routine does not guarantee that V p does not affect the tunneling rates.It rather ensures that V t does not (strongly) affect the quantum dot level ξ.

D. Larger tunnel voltage range
In Fig. S15 we show the behaviour of the singlet and doublet regions beyond the V t range investigated in Fig. 5 of the main text.At φ ext = π we do not observe the doublet phase boundary fully closing for any V t .According to theory, this should only occur if ξ ≈ 0 and Γ L ≈ Γ R are maintained at each gate setting in the experiment.That this condition would remain satisfied for any V t is implausible given the cross-coupling present in the system.We instead speculate that at higher gate voltages the tunnel rates cease to be a monotonically increasing function, which is substantiated by the tunnel gate dependence at φ ext = 0.Here we observe a temporary recovery of the doublet region at higher V t , which should not occur for increasing values of Γ.We further speculate that in this regime of increasingly large Γ/U the dot can eventually be tuned to a different charge configuration, involving energy levels not captured by the single-level model of main text Eq. ( 6).
We note that for these measurements only single tone spectroscopy was performed.We therefore plot ∆f r = f res − f 0 res , where f 0 res denotes the resonator frequency with the quantum dot junction pinched off.Its qualitative interpretation is the same as that of ∆f 01 used in the main text.

E. State selective spectroscopy
For the measurements performed close to singlet-doublet transitions, single-tone spectroscopy simultaneously shows two resonances whose relative depth varies with the distance from the transition.This is once more illustrated in Fig. S16, which shows single-tone spectroscopy at several different V p regions while φ ext is varied.It corresponds to the measurements of Fig. 4 of the main text.In panels (a) and (d) we observe only a single resonance; at these plunger gate values the quantum dot junction is sufficiently deep in the singlet and doublet parity sector respectively that only one state is occupied.However, at the plunger gate values between these two regimes (panels (c-d)) the behaviour is more complex.We simultaneously observe two resonances and their depth becomes a function of flux.
For the two-tone spectroscopy measurements in the main text we make use of the averaged occupation of the states captured in the single-tone spectroscopy measurement to identify most occupied state.This can be inferred from the relative depth of the resonances: for example in Fig. S16(e) the most occupied state is the singlet, albeit by a small margin.This in turn allows us to do state selective two-tone spectroscopy, revealing the transmon transition that  corresponds to the most occupied state of the system.To do so we fix the frequency of the first tone f r at the bottom of the deepest resonance, corresponding to the most populated sector of the system.We illustrate this in Figs.S16(  and (g), where by fixing f r at the bottom of the resonance corresponding to the singlet (doublet) state we observe a peak only when f t is equal to the transmon frequency corresponding to the singlet (doublet) state.It is this peak position that we report as f 01 .

IV. MAGNETIC FIELD DEPENDENCE OF DEVICE A
In this section we elaborate on the analysis of the data shown in Fig. 6(c) in the main text.When varying both φ ext and B z in a measurement, one has to consider the possibility of an unwanted misalignment of the magnetic field with respect to the nanowire axis.This, in combination with the multiple orders of magnitude difference between the applied B z (hundreds of mT) and the B x (less than a µT) or B y (several mT) needed to thread a flux quantum through the SQUID loop, can result in big changes of the φ ext = 0 point for different values of B z .Therefore, one has to re-calibrate the value of B y that corresponds to φ ext = 0 for each B z value.To do so, we use the flux dependence of f 01 at a gate point for which the quantum dot junction ground state remains a singlet for the whole B z range as a reference for identifying φ ext = 0.This gate point is indicated with a grey cross in Fig. S17(a).
The measurement shown in Fig. 6(c) is therefore performed as follows: for each B z value do apply B z for each B y value do apply B y measure f 01 at the grey gate point measure f 01 at the green gate point For each B z value we then reconstruct the B y dependence of φ ext through the dependence of the reference gate point (grey).Furthermore, we use this method to identify points in B y where flux jumps happen and correct for them.While they almost never occur for small magnetic fields, and none of the other data required such a correction, we found that at increasing B z jumps would occur more often.We believe this is due to a small misalignment between B z and the plane of the chip.The resulting corrected φ ext reference is shown in Fig. S17   In this section we elaborate on the analysis method for extracting the characteristic lifetimes of the singlet and doublet states, T s and T d .We start with a continuous measurement at a fixed readout frequency where we monitor the demodulated output signal integrated in time bins of t int = 2.3 µs.This reveals a complex random telegraph signal jumping between two states in the (I,Q)-plane.The histogram of the acquired (I,Q) points shows two states (Fig. S18(a)) whose centers define an axis X.A segment of the measured telegraph signal, projected onto this X axis, is shown in Fig. S18(c).Taking the histogram along this axis results in a double Gaussian distribution (Fig. S18(d)) that is well-described by From the time domain information of the signal we construct its power spectral density (PSD), which is its squared discrete Fourier transform (Fig. S18(b)) where X(t) is the measured signal (as projected onto the previously defined X-axis), ∆t = 2.3 µs is the discrete time bin in which the data is measured, N = T ∆t is the number of points and T is the total signal length.In practice we use Welch's method with a Hanning window [24] to calculate the power spectral density, dividing the trace into 50 sections of length 40 ms that overlap by 20 ms and averaging the power spectral density of all segments.This results in a spectrum that is well fit by a single Lorentzian of the form Combining the amplitude ratio R = A 1 /A 2 obtained from the Gaussian fit of the two quadratures and the Γ value obtained from the Lorentzian fit of the PSD, we calculate Fig. S19 shows the flux dependence of the lifetimes of the singlet and doublet states at V p = 554.4mV, which accompanies main text Fig. 7.We find that both singlet and doublet lifetimes show an approximate sinusoidal dependence on the applied flux.As discussed in the main text, this flux dependence most likely originates from the oscillation of the singlet-doublet energy gap with flux.However it could also be indicative of a coherent suppression of the tunneling rates [25].We further note that the sudden drops in SNR are due to crossings of the transmon frequencies of the singlet and doublet states.At these points both resonator frequencies become indistinguishable and their lifetimes can not be measured.

B. Power and temperature dependence of parity lifetimes
Here we present additional data on the readout power and temperature dependence of the parity lifetimes shown in Fig. 7 of the main text.The power dependence at four selected points across a phase boundary is shown in Figs.S20(cf).Away from the transition (purple) and right on top of the transition (green) the readout power does not have a strong effect on the extracted lifetimes in the investigated range.For plunger gate values V p closer to the transition, however, the asymmetry of the lifetimes decreases with power (blue).Although the origin of this dependence is not clear, we conjecture it is related to parity pumping effects [26].
Temperature dependencies at the same gate points, measured at a readout power of -22 dBm at the fridge input, are shown in Figs.S20(g-j).Here the mixing chamber temperature of the dilution refrigerator is measured with a ruthenium oxide resistance thermometer and increased in a controlled step-wise fashion with a variable-output heater mounted on the mixing chamber plate.We observe different effects of temperature for each of the gate points.In general, there is a temperature independent regime at low temperatures, followed by a temperature dependent drop above a certain characteristic temperature, which varies over tens of mK for different gate points.For some of the gate points, however, the temperature independent contribution is absent and the effect of increased mixing chamber temperature starts immediately at base temperature (Fig. S20(i)).These results are indicative of non-equilibrium effects playing a role in the physics of the devices under study, their exact behaviour dependent on the energy level configuration of the quantum dot junction.

C. Parity lifetimes versus tunnel gate
To complement the data shown in Fig. 7 of the main text, taken at V t = −60 mV, we also show the V t dependence of the parity lifetimes at φ ext = 0 in Fig. S21.As for device A, the doublet ground state region exhibits a dome shape in V p and V t space, and at the transition between singlet and doublet ground states the lifetimes for both states become equal.Away from the transition, the lifetime asymmetry increases and the lifetimes differ by more than one order of magnitude.We note that the gate compensation of device B was not ideal, resulting in a small tilt of the dome.
Similarly to the behaviour shown in the main text for φ ext and V p , in this case we also observe contours of equal ratio where the lifetime asymmetry abruptly increases or decreases.For higher readout power these contours become accentuated, as shown in Fig. S21(c).Furthermore, for higher power the region with similar lifetimes around the ground state transition becomes wider.This is due to the parity lifetimes having a different dependence on power for different regions in gate space.For most regions in gate space there is again almost no dependence on readout power in the range explored (Fig. S21(e,f)).However, on special gate regions, such as close to ground state transitions (Fig. S21(g)) and on top of the observed contours (Fig. S21(d)), the lifetime asymmetry decreases rapidly with power, similar to the effect shown in Fig. S20.

FIG. 1 .
FIG. 1.(a) Schematic diagram of a quantum dot junction incorporated into a transmon circuit.The transmon island with charging energy Ec is connected to ground by a SQUID formed by the parallel combination of a quantum dot junction and a reference junction.In this panel φ and δ denote the superconducting phase difference across the quantum dot and reference junctions respectively.Φext is the externally applied magnetic flux through the SQUID loop.(b) Model diagram of the quantum dot junction in the excitation picture.Two s-wave superconductors are connected via tunnel barriers to a single level quantum dot.(c) Level diagram of the quantum dot hosting 0, 1, or 2 electrons when disconnected from the leads (ΓL = ΓR = 0).(d) Phase dependence of the Josephson potential of the quantum dot junction in the singlet (orange) and doublet (purple) state.(e) Josephson potential of the reference junction.(f) Josephson potential of the DC SQUID for φext = (2e/ )Φext = 0, with the quantum dot junction in the singlet (orange) and doublet (purple) state.The dashed lines represented the two lowest transmon energy levels in each branch of the Josephson potential, with the arrow denoting the resulting transition frequency, which can differ for the two quantum dot junction states (orange and purple arrows for singlet and doublet, respectively).

FIG. 2 .
FIG. 2. Device overview.(a) Diagram of the microwave circuit.A coplanar waveguide transmission line (green center conductor) is capacitively coupled to a grounded LC resonator.The resonator consists of an island (yellow) capacitively and inductively (pink) shunted to ground (blue).The resonator is in turn capacitively coupled to a transmon island (red), which is shunted to ground capacitively as well as via two parallel Josephson junctions.(b) False-colored optical microscope image of device A showing the qubit island, the resonator island, the resonator inductor, the transmission line, the electrostatic gates and ground.(c) False-colored scanning electron micrograph (SEM) of the transmon's Josephson junctions, showing the InAs/Al nanowire into which the junctions are defined.The By component of the magnetic field is used to tune Φext [60].Bz is the magnetic field component parallel to the nanowire.(d) False-colored SEM of the quantum dot junction in which the quantum dot is gate defined.The three bottom gates have a width and spacing of 40 nm, although this is obfuscated by the dielectric layer placed on top [60].
. The doublet sector includes the states of odd parity, which have S = 1/2: |↑ = d † ↑ |0 and |↓ = d † ↓ |0 .It is convenient to introduce the energy ξ = + U/2, corresponding to half of the energy gap in the singlet sector, so that ξ = 0 corresponds to the electron-hole symmetry point, where |0 and |2 are degenerate in energy.The ground state of H dot belongs to the doublet sector for |ξ/U | < 1/2.

FIG. 3 .
FIG. 3. Resonator and transmon spectroscopy.(a) Vp dependence of single-tone spectroscopy for φext = 0, showing the resonator's transition frequency.Vp is a virtual gate voltage defined as a linear combination of Vc and Vr (see text).(b) Zoom-in of (a) in the plunger gate range indicated with dashed lines in (a).(c) Same as (b) but for φext = π.(d) Vp dependence of two-tone spectroscopy for φext = 0, showing the transmon's transition frequency.The black arrow indicates f 0 01 , the transmon frequency set by the reference junction when the quantum dot is pinched off.(e) Same as panel (d) but for φext = π.For panels (a-e) Vt = 182 mV and V l = 470 mV.(f) Theoretical estimates of the singlet (orange), doublet (purple) and reference junction-only (dotted, grey) transmon frequencies as ξ is varied for φext =0.Solid (dashed) lines indicate which quantum dot occupation corresponds to the ground (excited) state.(g) Same as panel (f) but for φext = π.For panels (f-g) ∆/h = 46 GHz, U/∆ = 12.2, ΓL/∆ = 1.05 and ΓR/∆ = 1.12.ft and f01 denote, respectively, the frequency of the second tone in two-tone spectroscopy and the first transmon transition frequency (see text).

FIG. 4 .
FIG. 4. Flux and plunger gate dependence.(a) ∆f01 = f01 − f 0 01 versus Vp and φext as extracted from two-tone spectroscopy.The dashed line is a sinusoidal guide for the eye, denoting the transition boundary in line with the theoretical expectation [60].(b)-(d) Three linecuts of f01 versus φext at representative Vp values, indicated in panel (a) and Fig. 3(be).The dotted line indicates f 0 01 .For all panels Vt = 182 mV.

FIG. 5 .
FIG. 5. Tunnel gate dependence.(a) ∆f01 versus Vp and Vt at φext = 0, where Vt is a virtual gate voltage defined as a linear combination of Vcand Vr(see text).The blue region corresponds to a negative supercurrent contribution from the quantum dot junction, while the red region corresponds to a positive contribution.(b) The same measurement as (a) repeated for φext = π.(c) Linecuts of (a) and (b) at Vp = 395 mV overlayed with best-fits based on NRG calculations.(d) Extracted dependence of ΓL,R on Vt.(e) Calculated transmon frequencies based on NRG calculations at φext = 0 as matched to the measured data, with the Vt axis as given in figure (d).The color bar is shared with panel (a).(f) Same as (e) but for φext = π, with the same color bar as (b).For the NRG calculations in panels (c-f) we fix ∆/h = 46 GHz and U/∆ = 12.2.

FIG. 6 .
FIG. 6. Parallel magnetic field dependence.(a)Borders between singlet and doublet regions for Bz = 10 mT (black) and Bz = 200 mT (blue).φext = 0 (b) Same as (a) but for φext = π.The grey marker denotes the gate point used to re-calibrate the flux axis after varying Bz [60].(c) ∆f01 versus Bz and φext, measured at the gate point indicated in (a) and (b) with a green marker.The sinusoidal dashed line serves as a guide for the eye, in line with the transition boundary expect from theory [60].
FIG. 7. Dependence of parity lifetimes on Vp and φext for device B. (a) A 18 ms cut of a continuously measured time trace integrated in time bins of tint = 11.4 µs, revealing jumps between two distinct states.X is the common axis onto which the quadratures of the outgoing microwave field are rotated to obtain the highest SNR, which takes a value of 3.3 in this panel.(b) 1D histogram of the response in (a) (black) and the best fit of a double Gaussian line-shape (gray).The separation of their centers δx and their width σ together define the SNR.The ratio of their amplitudes determines the ratio of the lifetimes.For panels (a-b) Vp = 551.4mV and φext = 0. (c) Vp dependence of |S21| at φext = 0. (d) Vp dependence of the extracted lifetimes at φext = 0. Markers indicate the mean while error bars indicate the maximum and minimum values of 10 consecutive 2 s time traces.The SNR is shown in greyscale in the background.For points where SNR < 1, the extracted lifetimes are discarded.(e) 2D map of log10(T d /Ts) versus Vp and φext, extracted from a 2 s time trace for each pixel.White pixels indicate points at which SNR < 1, while grey regions indicate where the resonator frequencies of singlet and doublet states overlap and thus cannot be distinguished.
., M.P.V., B.v.H. and A.K. conceived the experiment.Y.L. developed and provided the nanowire materials.A.B., M.P.V., L.S., L.G. and J.J.W prepared the experimental setup and data acquisition tools.L.S. deposited the nanowires.A.B. and M.P.V. designed and fabricated the device, performed the measurements and analysed the data, with continuous feedback from L.S., L.G., J.J.W, C.K.A, A.K. and B.v.H. R.A., J.A., B.v.H. and R.Z. provided theory support during and after the measurements, and formulated the theoretical framework to analyze the experiment.R.Z. performed the NRG calculations.A.B., M.P.V. and B.v.H. wrote the code to compute the circuit energy levels and extract experimental parameters.L.P.K., R.A., A.K. and B.v.H. supervised the work.A.B., M.P.V., R.Z. and B.v.H. wrote the manuscript with feedback from all authors.

Figure S2 .
Figure S2.Numerical matching of model parameters (a) Calculation of the value of ΓR that leads to a singlet-doublet transition with other model parameters held fixed.Here we fix φext = 0 and ξ = 0.A value of zero indicates that no such transition occurs.(b) Calculation Eq. (S5) in the U − ΓL plane evaluated at ∆ = 46 GHz.(c) Same as (b) in the ∆ − ΓL plane for U/∆ = 12.2.(d) Same as (b) in the U − ∆ plane for ΓL/∆ = 1.19.

Figure S3 .
Figure S3.Numerical matching in ΓL − U plane (a) Difference in the calculated and measured singlet qubit frequency at φext = 0 evaluated at ∆ = 46 GHz.(b) Same as (a) for the doublet qubit frequency at φext = 0. (c) Same as (a) for the doublet qubit frequency at φext =π.(d) The absolute sum of the differences in panels (a-c).

Figure S4 .
Figure S4.Numerically calculated transmon frequency maps (a,e,i) Boundaries between singlet and doublet ground states extracted from NRG calculations for φext = 0, φext = π, and ΓR = 1.23ΓL respectively.Panels (b-d), (f-h), (j-l) show how the singlet qubit frequency, the doublet qubit frequency, and the combined result conditioned on the ground state of panels (a,e,i) respectively depend on the parameters.Each row shares the same color map.This leads to saturation of the color map in the panels corresponding to the unconditioned singlet and doublet qubit frequencies, but facilitates comparison to the experimental results.For all panels U/∆ = 12.2 and ∆ = 46 GHz.

Figure S5 .
Figure S5.State population versus ξ and φext (a) Ratio of the expected state population as calculated from Eq. (S6) for a temperature of 400 mK.The colormap is saturated to facilitate comparison to main text Fig. 7. (b) Difference between doublet and singlet energy.Each contour indicates a boundary of equal energy difference.Parameters are the same as that of Fig. S4.

Figure S6 .
Figure S6.Chip design.(a) The chip of device A, containing four nearly identical devices coupled to the same transmission line.The image is taken after wire-bonding onto a PCB.(b) The chip of device B, incorporating an input capacitor in the transmission line (enlarged in inset).The image is taken before wire-bonding onto a PCB.

Figure S8 .Figure S9 .
Figure S8.Flux jumps under |B| = 9 mT for device A. Multiple flux jumps and a distorted periodicity observed at low magnetic fields disappear when |B| > 9 mT.Here, Bz = Bx = 0

2 T
Figure S10.Reference junction characterization for device A. (a) Vj dependence of single-tone spectroscopy when the quantum dot junction is pinched-off (Vc = 52.4mV, V l = 470 mV, Vr = 373 mV).At low Vj values the reference junction is pinched-off and EJ ∼ 0, thus the resonator is at its bare resonance frequency.As Vj increases, the resonator frequency increases non-monotonically due to mesoscopic fluctuations of the overall increasing transmission of different junction channels.(b) Vj-dependence of two-tone spectroscopy for the Vj range indicated in (a) with a dashed line rectangle.The black lines in (a) and (b) indicate the Vj = 640 mV set-point which sets the transmon frequency to its set-point used for the main text figures, f01 = f 0 01 = 4.4 GHz.(c) Line-cut of (a) at the Vj set-point, showing a resonance.(d) Line-cut of (b) at the Vj set-point, showing two peaks.The highest peak, at higher frequency, appears when the second tone frequency matches the transmon frequency (ft = f 0 01 ).The lower peak corresponds to f02/2 and shows the anharmonicity of the transmon.For (d), the first tone frequency fr is fixed at the bottom of the resonance, indicated with a grey arrow in (c).(e) Bz evolution of f 0 01 and f02/2 at Vj = 640 mV.

Figure S11 .
Figure S11.Quantum dot gates characterization for device A. (a) Gate voltage dependence (V l = Vc = Vr = Vgate) of single-tone spectroscopy, showing how the quantum dot junction is pinched off at Vgate values lower than 300 mV.(b-d) Vc, V l and Vr dependence, respectively of single-tone spectroscopy.In each panel, the two unused gates are kept at 1250 mV.This shows how each of the three quantum dot gates can independently pinch off the quantum dot junction.For all panels, the reference junction is closed (Vj = −200 mV).The red line in (c) indicates the fixed value of V l = 470 mV at which all main text figures are taken.
Figure S12.Quantum dot tune up for device A. (a) Single-tone spectroscopy measured at V l = 250 mV, exhibiting two small resonances.Here the reference junction is fully closed (Vj = -200 mV).The red line indicates the readout frequency used in panel (b).(b) Single frequency readout of the resonator.Bright colors indicate a shift in the resonance frequency, marking the onset of supercurrent through the dot.The red line indicates the V l value of panel (a).(c) Same as panel (a) but with the reference junction opened to the Vj = 640 mV setpoint used throughout the manuscript.The two junctions in parallel form a SQUID, increasing the qubit frequency and in turn the resonance frequency.Measured at φext = 0. (d) Same as panel (b) but with the reference junction set to Vj = 640 mV and φext = 0, measured at the frequency indicated with a red line in (c).For (a-d), Vc = 100 mV (close to pinchoff), indicated with a black line in Fig. S11.(e) f01 versus φext at fixed Vj = 640 mV, for three quantum dot gates setpoints corresponding to a quantum dot junction which is fully closed (grey), slightly open (violet) or very open (blue) showing the DC SQUID behaviour of the two parallel Josephson junctions.

Figure S13 .
Figure S13.Quantum dot gate dependence for device A. (a) Single frequency readout of the resonator at the frequency indicated in Fig. S12(c) with a red line, performed versus Vc and Vr for fixed V l = 470 mV.(b) Same as (a) but versus Vc and V l and for fixed Vr = 425 mV.(c) Single-tone spectroscopy versus Vc, measured at V l = 470 mV and Vr = 425 mV, revealing a quantum dot resonance.For all panels φext = 0.

Figure S14 .
Figure S14.Gate compensation for device A. (a) Single-tone spectroscopy versus Vcat Vr = 427 mV.(b) Single frequency readout of the resonator measured versus the central (Vc) and right (Vr) quantum dot gate voltages, performed at a at fixed V l 470 mV.The red line indicates the Vr value of panel (a).(c) Resonator spectroscopy versus Vp at Vt = 180 mV.(d) Same as (b) but in the transformed coordinate frame, measured vs. the virtual plunger (Vp) and right tunnel (Vt) gate voltages.In (a) and (c), the red lines indicate the readout frequencies used in panels (b) and (d), respectively.For all panels φext = 0.

Figure S15 .
Figure S15.Extended Vt dependence.(a) ∆fr versus Vp and Vt at φext = 0, revealing singlet (red) and doublet (blue) ground state regions separated by sharp transitions.(b) Same as (a) but for φext = π.We note that the plunger gate axis is shifted by about 5 mV with respect to (a) and the data shown in the main text, which we speculate is due to an irreproducible gate jump.Dashed rectangles indicate the gate ranges in which the measurements of Fig. 5 of the main text are taken.

Figure S16 .
Figure S16.State selective spectroscopy.(a-d) φext dependence of single-tone spectroscopy at four representative Vp values, indicated in Fig. 3 in the main text.Traces at intermediate Vp values show two resonances simultaneously due to switches on timescales faster than the integration time.(e) Linecut of (c) at φext = 0, indicated by a black line in (c).(f) Two-tone spectroscopy at the same settings as in (e), with the first tone at the frequency of the singlet resonance.The measurement shows a peak at the transmon frequency of the singlet state.(g) Same as (f) but with the readout frequency corresponding to the doublet resonance, which shows a peak at the transmon frequency of the doublet state.

Figure S17 .
Figure S17.Data analysis for magnetic field dependence of device A. (a) Borders between singlet and doublet regions for Bz = 10 mT (black) and Bz = 200 mT (blue).Solid and empty markers correspond to φext = π and φext = 0, respectively.(b) and (c) show ∆f01 versus Bz and φext, measured at the two gate points indicated in (a) with grey and green markers, respectively.In (b), the singlet is the ground state for all Bz.This gate point is used to identify a flux reference for each Bz.For (c), there is a singlet-doublet ground state transition with Bz, where the sinusoidal dashed line serves as a guide for the eye.(d) f01 versus φext for the three Bz values indicated in (b) and (c).The dotted line indicates f 0 01 , which decreases with Bz as shown in Fig. S10(e).

Figure S18 .
Figure S18.Parity lifetime analysis.(a) Logarithmic-scale histogram of the resonator response in the (I,Q)-plane after integrating a 2 s time trace with time bins of tint = 11.4 µs.It exhibits two separate Gaussian distributions whose centers define an axis, X, indicated with a dashed line.(b) Power spectral density (black) of an unintegrated 2 s time trace projected onto the X axis.In grey, best fit of a Lorentzian lineshape with a white noise background (Eq.(S8)).(c) 18 ms cut of the integrated response projected onto the X axis, revealing jumps between two distinct states.(d) 1D histogram of the response in (a) projected onto the X axis (black) and the best fit of a double Gaussian line-shape (grey, Eq. (S7)).For all panels V l = 325 mV, Vt = −60 mV, Vp = 551.4mV, Bz = 0 and φext = 0.
1,2 are the relative populations of singlet and doublet occupation, x 1,2 are the centers of each Gaussian and σ is their standard deviation.For the data shown in Fig.S18, the fit results in A 1 = 2169σ, A 2 = 506σ, x 1 = 0.37σ and x 2 = −6.19σ,from which we determine the SNR = |x 1 − x 2 |/2σ = 3.28.

Figure S19 .
Figure S19.Flux dependence of parity lifetimes (a) φext dependence of single-tone spectroscopy at Vp = 551.4mV.(b) φext dependence of the parity lifetimes extracted following the analysis in Fig. S18 at Vp = 551.4mV.Markers indicate the mean and error bars indicate the maximum and minimum values of 10 consecutive 2 s time traces.SNR = |δx| 2σ is shown in greyscale in the background.For points where SNR < 1, the extracted parity lifetimes are not shown as we do not consider them reliable.Measured at the same Vt , V l and Bz as for Fig. S18.

Figure S20 .
Figure S20.Power and temperature dependence of parity lifetimes across the singlet-doublet transition (a) 2D map of log10(T d /Ts) versus Vp and φext, extracted from a 2 s time trace for each pixel.This is the same panel as Fig. 7(e) in the main text.(b) Vp dependence of single-tone spectroscopy at φext = 0, across a singlet/doublet transition.For (a) and (b), the mixing chamber temperature is 18 mK and the readout power is −22 dBm.(c-f) Readout power dependence at 18 mK of the extracted parity lifetimes at the plunger points indicated in (a) and (b).Markers indicate the mean and error bars indicate the maximum and minimum values of 10 consecutive 2 s time traces.The SNR is shown in greyscale in the background.For points where SNR < 1, the extracted parity lifetimes are discarded.(g-j) Same as (c-d) but versus temperature and at a power of −22 dBm.All powers are given at the fridge input.

Figure S21 .
Figure S21.Tunnel, plunger and power dependence of parity lifetimes (a) ∆fr = fres − f 0 res versus Vp and Vt measured at φext = 0.It shows a regions of constructive and destructive interference, separated by sharp dome-like boundary.(b) Two-dimensional map of log10(T d /Ts) versus Vp and Vt, measured at a power of −22 dBm.(c) Same as (b) but for a power of −14 dBm.(d-e) Power dependence of the extracted parity lifetimes at the gate points indicated in (a-c).All powers are given at the fridge input.