Phase-engineering the Andreev band structure of a three-terminal Josephson junction

TL;DR

This study demonstrates the control and measurement of a synthetic Andreev band structure in a three-terminal Josephson junction, revealing hybridized Andreev states and their phase-dependent properties.

Contribution

It introduces a method to engineer and observe the Andreev band structure in a three-terminal junction with independent phase control, advancing understanding of Andreev matter.

Findings

Observation of phase-controlled Andreev molecule formation

Detection of avoided crossings indicating hybridization

Identification of band structure anisotropies in phase space

Abstract

In hybrid Josephson junctions with three or more superconducting terminals coupled to a semiconducting region, Andreev bound states may form unconventional energy band structures, or Andreev matter, which are engineered by controlling superconducting phase differences. Here we report tunnelling spectroscopy measurements of three-terminal Josephson junctions realised in an InAs/Al heterostructure. The three terminals are connected to form two loops, enabling independent control over two phase differences and access to a synthetic Andreev band structure in the two-dimensional phase space. Our results demonstrate a phase-controlled Andreev molecule, originating from two discrete Andreev levels that spatially overlap and hybridise. Signatures of hybridisation are observed in the form of avoided crossings in the spectrum and band structure anisotropies in the phase space, all explained by a numerical model. Future extensions of this work could focus on addressing spin-resolved energy levels, ground state fermion parity transitions and Weyl bands in multiterminal geometries.

Equations

$$E_{\pm }=\frac{E_{\mathrm{L}}\left({\varphi }_{\mathrm{L}}\right)+E_{\mathrm{R}}\left({\varphi }_{\mathrm{R}}\right)}{2}\pm \sqrt{{\left(\frac{\gamma }{2}\right)}^2+{\left[\frac{E_{\mathrm{L}}\left({\varphi }_{\mathrm{L}}\right)-E_{\mathrm{R}}\left({\varphi }_{\mathrm{R}}\right)}{2}\right]}^2},$$

$$E_{\alpha }=-\Delta \sqrt{1-{\tau }_{\alpha }^2{\sin }^2\frac{{\varphi }_{\alpha }}{2}},\alpha \in \left\{\mathrm{L},\mathrm{R}\right\}$$

$$E_{\alpha }\approx -\Delta \sqrt{1-{\tau }_{\alpha }^2\left[1-\frac{{\left({\varphi }_{\alpha }-\pi \right)}^2}{4}\right]}\approx -\frac{\Delta }{2}\sqrt{{\left({\varphi }_{\alpha }-\pi \right)}^2+{\left(2r_{\alpha }\right)}^2},$$

$$E_{\pm }\approx -\frac{\Delta }{2}\sqrt{{\delta }^2+\frac{{\left({\phi }_{\mathrm{L}}+{\phi }_{\mathrm{R}}\right)}^2}{4}+\frac{{\left({\phi }_{\mathrm{L}}-{\phi }_{\mathrm{R}}\right)}^2}{4}\pm \left({\phi }_{\mathrm{L}}-{\phi }_{\mathrm{R}}\right)\sqrt{{\delta }^2+\frac{{\left({\phi }_{\mathrm{L}}+{\phi }_{\mathrm{R}}\right)}^2}{4}}},$$

$$\begin{gathered} E_{\pm }\approx -\frac{\Delta }{2\sqrt{2}}\times \\ \times \sqrt{2{\delta }^2+{\phi }_{\mathrm{L}}^2+{\left(2r_{\mathrm{L}}\right)}^2+{\phi }_{\mathrm{R}}^2+{\left(2r_{\mathrm{R}}\right)}^2\pm \sqrt{{\left[{\phi }_{\mathrm{L}}^2+{\left(2r_{\mathrm{L}}\right)}^2-{\phi }_{\mathrm{R}}^2-{\left(2r_{\mathrm{R}}\right)}^2\right]}^2+4{\delta }^2\left[{\left({\phi }_{\mathrm{L}}-{\phi }_{\mathrm{R}}\right)}^2+{\left(2r_{\mathrm{L}}\right)}^2+{\left(2r_{\mathrm{R}}\right)}^2\right]}}. \end{gathered}$$

$$H=H_{\text{DD}}+\sum _{\alpha =\mathrm{M},\mathrm{L},\mathrm{R}}{H}_{\text{S},\alpha }+H_{\text{DS}}.$$

$$H_{\text{DD}}={ϵ}_1\sum _{\sigma }{d}_{1\sigma }^{†}{d}_{1\sigma }+{ϵ}_2\sum _{\sigma }{d}_{2\sigma }^{†}{d}_{2\sigma }+t\sum _{\sigma }{d}_{1\sigma }^{†}{d}_{2\sigma }+\text{h.c.},$$

$$H_{\text{S},\alpha }=\sum _{\mathbf{k}\sigma }{\xi }_{\mathbf{k}}{c}_{\mathbf{k}\sigma ,\alpha }^{†}{c}_{\mathbf{k}\sigma ,\alpha }+\sum _{\mathbf{k}}\left(\Delta {\mathrm{e}}^{i{\varphi }_{\alpha }}{c}_{\mathbf{k}\uparrow ,\alpha }^{†}{c}_{-\mathbf{k}\downarrow ,\alpha }^{†}+\Delta {\mathrm{e}}^{-i{\varphi }_{\alpha }}{c}_{-\mathbf{k}\downarrow ,\alpha }{c}_{\mathbf{k}\uparrow ,\alpha }\right),$$

$$H_{\text{DS}}=v_{\mathrm{M},1}\sum _{\mathbf{k},\sigma }{d}_{1\sigma }^{†}{c}_{\mathbf{k}\sigma ,\mathrm{M}}+v_{\mathrm{M},2}\sum _{\mathbf{k},\sigma }{d}_{2\sigma }^{†}{c}_{\mathbf{k}\sigma ,\mathrm{M}}+v_{\mathrm{L},1}\sum _{\mathbf{k},\sigma }{d}_{1\sigma }^{†}{c}_{\mathbf{k}\sigma ,\mathrm{L}}+v_{\mathrm{R},2}\sum _{\mathbf{k},\sigma }{d}_{2\sigma }^{†}{c}_{\mathbf{k}\sigma ,\mathrm{R}}+\text{h.c.},$$

$${\mathbf{G}}_{\text{DD}}={\mathbf{g}}_{\text{DD}}+{\mathbf{g}}_{\text{DD}}{\mathbf{V}}_{\text{DS}}{\mathbf{g}}_{\text{SS}}{\mathbf{V}}_{\text{DS}}^{†}{\mathbf{G}}_{\text{DD}},$$

$$T_1=\frac{4{\Gamma }_{\mathrm{M},1}{\Gamma }_{\mathrm{L},1}}{{\left({\Gamma }_{\mathrm{M},1}+{\Gamma }_{\mathrm{L},1}\right)}^2},T_2=\frac{4{\Gamma }_{\mathrm{M},2}{\Gamma }_{\mathrm{R},2}}{{\left({\Gamma }_{\mathrm{M},2}+{\Gamma }_{\mathrm{R},2}\right)}^2},$$

$$\left(\begin{array}{c}{\varphi }_{\mathrm{L}}^{\ast }\\ {\varphi }_{\mathrm{R}}^{\ast }\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\cdot \left(\begin{array}{c}{\varphi }_{\mathrm{L}}\\ {\varphi }_{\mathrm{R}}\end{array}\right),$$

$$\left(\begin{array}{c}{\Phi }_{\mathrm{L}}\\ {\Phi }_{\mathrm{R}}\end{array}\right)=\mathbf{M}\cdot \left(\begin{array}{c}{I}_{\mathrm{L}}\\ I_{\mathrm{R}}\end{array}\right)=\left(\begin{array}{cc}{M}_{\mathrm{LL}}& M_{\mathrm{LR}}\\ M_{\mathrm{RL}}& M_{\mathrm{RR}}\end{array}\right)\cdot \left(\begin{array}{c}{I}_{\mathrm{L}}\\ I_{\mathrm{R}}\end{array}\right),$$

$$\mathbf{M}=\left(\begin{array}{cc}6.96\mathrm{pH}& -1.43\mathrm{pH}\\ -1.69\mathrm{pH}& 5.79\mathrm{pH}\end{array}\right).$$

$$\left(\begin{array}{c}\Delta {\Phi }_{\mathrm{L}}\\ 0\end{array}\right)=\mathbf{M}\cdot {\mathbf{\Delta I}}_{\mathrm{L}}^{\ast },\left(\begin{array}{c}0\\ \Delta {\Phi }_{\mathrm{R}}\end{array}\right)=\mathbf{M}\cdot {\mathbf{\Delta I}}_{\mathrm{R}}^{\ast }.$$

$$L_{\mathrm{k}}=\frac{l}{w}\frac{h}{2{\pi }^2}\frac{R_{\square }}{\Delta }\approx 44\mathrm{pH},$$

Full text

Introduction

In a normal conductor interfacing two or more superconductors, charge carriers at energies within the superconducting gap are confined by Andreev reflection processes occurring at the interfaces [1]. As a result, resonant sub-gap electronic excitations known as Andreev bound states (ABSs) arise in the normal region, enabling transport of a Josephson supercurrent between the superconducting terminals [2, 3]. These discrete ABS levels were proposed as a basis for quantum computing applications [4, 5, 6, 7]. More recently, ABSs have been the subject of intense experimental investigation in various material platforms [8, 9, 10, 11, 12, 13, 14, 15], culminating in the coherent control of Andreev pair [16, 17] and spin [18, 19] qubits. While these studies focused on Josephson junctions (JJs) with two superconducting terminals, where the ABS energies depend on a single superconducting phase difference, multiterminal JJs (MTJJs) have also emerged as a promising alternative. In the presence of $N\geq 3$ terminals, the ABS band structure spanned by the $N-1$ independent phase differences was predicted to exhibit a plethora of phenomena, including lifting of the spin degeneracy [20], ground state fermion parity transitions [20], Weyl singularities [21, 22, 23, 24, 25] and other topological properties [26, 27, 28, 29, 30, 31]. Moreover, MTJJs were proposed to realise Andreev molecules [32, 33, 34, 35]—a system where single ABSs overlap and form hybridised energy levels—and explored as a platform to generate Cooper quartets [36, 37, 38, 39, 40]. Extensive experimental work was also conducted on MTJJs in the presence of DC current bias [41, 42, 43, 44] and microwave irradiation [45]. While tunnelling spectroscopy of metallic three-terminal JJs (3TJJs) was performed previously [46], establishing phase control over the superconducting proximity gap, further investigation is needed to understand the properties of ABS bands in multiterminal devices. The perspective of engineering synthetic band structures by exploiting the higher dimensionality of the phase space is particularly attractive, as it could enable effects unattainable in two-terminal geometries.

Here, we report on an experimental realisation of superconductor–semiconductor 3TJJs and study ABSs in the system with tunnelling spectroscopy. Owing to the independent control over two superconducting phase differences, we probe the Andreev band structure in the two-dimensional (2D) phase space and find signatures of hybridisation between highly transmissive ABSs, resulting in the formation of an Andreev molecule. Our measurements are supported by a theoretical model and demonstrate the feasibility of Andreev matter and phase-engineering of Andreev bands in hybrid nanostructures.

Results

Implementation of phase-controllable 3TJJ

The device under investigation, shown in Fig. 1a–c, was realised in an InAs/Al heterostructure [47, 48]. Selective etching of the Al layer defined three superconducting terminals (labelled L, M and R) coupled to a normal region, constituting a hybrid 3TJJ. The three terminals were connected to a common node (D) forming two closed loops; while leads L and M were directly connected to D via Al strips, a superconductor–normal–superconductor (SNS) JJ was integrated on terminal R. The junction, with a length of 40 nm and a width of 5 $\mu$m, was designed to have a critical current much larger than that existing between any pairs of L, M and R, hence the superconducting phase difference across the junction is neglected for the following discussion. A fourth superconducting lead (S) was employed as a probe to enable DC tunnelling spectroscopy of the sub-gap states in the 3TJJ. Metallic gate electrodes and flux-bias lines were patterned on top of an insulating layer uniformly deposited across the entire sample. The transmission between the probe and the 3TJJ was tuned via two gates energised by the common voltage $V_{\mathrm{T}}\equiv V_{\mathrm{TL}}=V_{\mathrm{TR}}$, which was set to $-1.07~{}\mathrm{V}$ to enter the tunnelling regime: in this configuration, the probe was weakly coupled to the 3TJJ and its influence on the rest of the circuit was limited. Gate tuning of the SNS junction enabled its operation as a switch, with the ON (OFF) state defined by $V_{\mathrm{Switch}}=0$ ($V_{\mathrm{Switch}}=-1.5~{}\mathrm{V}$). This allowed for the connection or disconnection of terminal R from D, hence electrostatically selecting between a three-terminal (switch ON) and a two-terminal configuration (switch OFF). Two additional gate voltages were kept to $V_{\mathrm{Probe}}=150~{}\mathrm{mV}$ and $V_{\mathrm{G}}=50~{}\mathrm{mV}$. A current $I_{\mathrm{L(R)}}$ injected into the left (right) flux-bias line generated an external magnetic flux $\mathit{\Phi}_{\mathrm{L(R)}}$ threading the left (right) superconducting loop, thus tuning the phase difference between L (R) and M, and enabling control over the whole 2D phase space. As schematically depicted in Fig. 1b, a DC voltage bias $V_{\mathrm{SD}}$ was applied to the probe S and the differential tunnelling conductance $G$ was measured between S and D with standard lock-in techniques. Experiments were performed in a dilution refrigerator with base temperature below $10~{}\mathrm{mK}$. Further details about materials, fabrication and measurement are provided in the Methods section.

Probing the Andreev band structure in the 2D phase space

In Fig. 1d,e, the voltage bias was set to $-170~{}\mu\mathrm{V}$ and the tunnelling conductance was measured as a function of the currents $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ injected into the flux-bias lines, resulting in a scan over an extended region of the 2D phase space at constant energy. Resonances in conductance correspond to peaks in the density of states (DOS) of the normal region under study [8, 15] and represent ABSs in the 3TJJ. Each state is intersected twice per period at $V_{\mathrm{SD}}=-170~{}\mu\mathrm{V}$ (i.e., at an energy below its maximum), giving rise to the characteristic appearance of pairs of lines in these maps. When the switch is set to the ON state (Fig. 1d), we observe periodic features as a function of both $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$, attributed to the presence of two distinct ABSs whose energies disperse with $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$, respectively. The cross-dependence between the flux-bias lines accounts for the finite slope of the $(\mathit{\Phi}_{\mathrm{L}},\mathit{\Phi}_{\mathrm{R}})$ axes, as indicated by the black arrows. Remarkably, resonances associated to different ABSs are connected to each other in proximity of the intersections, forming closed diamond-like loops and avoided crossings at the corners of the diamonds. Each state undergoes a phase shift when intersecting the other, defining a zig-zag trend. Next, we set the switch junction voltage to $-1.5~{}\mathrm{V}$ (OFF, Fig. 1e) and observe that the ABS resonances depending on $\mathit{\Phi}_{\mathrm{R}}$ disappear, while the complex 2D periodic pattern is transformed into a simple structure depending on a single flux. Notably, the slope of the $\mathit{\Phi}_{\mathrm{L}}$-dependent resonances differs between panels d and e (see dotted yellow lines), which is compatible with the periodic phase shift present only when the switch is ON.

Andreev dispersion along phase space linecuts

With the phase space overview acquired at constant voltage bias, we select cut lines $\gamma_{\mathrm{1-6}}$ (arrows in Fig. 1d,e) along which the tunnelling conductance is measured as a function of $V_{\mathrm{SD}}$. The full datasets along $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are displayed in Fig. 2a–c. Each shows a transport gap $2\mathit{\Delta}/e\approx 310~{}\mu\mathrm{V}$ ($e$ is the elementary charge), consistent with a superconducting gap of the Al probe $\mathit{\Delta}\approx 155~{}\mu\mathrm{eV}$, and an electron–hole-symmetric spectrum revealing phase-dependent ABSs. The presence of regions with finite conductance at $e|V_{\mathrm{SD}}|\lesssim\mathit{\Delta}$ is ascribed to broadened features in the DOS of the superconducting probe for energies $|E|\sim\mathit{\Delta}$. This could be due to a combination of quasiparticle-lifetime broadening [49] and additional subgap bound states forming between probe and 3TJJ. On either side of the spectrum, we notice two differential conductance peaks at $V_{\mathrm{SD}}=\pm 155~{}\mu\mathrm{V}$ and $\pm 175~{}\mu\mathrm{V}$, whose position in bias does not vary appreciably with $\gamma_{i}$. The first is attributed to the multiple Andreev reflection peak at $\pm\mathit{\Delta}/e$, while the second, specific to this device, might be related to mesoscopic defects in the tunnelling probe or to a region in the device with a larger superconducting gap. These peaks provide a contribution to the measured differential conductance, adding to that of ABS resonances and spectroscopy background. This accounts for the intensity modulation of these peaks depending on $\gamma_{i}$. Further supported by measurements at different tunings of the probe (see Supplementary Note 4), we do not observe a distortion of the ABS dispersion when the states cross the peaks. In Fig. 2d–i, all six linecuts are plotted in restricted $V_{\mathrm{SD}}$ and $\gamma_{i}$ ranges for better clarity. In each case, ABSs approach the transport gap edge very closely, which indicates near-unity transmission. We also note that such highly transmissive ABSs intersect $V_{\mathrm{SD}}=-170~{}\mu\mathrm{V}$ twice per period, thus accounting for the pairs of resonances in Fig. 1d,e. Interestingly, when comparing the $\gamma_{1}$ and $\gamma_{2}$ linecuts (Fig. 2d,e), we find anisotropic ABS phase dispersion in the vicinity of $(\pi,\pi)$ phase ($\gamma_{\mathrm{1,2}}=0$), with a narrow, cusp-like shape versus a broader and flatter peak, respectively. In the spectroscopy measurements performed along $\gamma_{\mathrm{3-5}}$, i.e., lines parallel to $\gamma_{2}$ but offset from an intersection point of Fig. 1d, two distinct highly transmissive ABSs appear. When their separation in phase is small, the states partially mix and an avoided crossing is observed (Fig. 2f), an effect which is weaker at larger separation (Fig. 2g), until it is completely suppressed (Fig. 2h). Finally, linecut $\gamma_{6}$, where the switch junction is kept in the OFF state, reveals a single highly transmissive ABS (Fig. 2i).

Andreev molecule model and numerical simulations

To simulate ABSs arising in a 3TJJ, we study a minimal theoretical model comprising three superconducting terminals with phases $\phi_{\mathrm{L}}$, $\phi_{\mathrm{R}}$ and $\phi_{\mathrm{M}}$, as schematically sketched in the inset of Fig. 3a. Due to gauge invariance, only two phases are independent, therefore we set $\phi_{\mathrm{M}}\equiv 0$. Between lead L (R) and M, we assume a highly transmissive channel, fully described by two coupling parameters to the leads and hosting a spin-degenerate ABS, whose energy $E_{\mathrm{L(R)}}$ depends on $\phi_{\mathrm{L(R)}}-\phi_{\mathrm{M}}\equiv\phi_{\mathrm{L(R)}}$. The choice of the coupling parameters led to transmissions $T_{1}\approx 0.998$ and $T_{2}\approx 0.992$ for the left and right ABS, respectively. Coupling between the ABSs is described by the parameter $t$, resulting in a minimal set of five numerical parameters, plus a broadening parameter (more details of the model are explained in Supplementary Note 1). When the two ABSs are isolated ($t=0$), their energy dispersions solely depend on the transmission of the respective channels [2]. However, for finite coupling $t\neq 0$, hybridisation between the states is enabled in the regions of $(\phi_{\mathrm{L}},\phi_{\mathrm{R}})$ where the isolated energies $E_{\mathrm{L}}$ and $E_{\mathrm{R}}$ are comparable, resulting in an Andreev molecule [32, 33]. To study the experimentally relevant case, we assume $t=1.1\mathit{\Delta}$, indicating ABSs strongly coupled to each other. In Fig. 3a, we display a numerical simulation of the Andreev energy bands as a function of $\phi_{\mathrm{L}}$ and $\phi_{\mathrm{R}}$ for negative energies (i.e., below the Fermi level $E_{\mathrm{F}}\equiv 0$), noting that specular bands exist at positive energies owing to electron–hole symmetry. The finite coupling between the ABSs results in hybridised bands with a pronounced splitting. In the region with both phases tuned near to $\pi$, the band closer to zero energy shows a striking anisotropy in the phase space and the dispersion is strongly modified compared to that of high-transmission ABSs in a ballistic two-terminal JJ. To better compare experimental and numerical results, we calculate the DOS at fixed energy $E=-0.09\mathit{\Delta}$ as a function of the cross-coupled phases $\phi_{\mathrm{L}}^{*}$ and $\phi_{\mathrm{R}}^{*}$ (Fig. 3b), and as a function of energy along the three phase space linecuts $\gamma_{1-3}$ (Fig. 3c–e). All simulations qualitatively reproduce the key features observed in the measurements of Figs. 1(d) and 2, and lead to the same order of magnitude for the avoided crossings. In the constant-energy simulation (Fig. 3b), we confirm the presence of a periodic pattern characterised by avoided crossings and phase shifts of the ABS resonances near the intersection points, where individual ABSs are connected to each other forming closed loops. Further, the spectra presented in Fig. 3c–e resemble the results shown in Fig. 2d–f: the ABS dispersion approaching zero energy has a sharp peak along $\gamma_{1}$ and is broader along $\gamma_{2}$. The individual ABSs reappear with a significant separation when probed along $\gamma_{3}$, while maintaining a sizeable avoided crossing. The experimentally measured phase shifts are larger than those calculated theoretically. The enhanced shift is likely produced by the finite inductive coupling between the loops, which is not included in the numerical model and implies a coupling between the fluxes $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$. The change of slope of the ABS resonances between Figs. 1d and 1e is consistent with the phase shift and is related to the same coupling mechanism [50, 51]. These effects are discussed in more detail in Supplementary Note 6.

Tomography of the Andreev band structure

Another visualisation of the ABS energy bands is provided by combining multiple constant-energy cut planes—each showing the dependence on both superconducting phase differences—to achieve a tomographic representation of the band structure. For this purpose, we measured the tunnelling conductance as a function of $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ for several values of $V_{\mathrm{SD}}$. The outcome is summarised in Fig. 4, where panels a–c display three theoretically calculated planes in the low-energy spectrum and d–l report experiments for nine values of $V_{\mathrm{SD}}$. At $V_{\mathrm{SD}}\approx-\mathit{\Delta}/e$ (i.e., near zero energy), where the ABSs are probed in proximity of their maximum, we observe single resonances (Fig. 4a,d), which split into pairs at more negative bias (Fig. 4b,e). We confirm the presence of the features described for Fig. 1d and Fig. 3b. Notably, in Fig. 4c,f resonances related to different states are connected to each other by arcs enclosing an oval region. At even more negative bias, additional resonances arise, compatible with low-transmission ABSs appearing in the spectrum from $V_{\mathrm{SD}}\approx-200~{}\mu\mathrm{V}$ (see Fig. 2). Similar to the high-transmission ABSs discussed previously, low-transmission states first occur as single lines (Fig. 4g) and then split into pairs (Fig. 4h,i). Several additional modes emerge at larger $|V_{\mathrm{SD}}|$, making it difficult to resolve individual states while approaching the continuum (Fig. 4j–l). We remark that the theoretical model assumes only two high-transmission modes and does not include any additional states with lower transmission, unlike the experimentally measured spectrum. Therefore, a direct comparison of the spectrum at high $|E|$ is beyond the scope of the model.

Our main experimental findings, including avoided crossings and phase shifts in the constant-bias planes of Figs. 1d and 4, as well as the anisotropy highlighted in Figs. 2d,e and 3c,d, were qualitatively reproduced on a second device (see Supplementary Note 3). These measurements suggest the generality of the phenomena observed.

Discussion

Supported by our theoretical model of coupled ABSs, we interpret the experimental results summarised in Figs. 1, 2 and 4 as evidence of coupling and hybridisation between two highly transmissive ABSs in the normal region of the 3TJJ. Thus, our devices constitute an implementation of an Andreev molecule, comparable to Refs. [32, 33]. Avoided crossings in both the phase space and the energy spectrum, together with the anisotropic ABS dispersion motivate our interpretation. Unlike previous experimental studies of Andreev molecules in two-terminal geometries [52, 53], this realisation is in an open system and is a direct manifestation of the phase-controlled, multidimensional Andreev band structure.

Our system can be considered the magnetic dual of a double quantum dot [54], where electric fields controlled by gate voltages, capacitance and charge on the dots (quantised in units of $e$) are substituted by magnetic fields controlled by currents in flux-bias lines, inductance and magnetic flux threading superconducting loops (quantised in units of the superconducting flux quantum $\mathit{\Phi}_{0}=h/2e$, with $h$ the Planck constant). Both in a double quantum dot and in the present Andreev molecule, overlapping wave functions of two discrete and localised states, coupling in a middle region, result in avoided crossings between their otherwise degenerate energy levels. In our device structure, two discrete levels, namely high-transmission ABSs, form in the short L–M and R–M junctions (whose minimum length is lithographically $30~{}\mathrm{nm}$) and are coupled to each other by their close proximity.

Spin-resolved ABSs are not observed in the experiments, well described by a spin degenerate model, despite the presence of spin–orbit coupling in our system. We note that here the spin–orbit length, $l_{\mathrm{SO}}\sim 150~{}\mathrm{nm}$ for InAs [55], is larger than the separation between pairs of terminals of the 3TJJ, resulting into a relatively weak strength of spin–orbit coupling. Enlarging the size of the 3TJJ would thus be required to resolve spin–orbit splitting of ABSs.

In conclusion, ABSs in hybrid 3TJJs were investigated with tunnelling spectroscopy measurements. Owing to the individual control over two superconducting phase differences, we explored a synthetic Andreev band structure and found signatures of coupling and hybridisation between two highly transmissive ABSs, consistent with their overlap in the 3TJJ region and the formation of an Andreev molecule. In the 2D phase space probed at constant voltage bias, we observed periodic patterns with avoided crossings and phase shifts near the intersections between ABS resonances. We measured the spectrum along selected linecuts of the phase space, finding a strong anisotropy of the ABS band structure and avoided crossings between the states. The experiments are well described by a theoretical model of two coupled ABSs. Our results provide new insights into the physics of multiterminal devices, establish phase control over the ABS band structure and demonstrate the feasibility of realising exotic Andreev matter. Future studies of multidimensional band structures could focus on phase-engineering spin-resolved Andreev levels [20], ground state fermion parity transitions [20, 56, 57] and topological bands, including Weyl singularities [21, 22, 23, 24, 25].

Note. We recently became aware of the unpublished data of Refs. [58] and [59], where three-terminal devices were investigated.

Methods

Materials and fabrication

Devices were fabricated in a III–V heterostructure grown with molecular beam epitaxy techniques on an InP (001) substrate. The stack consisted of a step-graded InAlAs buffer layer covered by an $\mathrm{In_{0.75}Ga_{0.25}As}$/InAs/$\mathrm{In_{0.75}Ga_{0.25}As}$ quantum well and two monolayers of GaAs. The InAs layer, hosting a two-dimensional electron gas (2DEG), was 8 nm thick and buried 13 nm below the surface. On top of the III–V stack, a 15 nm thick Al layer was deposited in situ without breaking vacuum. Characterisation of the 2DEG in a gated Hall bar revealed a peak mobility of $18000~{}\mathrm{cm^{2}V^{-1}s^{-1}}$ at an electron sheet density of ${8\cdot 10^{11}~{}\mathrm{cm^{-2}}}$. This resulted in an electron mean free path $l_{e}\gtrsim 260~{}\mathrm{nm}$, indicating that both the three-terminal Josephson junction and the two-terminal switch junction were in the ballistic regime.

First, large mesa structures were isolated, suppressing parallel conduction between devices and across the middle regions of the superconducting loops. This was done by selectively etching the Al layer with Transene type D, followed by a second chemical etch to a depth of $\sim 380~{}\mathrm{nm}$ into the III–V material stack, using a $220:55:3:3$ solution of $\mathrm{H_{2}O:C_{6}H_{8}O_{7}:H_{3}PO_{4}:H_{2}O_{2}}$. Next, Al was defined by wet etching with Transene type D at $50\mathrm{{}^{\circ}C}$ for $4~{}\mathrm{s}$. The dielectric, deposited on the entire chip by atomic layer deposition, consisted of a $3~{}\mathrm{nm}$ thick layer of $\mathrm{Al_{2}O_{3}}$ and a $15~{}\mathrm{nm}$ thick layer of $\mathrm{HfO_{2}}$. Gate electrodes and flux-bias lines were defined by evaporation and lift-off. In a first step, $5~{}\mathrm{nm}$ of Ti and $20~{}\mathrm{nm}$ of Au were deposited to realise the fine features of the gates; in a second, a stack of Ti/Al/Ti/Au with thicknesses $5~{}\mathrm{nm}$, $340~{}\mathrm{nm}$, $5~{}\mathrm{nm}$ and $100~{}\mathrm{nm}$ was deposited to connect the mesa structure to the bonding pads and to define the flux-bias lines.

Measurement techniques

Experiments were performed in a dilution refrigerator with base temperature at the mixing chamber below $10~{}\mathrm{mK}$. The sample was mounted on a QDevil QBoard sample holder system, without employing any light-tight enclosure. Electrical contacts to the devices, excepts for the flux-bias lines, were provided via a resistive loom with QDevil RF and RC low-pass filters at the mixing chamber stage, and RC low-pass filters integrated on the QBoard sample holder. Currents in the flux-bias lines were injected via a superconducting loom with only QDevil RF filters at the mixing chamber stage. Signals were applied to all gates and flux-bias lines via home-made RC filters at room temperature. Electrical transport measurements were performed with low-frequency AC lock-in techniques. A fixed AC voltage $\delta V_{\mathrm{SD}}=5~{}\mu\mathrm{V}$ at frequency $211~{}\mathrm{Hz}$ and a variable DC voltage $V_{\mathrm{SD}}$ were applied to a contact at the superconducting probe (labelled S in Fig. 1a). The AC current $\delta I$ and the DC current $I_{\mathrm{SD}}$ flowing in the grounded terminal D were measured via a current-to-voltage (I–V) converter. By measuring the AC voltage $\delta V$ between terminals S and D in a four-terminal configuration, the differential conductance $G\equiv\delta I/\delta V$ was determined. The refrigerator was equipped with a vector magnet which, despite not being utilised for the experiments, produced a small magnetic field offset. Hence, arbitrary offsets in the flux-bias line currents $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ of $-18~{}\mu\mathrm{A}$ and $74~{}\mu\mathrm{A}$ were considered in datasets, in such a manner that the point where $I_{\mathrm{L}}=I_{\mathrm{R}}=0$ was at the centre of a diamond-like region in the constant-bias maps.

Data availability

The data presented in this study have been deposited in Zenodo [https://zenodo.org/record/8360770]. Further data that support the findings of this study are available upon request from the corresponding author.

Code availability

Computer code used to perform the numerical simulations presented in this work has been deposited in Zenodo [https://zenodo.org/record/8360770].

Supplementary Information: Phase-engineering the Andreev band structure of a three-terminal Josephson junction

Supplementary Note 1: Theory

In the following, we discuss hybridisation of two Andreev bound states (ABSs) between three superconducting leads. The problem has been investigated extensively in the context of the Andreev molecule [32, 33, 34, 35]. Here, we treat hybridisation empirically with a coupling parameter $\gamma$ introducing an avoided crossing in the spectrum. Assuming relatively weak coupling between the two ABSs, the spectrum resembles the energy levels of the independent states except for the points where they would cross, where hybridisation leads to the anticrossing. A simple perturbation theory can be applied if the energies $E_{\mathrm{L,R}}$ are sufficiently larger than $\gamma$:

[TABLE]

where

[TABLE]

is the energy of an independent ABS between two superconducting leads with phase difference $\phi_{\alpha}$ and energy gap $\mathit{\Delta}$ in a channel with transmission $\tau_{\alpha}^{2}$. We recall that the perturbation theory is applicable for $\left|E_{\alpha}(\phi_{\alpha})\right|\gg|\gamma|$: this condition does not hold if the transmissions of the junctions are large ($1-\tau_{\alpha}\ll 1$) and the phase differences are close to $(2n+1)\pi$, which is the regime studied in the experiment. Nevertheless, in such experimentally relevant limit we follow the approach of Refs. [33, 34] and expand the independent ABS energies around $\phi_{\alpha}=\pi$:

[TABLE]

with $r_{\alpha}^{2}=1-\tau_{\alpha}^{2}\ll 1$ the reflection amplitude squared. In the limit of perfect transmission $r_{\mathrm{L}}=r_{\mathrm{R}}=0$, the equations derived in Refs. [33, 34] for the ballistic regime can be employed:

[TABLE]

where $\varphi_{\alpha}=\phi_{\alpha}-\pi$ and $\delta=\gamma/\mathit{\Delta}$ is the dimensionless coupling parameter. The formula is applicable for any points close to $(\phi_{\mathrm{L}},\phi_{\mathrm{R}})=\left[\pi(2n+1),\,\pi(2m+1)\right]$, with $n$ and $m$ integers. The corresponding constant-energy cut plane of the spectrum at $E=-0.01\mathit{\Delta}$ is shown in Supplementary Fig. 1a as a function of the phase differences. Here, hybridisation of purely ballistic channels introduces avoided crossings along the anti-diagonal direction $\phi_{\mathrm{L}}=2\pi-\phi_{\mathrm{R}}$.

Moreover, we consider small but finite reflection amplitudes $r_{\alpha}\ll 1$ for each channel and use the effective Hamiltonian introduced in [33] (by projecting on the low-energy states) to derive the spectrum:

[TABLE]

The constant-energy plane at $E=-0.01\mathit{\Delta}$ is displayed in Supplementary Fig. 1b, where we observe that non-zero reflection amplitudes introduce additional avoided crossings along the diagonal direction $\phi_{\mathrm{L}}=\phi_{\mathrm{R}}$.

In order to go beyond the limiting case accessible analytically, valid for weak coupling between the ABSs, we perform numerical simulations of a comparable system. Andreev states are modelled as two single-level quantum dots (QDs) coupled to three superconducting terminals. A coupling between the QDs accounts for hybridisation of the ABSs in the three-terminal junction. The model is schematically shown in Supplementary Fig. 2. For each ABS, the only relevant parameter is the transmission of the corresponding junction [60], hence we can arbitrarily opt for any microscopic model. The QD model offers a flexible tool to describe a scattering region and is well suited for numerical analyses of the system, even for strong coupling between the ABSs, where the analytical approach cannot be used. The coupling between each QD and the superconducting terminals can be arbitrarily strong, hence this model is also suited to describe ABSs in the open regime, rather than isolated dots weakly coupled to the leads [61].

The total Hamiltonian of the system is expressed as:

[TABLE]

Here, $H_{\text{DD}}$ is the Hamiltonian of the double QD:

[TABLE]

where the creation operator $d_{j\sigma}^{\dagger}$ of the $j$-th QD, its energy $\epsilon_{j}$ and the coupling parameter $t$ are defined. The superconducting terminals are described by Bardeen–Cooper–Schriffer Hamiltonians [62]:

[TABLE]

where $c_{\bm{k}\sigma,\alpha}^{\dagger}$ is the creation operator of an electron with momentum $\bm{k}$ and spin $\sigma$ in lead $\alpha\in\left\{\mathrm{M,L,R}\right\}$, $\mathit{\Delta}$ is the superconducting gap (assumed to be the same for all leads), $\phi_{\alpha}$ the superconducting phase of lead $\alpha$ and $\xi_{\bm{k}}$ the normal state dispersion in the leads. The tunnel coupling between each QD and the superconductors is expressed by the term $H_{\text{DS}}$, having the form:

[TABLE]

where $v_{\alpha,j}$ denotes the coupling of the $j$-th QD to lead $\alpha$.

Two superconductor–QD–superconductor junctions are identified: between leads L and M through QD1 and between leads R and M through QD2. In the absence of coupling between the QDs ($t=0$), the ABS energy of QD1(2) depends only on $\phi_{\mathrm{L(R)}}-\phi_{\mathrm{M}}$. For finite coupling ($t\neq 0$) the ABSs hybridise, which is the case presented in the Main Text. We do not consider a direct ballistic channel between terminals L and R due to their larger separation in the device compared to the L–M and R–M junctions. Nevertheless, the energy of a hypothetical ABS in such channel would depend on $\phi_{\mathrm{L}}-\phi_{\mathrm{R}}$, hence it would be close to the superconducting gap in proximity of the points where $\phi_{\mathrm{L}}=\phi_{\mathrm{R}}=\pi$. As a consequence, any hybridisation of the additional state with the two included in the model would be suppressed near those regions of phase space. To evaluate this model in the limit of strong coupling between superconducting leads and QDs, i.e., $\mathit{\Gamma}_{\alpha,j}\equiv\pi N_{0}v_{\alpha,j}^{2}>\mathit{\Delta}$ (where $N_{0}$ is the normal density of states in the leads), we calculate the ABS density of states by determining the Green’s function of the coupled system with the Dyson equation:

[TABLE]

where $\bm{G}_{\text{DD}}$ is the dressed Green’s function of the double QD, $\bm{g}_{\text{DD}}$ the unperturbed Green’s function of the two QDs, $\bm{V}_{\text{DS}}$ the coupling between the QDs and the superconducting leads and $\bm{g}_{\text{SS}}=\text{diag}\left(\bm{g}_{\text{SS,M}},\,\bm{g}_{\text{SS,L}},\,\bm{g}_{\text{SS,R}}\right)$ the unperturbed Green’s function of the three leads. From this expression, the ABS density of states is calculated as $\rho=-\frac{1}{\pi}\text{Im}\left\{\text{tr}\left(\bm{G}_{\text{DD}}\right)\right\}$ and the eigenenergies are obtained from the poles of $\bm{G}_{\text{DD}}$, yielding the results shown in Figs. 3 and 4 of the Main Text. As previously discussed, the only relevant parameters for the hybridised ABSs are the inter-QD coupling $t$ and the transmission amplitudes, which are related to the QD–lead couplings $\mathit{\Gamma}_{\alpha,j}$ via the expressions:

[TABLE]

valid in the limit of strong lead–QD coupling. These parameters influence the region of the avoided crossing between the two ABS, as illustrated in Supplementary Fig. 3. In all simulations, a broadening parameter $\eta$ was assumed in the Green’s functions. We observe qualitative agreement between the numerical simulations and the analytical study previously discussed (see Supplementary Fig. 1).

The parameters used for the simulations in the Main Text were: $t=1.1\mathit{\Delta}$, $\mathit{\Gamma}_{\mathrm{L},1}=5.5\mathit{\Delta}$, $\mathit{\Gamma}_{\mathrm{M},1}=\mathit{\Gamma}_{\mathrm{R},2}=6\mathit{\Delta}$, $\mathit{\Gamma}_{\mathrm{M},2}=5\mathit{\Delta}$, corresponding to transmissions $T_{1}\approx 0.998$ and $T_{2}\approx 0.992$, and broadening ${\eta=0.02\mathit{\Delta}}$. Moreover, we recall that the constant-energy planes of Figs. 3b and 4a–c were plotted by introducing a cross-dependence between the phase differences $\phi_{\mathrm{L}}-\phi_{\mathrm{M}}\equiv\phi_{\mathrm{L}}$ and $\phi_{\mathrm{R}}-\phi_{\mathrm{M}}\equiv\phi_{\mathrm{R}}$ to better represent the experimental data, where a cross-dependence between the two flux-bias lines was present (see also discussion in Section Supplementary Note 5: Mutual inductance matrix). In particular, the cross-coupled phases $\phi_{\mathrm{L}}^{*}$ and $\phi_{\mathrm{R}}^{*}$ were defined as linear combinations of $\phi_{\mathrm{L}}$ and $\phi_{\mathrm{R}}$ with the transformation:

[TABLE]

where the coefficients $a=0.9708$, $b=0.2400$, $c=0.2832$ and $d=1.122$ were used.

Supplementary Note 2: Gate dependence of differential conductance

The transmission between the superconducting probe and the three-terminal Josephson junction (3TJJ) was controlled with two gate electrodes, denoted tunnel gates and set to the same voltage $V_{\mathrm{T}}\equiv V_{\mathrm{TL}}=V_{\mathrm{TR}}$, and had a finite dependence on the gates energised by the voltages $V_{\mathrm{Probe}}$ and $V_{\mathrm{G}}$, denoted probe gate and global gate respectively (see Fig. 1a,c). We always kept $V_{\mathrm{Probe}}=150~{}\mathrm{mV}$. The effect of the gate at voltage $V_{\mathrm{Switch}}$ (switch gate) on the probe transmission was negligible. The differential conductance $G$ measured as a function of $V_{\mathrm{T}}$ and $V_{\mathrm{SD}}$, with $V_{\mathrm{SD}}$ the DC voltage bias applied to the probe, is shown in Supplementary Fig. 4a,b for $V_{\mathrm{G}}=50~{}\mathrm{mV}$ and $V_{\mathrm{G}}=-150~{}\mathrm{mV}$ respectively. In both cases the conductance, hence the probe transmission, decreases with $V_{\mathrm{T}}$, transitioning from open to tunnelling regime. In the former, we observe a zero-bias conductance peak, corresponding to a remnant of supercurrent between the probe and the superconducting leads, and several peaks at finite bias, related to multiple Andreev reflection (MAR) processes. In the tunnelling regime, where the conductance at large $|V_{\mathrm{SD}}|$ is substantially lower than the conductance quantum $G_{0}=2e^{2}/h$, a transport gap of $\approx 310~{}\mu\mathrm{V}$ is present in the spectrum, consistent with a superconducting gap of Al $\mathit{\Delta}\approx 155~{}\mu\mathrm{eV}$. Pronounced features appearing at the edges of the gap correspond to Andreev bound states (ABSs). As shown in Supplementary Fig. 4b, by further lowering $V_{\mathrm{T}}$ the probe could be completely pinched off. We remark that tunnelling spectroscopy in our devices was performed in a superconductor–insulator–superconductor (SIS) configuration, hence the $G(V_{\mathrm{SD}})$ traces result from a convolution product between two relatively complex densities of states [8, 15] and do not provide a direct measurement of the gap hardness. Nevertheless, since $G\approx 0$ across a significant voltage range around $V_{\mathrm{SD}}=0$ (for any flux-line currents $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$, as seen in Fig. 2 of the Main Text), we consider the gap hardness to be comparable to previous reports [63], where also a density-of-states broadening to what we noted in the Main Text was observed.

The dependence on the global gate is presented in Supplementary Fig. 4c for $V_{\mathrm{T}}=-1.07~{}\mathrm{V}$. We find that its main effect on the spectrum is also to change the transmission of the probe, which is progressively reduced for decreasing $V_{\mathrm{G}}$ until the pinch off is reached. Notably, in both the $V_{\mathrm{T}}$ and the $V_{\mathrm{G}}$ dependence, no sharp conductance peaks relatable to resonant transport via spurious quantum dots are observed.

Tunnelling spectroscopy measurements were performed on a second device, fabricated on the same chip of the first and measured in the same cool down. Device 2 was lithographically similar to Device 1, except for the width of the superconducting probe and the shape of the tunnel and probe gates (see Supplementary Fig. 4d). The dependences on the tunnel gate voltage $V_{\mathrm{T}}$ and on the global gate voltage $V_{\mathrm{G}}$ are plotted in Supplementary Fig. 4e and f respectively, showing features qualitatively very similar to Device 1.

Supplementary Note 3: Results for Device 2

The main experimental results shown in the Main Text for Device 1 were qualitatively reproduced in Device 2, as illustrated in Supplementary Figs. 5 and 6. Here, the tunnel gates were set to $V_{\mathrm{T}}=-1.395~{}\mathrm{V}$ and $V_{\mathrm{T}}=-1.42~{}\mathrm{V}$ respectively, with the global gate voltage set to $V_{\mathrm{G}}=-150~{}\mathrm{mV}$ and the probe gate voltage to $V_{\mathrm{Probe}}=100~{}\mathrm{mV}$. In these configurations, the probe was in the tunnelling regime and its transmission was comparable to that of Device 1 for the measurements presented in the Main Text.

In Supplementary Fig. 5, we show the constant-bias planes (i.e., $G$ as a function of $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ at fixed values of $V_{\mathrm{SD}}$) corresponding to Figs. 1d,e, and 4d–l of the Main Text. When the switch junction is in the ON state (defined by $V_{\mathrm{Switch}}=0$), we confirm the presence of avoided crossings between a $\mathit{\Phi}_{\mathrm{L}}$-dependent ABS and a $\mathit{\Phi}_{\mathrm{R}}$-dependent ABS, with the resonances associated to one state connecting to those of the other, and phase shifts occurring near the avoided crossings. The 2D pattern is strongly simplified when the switch is OFF ($V_{\mathrm{Switch}}=-1.5~{}\mathrm{V}$, Supplementary Fig. 5b), as the $\mathit{\Phi}_{\mathrm{R}}$-dependent ABS disappears and no sign of hybridization is observed. The band structure tomography is displayed in Supplementary Figs. 5c–k and is compatible with the results presented in the Main Text (see Fig. 4).

Further, we select linecuts of the phase space, indicated by the coloured arrows in Supplementary Figs. 5a,b, along which we perform bias-dependent spectroscopy (Supplementary Fig. 6, to be compared with Fig. 2). Again, we observe strong dispersion anisotropy when comparing linecut $\gamma_{1}$ to $\gamma_{2}$. Incidentally, while we still note a conductance peak at $V_{\mathrm{SD}}=\pm 155~{}\mu\mathrm{V}$, whose position in bias does not vary appreciably with $\gamma_{i}$ and which is attributed to MAR processes, this device does not reveal a second peak at larger $|V_{\mathrm{SD}}|$. This shows that the peak at $V_{\mathrm{SD}}=\pm 175~{}\mu\mathrm{V}$ in Device 1 is a device-specific feature. Since the ABS dispersion of Device 1 is qualitatively reproduced in Device 2, we corroborate that the peak at $\pm 175~{}\mu\mathrm{V}$ in the former does not interact with the ABSs or affect their main properties.

Supplementary Note 4: Results for different $V_{\mathrm{G}}$

The measurements of Device 1 shown in the Main Text were obtained in a single gate configuration (except for the switch gate voltage, set to $-1.5~{}\mathrm{V}$ for Figs. 1e, 2i and to $0~{}\mathrm{V}$ elsewhere). The main results were reproduced in different gate configurations, defined by setting a new value of $V_{\mathrm{G}}$ and adjusting $V_{\mathrm{T}}$ to maintain a comparable probe transmission, crucially remaining in the tunnelling regime. In Supplementary Figs. 7–12 we show the three cases $V_{\mathrm{G}}=-66~{}\mathrm{mV}$, $V_{\mathrm{G}}=-150~{}\mathrm{mV}$ and $V_{\mathrm{G}}=-250~{}\mathrm{mV}$, as described in the following. Remarkably, all the key features discussed in the Main Text were qualitatively reproduced in these configurations, highlighting the generality of our results and further supporting our conclusions. Incidentally, we did not observe a variation in the transmission of the highly transmissive ABSs as a function of the global gate voltage, namely they approached the edges of the transport gap at $|V_{\mathrm{SD}}|\approx 155~{}\mu\mathrm{V}$ very closely for any $V_{\mathrm{G}}$. This was verified in the range $-300~{}\mathrm{mV}\leq V_{\mathrm{G}}\leq 150~{}\mathrm{mV}$, noting that at $V_{\mathrm{G}}=-300~{}\mathrm{mV}$ the sub-gap states were at a limit of visibility but still preserved the same dispersion with phase.

• •

$V_{\mathrm{G}}=-66~{}\mathrm{mV}$: Supplementary Figs. 7a,b, 7c–k and 8 correspond to Figs. 1d,e, 4 and 2a–h respectively. $V_{\mathrm{T}}$ is adjusted to $-1.038~{}\mathrm{V}$. Interestingly, we observe that the two ABSs have very similar visibility in this regime (namely, the resonances associated with each state have similar conductance), in contrast to the regime shown in the Main Text ($V_{\mathrm{G}}=50~{}\mathrm{mV}$), where the ABS dispersing with $\mathit{\Phi}_{\mathrm{L}}$ had higher conductance. As a possible explanation, we hypothesise that, when $V_{\mathrm{G}}$ is reduced, the electrostatic potential profile at the tunnelling barrier shifts slightly, in such a manner that its maximum displaces from left to right. At the optimised value $V_{\mathrm{G}}=-66~{}\mathrm{mV}$, transport between the probe and either side of the three-terminal region is symmetric, resulting in similar conductance when tunnelling into either ABS. Furthermore, in the constant-bias planes of Supplementary Fig. 7, and particularly in panel b, we notice an additional set of resonances with the same dependence on $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ as the left ABS. This is attributed to a very small phase modulation of the MAR peak, i.e., the nearly flat line at $V_{\mathrm{SD}}=-155~{}\mu\mathrm{V}$ that is visible in Supplementary Fig. 8. These peak, particularly prominent in the present regime, accounts for a conductance background in the constant-bias maps and its modulation translates into a phase-periodic background.

• •

$V_{\mathrm{G}}=-150~{}\mathrm{mV}$: Supplementary Figs. 9a,b, 9c–k and 10 correspond to Figs. 1d,e, 4 and 2 respectively. $V_{\mathrm{T}}$ is adjusted to $-1.02~{}\mathrm{V}$, where the probe transmission is smaller than in the previous regimes: the differential conductance at large $V_{\mathrm{SD}}$ (near the range limits) is reduced by approximately a factor $2$. This results in a lower signal-to-noise ratio, due to which periodic noise features become visible in the constant-bias measurements (Supplementary Fig. 9), as long as $G$ is relatively small. At the same time, the conductance resonances at $V_{\mathrm{SD}}=\pm 155~{}\mu\mathrm{V}$ and $V_{\mathrm{SD}}=\pm 175~{}\mu\mathrm{V}$ in the spectra of Supplementary Fig. 10 become substantially less prominent, and so is their modulation as a function of $\gamma_{i}$. Visibly, the ABS dispersion is not distorted when the states cross the horizontal peaks, remaining qualitatively very similar to that observed in Fig. 2. This provides further support to the absence of an interaction between the peaks at $V_{\mathrm{SD}}=\pm 155~{}\mu\mathrm{V},\pm 175~{}\mu\mathrm{V}$ and the ABSs. Finally, in this gate configuration, the ABS depending on $\mathit{\Phi}_{\mathrm{R}}$ shows the highest conductance, consistent with the argument presented for the case $V_{\mathrm{G}}=-66~{}\mathrm{mV}$.

• •

$V_{\mathrm{G}}=-250~{}\mathrm{mV}$: Supplementary Figs. 11a,b, 11c–k and 12 correspond to Figs. 1d,e, 4 and 2a,b,d,e,i respectively. $V_{\mathrm{T}}$ is adjusted to $-984~{}\mathrm{mV}$, where the probe transmission is comparable to the regime $V_{\mathrm{G}}=-150~{}\mathrm{mV}$.

Supplementary Note 5: Mutual inductance matrix

As described in the Main Text, our devices feature two flux bias lines where currents $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ are injected. Since $I_{\mathrm{L}}$ generates a magnetic field that is stronger over the left superconducting loop than over the right one, it controls mainly the external magnetic flux $\mathit{\Phi}_{\mathrm{L}}$ threading the left loop, associated to the superconducting phase difference between the terminals L and M (see Fig. 1a). Similarly, $I_{\mathrm{R}}$ tunes mostly $\mathit{\Phi}_{\mathrm{R}}$, thus the phase difference between R and M, in such a manner that the combination of the two flux lines enables full phase control over a two-dimensional space. Nevertheless, each flux line has also a finite coupling to the opposite loop, thus $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$ depend on both $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ according to the linear relation:

[TABLE]

where $\mathbf{M}$ is the mutual inductance matrix. We calculate $\mathbf{M}$ from the constant-bias conductance measurement of Fig. 1d, plotted again in Supplementary Fig. 13a. First, we observe that the origin $(I_{\mathrm{L}},I_{\mathrm{R}})=(0,0)$ corresponds to the origin of the flux space $(\mathit{\Phi}_{\mathrm{L}},\mathit{\Phi}_{\mathrm{R}})=(0,0)$. We associate the centres of the adjacent diamond-like regions to the addition or subtraction of one superconducting flux quantum $\mathit{\Phi}_{0}=h/(2e)$ to either $\mathit{\Phi}_{\mathrm{L}}$ or $\mathit{\Phi}_{\mathrm{R}}$. By evaluating Eq. 13 in two centre points, for example those related to $(\mathit{\Phi}_{0},0)$ and $(0,\mathit{\Phi}_{0})$, we write a $4\times 4$ equation system and determine:

[TABLE]

The ratio $M_{\mathrm{RR}}/M_{\mathrm{LL}}\approx M_{\mathrm{LR}}/M_{\mathrm{RL}}\approx 0.8$ approximately corresponds to the ratio between the length of the vertical segments of the right and left flux line. Knowing $\mathbf{M}$, we can apply the linear transformation of Eq. 13 to convert the $(I_{\mathrm{L}},I_{\mathrm{R}})$ axes to the $(\mathit{\Phi}_{\mathrm{L}},\mathit{\Phi}_{\mathrm{R}})$ axes, as illustrated in Supplementary Fig. 13b for the dataset of panel a.

Supplementary Note 6: Phase shifts and inductances in the system

In the constant-bias measurements showing the dependence on $I_{\mathrm{L}}$ and $I_{\mathrm{R}}$ with the switch junction ON, such as in Figs. 1d and 4d–i, we remarked the presence of phase shifts occurring when resonances associated to different states (i.e., the $\mathit{\Phi}_{\mathrm{L}}$- and the $\mathit{\Phi}_{\mathrm{R}}$-dispersing ABS) intersect each other. Concomitantly, we noted a slope variation of the $\mathit{\Phi}_{\mathrm{L}}$-dispersing ABS depending on the ON/OFF state of the switch junction (yellow dashed line in Fig. 1e). In this section, we quantify the shifts and relate them to the different inductive contributions in the device. We find that a Josephson-like mutual inductive coupling between the two loops accounts for the presence of both the shifts and the slope difference.

Since a shift is present for both ABSs, it corresponds to a vector $(\Delta\mathit{\Phi}_{\mathrm{L}},~{}\Delta\mathit{\Phi}_{\mathrm{R}})$ in the space of the two external magnetic fluxes $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$. For the extraction of this vector, we consider again the dataset of Fig. 4e from the Main Text, shown also in Supplementary Fig. 13c and, upon applying the basis transformation described in Section Supplementary Note 5: Mutual inductance matrix, in Supplementary Fig. 13d as a function of $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$. As a guide for the eye, we overlay lines following the dips between pairs of ABS resonances. In panel c, we mark the (vector) shifts $\mathbf{\Delta I}^{*}_{\mathrm{L}}$ and $\mathbf{\Delta I}^{*}_{\mathrm{R}}$ along the periodicity directions, while in panel d the shifts $\Delta\mathit{\Phi}_{\mathrm{L}}$ and $\Delta\mathit{\Phi}_{\mathrm{R}}$ are by construction parallel to the axes. These quantities are linked by the following relations:

[TABLE]

From the data, we find $\Delta\mathit{\Phi}_{\mathrm{L}}\approx\Delta\mathit{\Phi}_{\mathrm{R}}\approx 0.092\cdot\mathit{\Phi}_{0}\equiv\Delta\mathit{\Phi}$, showing that the effect is symmetric for the two ABSs.

This phase shift is then expressed as $\Delta\mathit{\Phi}=M_{\mathrm{cpl}}\cdot I_{\mathrm{circ}}$, namely as the product between a mutual inductance, accounting for the coupling between two fluxes, and a circulating current in either loop. Neither $M_{\mathrm{cpl}}$ nor $I_{\mathrm{circ}}$ can be directly determined from the experimental data (only their product $\Delta\mathit{\Phi}$), therefore we estimate the inductive coupling term considering different origins. The geometric mutual inductance between the two loops of our devices is $M_{\mathrm{geom}}\approx 7.3~{}\mathrm{pH}$ [64]. Since the loops share a strip of epitaxial Al, its inductance, dominated by the kinetic inductance $L_{\mathrm{k}}$, has to be added to $M_{\mathrm{geom}}$ [65, 66]. The kinetic inductance of the strip is estimated as [67]:

[TABLE]

where $l=25~{}\mu\mathrm{m}$ and $w=1~{}\mu\mathrm{m}$ are the length and width of the strip, $R_{\square}\approx 1.5~{}\Omega$ the normal state resistivity of the heterostructure stack measured in a Hall bar geometry (where the Al film was not removed) and $\mathit{\Delta}\approx 180~{}\mu\mathrm{eV}$ the superconducting gap of Al. The geometric and kinetic contribution lead to a combined mutual inductance $M_{\mathrm{loops}}\approx 51.3~{}\mathrm{pH}$. If we assume that $M_{\mathrm{cpl}}=M_{\mathrm{loops}}$, we require a circulating current $I_{\mathrm{circ}}\approx 3.7~{}\mu\mathrm{A}$ for the shift $\Delta\mathit{\Phi}=0.092\cdot\mathit{\Phi}_{0}$. However, this value is much larger than any reasonable estimate of the critical current between each pair of terminals, which constitutes an upper bound to the supercurrent circulating in the loops. As a consequence, $M_{\mathrm{cpl}}$ must be substantially larger than $M_{\mathrm{loops}}$ to limit $I_{\mathrm{circ}}$.

The only contribution that we have not discussed so far is a Josephson-like coupling due to the 3TJJ. Based on its geometry, we expect our 3TJJ to have a critical current of the order of $100~{}\mathrm{nA}$ and a Josephson inductance of a few nH, which would explain the phase shift to a good degree. We deduce that, in our devices, the phase shifts derive from a Josephson-like coupling term rather than the geometric and kinetic inductive couplings, in analogy with existing couplers between superconducting qubits containing Josephson elements [50, 51].

This mutual coupling between $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$ is also responsible for the slope difference depending on the state of the switch junction (see Fig. 1d,e). In the OFF state, no circulating current can flow in the right loop, thus no mutual effect is present. In the $(I_{\mathrm{L}},I_{\mathrm{R}})$ plane (Fig. 1e of the Main Text), the ABS resonance is parallel to the $\mathit{\Phi}_{\mathrm{R}}$ axis (indicated in Fig. 1d). When the switch junction is ON, the slope of the ABSs deviates from the directions of $\mathit{\Phi}_{\mathrm{L}}$ and $\mathit{\Phi}_{\mathrm{R}}$. This is clearly visible on the $(\mathit{\Phi}_{\mathrm{L}},\mathit{\Phi}_{\mathrm{R}})$ axes (Supplementary Fig. 13b,d), as the resonances are not vertical or horizontal. In fact, due to the mutual coupling, the variation of a flux (for example $\mathit{\Phi}_{\mathrm{L}}$) from [math] to $\mathit{\Phi}_{0}/2$ causes a gradual change in the other ($\mathit{\Phi}_{\mathrm{R}}$), up to a maximum shift. Upon crossing $\mathit{\Phi}_{\mathrm{L}}=\mathit{\Phi}_{0}/2$, the circulating current in the left loop reverses direction, consequently the induced current in the right loop and the $\mathit{\Phi}_{\mathrm{R}}$-shift also flip sign. Therefore, the finite ABS slopes on the $(\mathit{\Phi}_{\mathrm{L}},\mathit{\Phi}_{\mathrm{R}})$ axes (corresponding to the slope difference depending on the state of the switch junction) and the phase shifts are consistent. Both result from the finite Josephson-like mutual inductive coupling between the two loops.

Supplementary Figures