# Observation of Josephson harmonics in tunnel junctions

Dennis Willsch, Dennis Rieger, Patrick Winkel, Madita Willsch, Christian Dickel, Jonas Krause, Yoichi Ando, Raphaël Lescanne, Zaki Leghtas, Nicholas Bronn, et al.

► To cite this version:

Dennis Willsch, Dennis Rieger, Patrick Winkel, Madita Willsch, Christian Dickel, et al. Observation of Josephson harmonics in tunnel junctions. Nature Physics, 2024, 20 (5), pp.815-821. 10.1038/s41567-024-02400-8. hal-04891568

HAL Id: hal-04891568

https://hal.science/hal-04891568v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Observation of Josephson harmonics in tunnel junctions

Received: 25 October 2023

Accepted: 14 January 2024

Published online: 14 February 2024

<figure><img src="image_3.png" /></figure>

Check for updates

Dennis Willsch Ⓔ¹,¹⁷, Dennis Rieger Ⓔ²,³,¹⁷, Patrick Winkel Ⓔ²,³,⁴,⁵,
Madita Willsch Ⓔ¹,⁶, Christian Dickel Ⓔ⁷, Jonas Krause Ⓔ⁷, Yoichi Ando Ⓔ⁷,
Raphaël Lescanne⁸,⁹, Zaki Leghtas⁸, Nicholas T. Bronn Ⓔ¹⁰, Pratiti Deb¹⁰,
Olivia Lanes¹⁰, Zlatko K. Minev Ⓔ¹⁰, Benedikt Dennig²,³, Simon Geisert²,
Simon Günzler Ⓔ², Sören Ihssen², Patrick Paluch Ⓔ²,³, Thomas Reisinger Ⓔ²,
Roudy Hanna Ⓔ¹¹,¹², Jin Hee Bae¹¹, Peter Schüffelgen Ⓔ¹¹,
Detlev Grützmacher Ⓔ¹¹,¹², Luiza Buimaga-larinca Ⓔ¹³, Cristian Morari¹³,
Wolfgang Wernsdorfer Ⓔ²,³, David P. DiVincenzo Ⓔ¹²,¹⁴, Kristel Michielsen Ⓔ¹,⁶,¹²,
Gianluigi Catelani Ⓔ¹⁵,¹⁶ & Ioan M. Pop Ⓔ²,³ ☑

Approaches to developing large-scale superconducting quantum processors must cope with the numerous microscopic degrees of freedom that are ubiquitous in solid-state devices. State-of-the-art superconducting qubits employ aluminium oxide (AlOₓ) tunnel Josephson junctions as the sources of nonlinearity necessary to perform quantum operations. Analyses of these junctions typically assume an idealized, purely sinusoidal current-phase relation. However, this relation is expected to hold only in the limit of vanishingly low-transparency channels in the AlOₓ barrier. Here we show that the standard current-phase relation fails to accurately describe the energy spectra of transmon artificial atoms across various samples and laboratories. Instead, a mesoscopic model of tunnelling through an inhomogeneous AlOₓ barrier predicts percent-level contributions from higher Josephson harmonics. By including these in the transmon Hamiltonian, we obtain orders of magnitude better agreement between the computed and measured energy spectra. The presence and impact of Josephson harmonics has important implications for developing AlOₓ-based quantum technologies including quantum computers and parametric amplifiers. As an example, we show that engineered Josephson harmonics can reduce the charge dispersion and associated errors in transmon qubits by an order of magnitude while preserving their anharmonicity.

The Josephson effect¹,² is the keystone of quantum information processing with superconducting hardware: it constitutes a unique source of low-loss nonlinearity, which is essential for the implementation of superconducting quantum bits, and it plays a similarly fundamental role as the nonlinear current-voltage relation of diodes in semiconductor circuitry. In particular, tunnel Josephson junctions (JJs), formed by

two overlapping superconducting films separated by a thin insulating barrier, have enabled superconducting hardware to become one of the leading platforms for the realization of fault-tolerant quantum computers³⁻⁶. JJs are also at the heart of quantum limited amplification⁷, metrological applications⁸ such as the definition of the voltage⁹ and a possible future current standard¹⁰, and they enable quantum detectors
<figure><img src="image_4.png" /><figcaption>Fig. 1 | Josephson harmonics are relevant for the CφR of tunnel junctions.</figcaption></figure>

The nonlinear CφR is the fingerprint of a JJ, which relates the supercurrent $I(\varphi)$ to the phase $\varphi$ (inset). For tunnel JJs, the CφR has been considered to be purely sinusoidal (dashed grey line; equation (1)), with the maximum given by the critical current $I_c$. However, as we show in this work, even in tunnel JJs, the underlying microscopic complexity of the charge transport can manifest in the contribution of higher harmonics to the CφR. As an example, the red line shows a CφR consistent with measured data (CDI of the KIT sample), which includes the harmonics expected from a mesoscopic model assuming an inhomogeneous AlO$_x$ barrier. The shaded red area shows the difference from the purely sinusoidal CφR. We provide CφRs for all other measured samples in Supplementary Fig. 7.

such as the microwave photon counter¹¹. With the advancement¹²⁻¹⁴ of superconducting artificial atom technology, the measurement and understanding of subtle features in the Josephson effect, similar to the fine structure discovered in natural atoms, is increasingly relevant in setting the accuracy of both circuit control and circuit models.

Although the mesoscopic dimensions of JJs imply the existence of many conduction channels, for tunnel junctions this complexity is usually condensed into a single effective parameter, the critical current $I_c$, in the well-known Josephson current-phase relation, CφR (grey line in Fig. 1):

$$I(\varphi) = I_c \sin \varphi, \qquad (1)$$

where $\varphi$ is the superconducting phase difference across the junction. This simplification is remarkable given the fact that other types of junctions, such as weak links, point contacts and ferromagnetic JJs, generally exhibit non-sinusoidal CφRs containing higher Josephson harmonics: $\sin(2\varphi)$, $\sin(3\varphi)$ and so on¹⁵⁻²¹. Here we show that Josephson harmonics are also relevant for tunnel JJs (Fig. 1).

To understand the limits of the approximation equation (1) for tunnel junctions, we have to take a closer look at commonly used Al-AlO$_x$-Al JJs, fabricated by shadow evaporation²² and schematized in Fig. 2a-c, which reveals a complex microscopic reality. The CφR of the junction is obtained by summing the supercurrents of N conduction channels, $I(\varphi) = \sum_{n=1}^{N} I_n(\varphi)$. Each channel (Fig. 2b) has a transparency-dependent CφR (refs. 16,23) that can be expressed as a Fourier series:

$$I_n(\varphi) \propto \frac{T_n \sin \varphi}{\sqrt{1 - T_n \sin^2(\varphi/2)}} = \sum_{m=1}^{\infty} c_m(T_n) \sin(m\varphi). \qquad (2)$$

The conduction channel transparency $T_n$ is defined as the tunnel probability for an electron impinging on the insulating barrier of channel n, and $c_m(T_n)$ are the order m Fourier coefficients for $I_n(\varphi)$. These coefficients alternate in sign and decay in magnitude with increasing order m (Fig. 2d). The ratio $|c_{m+1}/c_m|$ of successive coefficients increases with $T_n$ (Supplementary Section IA): the more transparent a channel, the more relevant the contribution of higher harmonics. To put it simply, in higher-transparency channels, it is more likely for Cooper pairs to

tunnel together in groups of m, which correspond to the $\sin(m\varphi)$ terms in the CφR.

In the limit $T_n \to 0$, only the $\sin \varphi$ term of equation (2) survives. If all channels in a JJ are in this limit, we recover the purely sinusoidal CφR of equation (1), with the critical current of the junction $I_c$ proportional to the sum of transparencies. Assuming a perfectly homogeneous barrier, for a typical junction with $-\mu m^2$ area and resistance comparable to the resistance quantum, one expects $N \approx 10^6$ and $T_n \approx 10^{-6}$ (refs. 24,25), leading to negligible (below $10^{-6}$) corrections to the purely sinusoidal CφR.

But is this the reality? Here we argue that in the presence of contaminants, atomic scale defects²⁶ and random crystalline orientations of the grains in contact, evidenced by scanning transmission electron microscope (STEM) images and molecular dynamics simulations (Fig. 2c and Supplementary Section IV), we have reasons to doubt it. In fact, about two decades ago, AlO$_x$ barrier inhomogeneity motivated the transition in magnetic junctions to more uniform oxides such as MgO (refs. 27–29). Consequently, we expect a distribution of transparencies in AlO$_x$ (refs. 30,31) with possibly a few relatively high-transparency channels³²,³³ introducing measurable corrections to the CφR (Fig. 1). The microscopic structure of each barrier is therefore imprinted on the CφR of the JJ, and the challenge is how to experimentally access this information.

For our study of tunnel JJs, we use transmon devices³⁴, in which a JJ is only shunted by a large capacitor to form a nonlinear oscillator with the potential energy defined by the CφR of the junction (Fig. 2e). The resulting individually addressable transition frequencies in the microwave regime can be measured using circuit quantum electrodynamics techniques³⁵. We compare the spectra of multiple samples to the prediction of the standard transmon Hamiltonian based on a sinusoidal CφR (equation (1)) and find increasing deviations for the higher energy levels of all samples, as sketched in Fig. 2e,f. Only by accounting for higher harmonics in the CφR are we able to accurately describe the entire energy spectrum. A similar methodology was used in ref. 18 to reconstruct the CφR of a semiconductor nanowire Josephson element. While our study focuses on transmon qubits, the conclusions we draw regarding the CφR of tunnel junctions should trigger a re-evaluation of the current models for tunnel-JJ-based devices used in quantum technology and metrology³⁵⁻³⁹.

Since transmons are widely available in the community, we are able to measure and model the spectra of multiple samples from laboratories around the globe: fixed-frequency transmons fabricated and measured at the Karlsruhe Institute of Technology (KIT; Supplementary Fig. 18) in three cooldowns (CDs; Supplementary Fig. 19) and Ecole Normale Supérieure (ENS) Paris (same device as in ref. 40), a tunable transmon subject to an in-plane magnetic field at the University of Cologne (Köln; identical setup and similar device as in ref. 41; Supplementary Fig. 23) and 20 qubits from the IBM Hanoi processor (IBM). All transmons are based on standard Al-AlO$_x$-Al tunnel junctions (Fig. 2) and are measured in either a three-dimensional architecture or a two-dimensional coplanar waveguide geometry (for detailed descriptions of each sample, see Supplementary Section III). The spectroscopy data consists of (1) transition frequencies $f_{0j}$ into transmon states $j=1, 2, ...$ up to $j=6$, each measured as j-photon transitions at frequencies $f_{0j}/j$, and (2) the resonator frequencies $f_{res,j}$ depending on the transmon state $j=0, 1$ (Methods).

In Fig. 3, we compare the measured transition frequencies to predictions $f_{0j}^{\text{model}}$, obtained by exact diagonalization of two different model Hamiltonians. The first model is the standard transmon model, which has served the community for over 15 years³⁴

$$H_{\text{std}} = 4E_C(n - n_g)^2 - E_J \cos \varphi + H_{\text{res}}, \qquad (3)$$

where $E_C$ is the charging energy, $E_J$ is the Josephson energy, $n_g$ is the offset charge and the operators n and $\varphi$ represent the charge
<figure><img src="image_5.png" /><figcaption>**Fig. 2 | Josephson harmonics result from junction barrier inhomogeneity.**
a, False-coloured scanning electron microscope image of a typical Al–AlOₓ–Al JJ fabricated at KIT. The bottom and top electrodes are coloured blue and yellow, respectively. Inset, circuit symbol for a JJ with phase difference φ across the barrier. b, Cross-section schematic of the superconductor–insulator–superconductor JJ at the location indicated by the dash-dotted line in a. The supercurrent Iₙ(φ) of each conduction channel n = 1, ..., N depends on its transparency Tₙ (equation (2)). We sketch a distribution of multiple low and a few high transparencies T₁, ..., Tₙ in green and red, respectively. c, False-coloured high-angle annular dark field STEM image centred on the AlOₓ tunnel barrier of a typical JJ fabricated at KIT, with average thickness d ≈ 2 nm as indicated by the white arrow. Individual columns of atoms of the Al grain in the top electrode are visible due to zone axis alignment, which is not the case for the bottom Al electrode (additional STEM images with thickness variations and structural</figcaption></figure>

defects such as grain boundaries are shown in Supplementary Fig. 27).

d, Normalized Fourier coefficients cₘ(Tₙ) of the JJ CφR (equation (2)) for a low (10⁻⁶, green) and high (10⁻², red) transparency channel. Note the alternating sign for even and odd order m and the fact that high-transparency channel coefficients (in red) remain relevant to higher order. e, Sketch of how the higher-order terms in the JJ Hamiltonian modulate the potential and shift the energy levels (red) of superconducting artificial atoms compared to a purely cosφ potential (grey). In this Article, we focus on transmon devices, which consist of a large capacitor in parallel to the JJ (refer to the circuit schematic inset). The discrepancy between the models generally increases at higher levels. f, The higher-order Josephson harmonics also influence the charge dispersion of the transmon levels versus offset charge nₗ. The two branches per energy level correspond to a change between even and odd charge parity (that is, quasiparticle tunnelling²⁹,⁸⁰; Supplementary Fig. 23 in Supplementary Section IIIC).

normalized by twice the electron charge and the phase difference across the junction, respectively. All models include the readout resonator Hamiltonian given by $H_{res} = Ωa†a + Gn(a + a†)$, where Ω is the bare resonator frequency, G is the electrostatic coupling strength and a†(a) is the bosonic creation (annihilation) operator. Including $H_{res}$ ensures that dressing of the states due to transmon-resonator hybridization is taken into account³⁴,³⁵,⁴²,⁴³.

We obtain the parameter set ($E_C$, $E_J$, Ω, G) of the standard transmon model in equation (3) by solving the inverse eigenvalue problem (IEP)⁴⁴⁻⁴⁷ for the measured spectroscopy data (Methods). For the Köln sample, these data include the offset charge dispersion (additional data for different magnetic fields are given in Supplementary Section IID). We note that the IEP is the very same science problem that was historically solved to model the energy spectra of natural atoms and molecules (see for example refs. 48–50), which led to the discovery of the fine structure.

In Fig. 3a, we show that the standard transmon model in equation (3) fails to describe the measured frequency spectra for all samples. The observed deviations are much larger than the measurement imprecision, for which we can set a conservative upper bound on the order of 1 MHz. While the standard transmon model with two parameters can trivially match the $f_{01}$ and $f_{02}$ transitions, the measured $f_{03}$ can already deviate by more than 10 MHz. The deviations are positive for the KIT, ENS and Köln samples, while the IBM transmons mostly show negative deviations (Supplementary Section IC5). It is important to remark that other corrections, such as the stray inductance in the JJ leads, hidden modes coupled to the qubit, the coupling between qubits as present on the IBM multi-qubit device, or an asymmetry in the superconducting

energy gaps, while being relevant, cannot fully account for the measured discrepancy (Supplementary Section ID). Notably, similar deviations can be found in previously published transmon spectra⁴¹,⁵¹⁻⁵³, as we detail in Supplementary Fig. 4 and Supplementary Sections IC2 and IC4.

In Fig. 3b, we demonstrate that orders of magnitude better agreement with our measured spectra can be achieved by using the Josephson harmonics model:

$$H_{har} = 4E_C(n - n_g)^2 - \sum_{m \ge 1} E_{jm} \cos(m\phi) + H_{res}. \quad (4)$$

In general, the values $E_{jm}$ are a fingerprint of each junction's channel-transparency distribution $ρ(T)$ with many degrees of freedom. Here we consider two simplified models (further models are discussed in Supplementary Section IC): (1) a phenomenological model truncated at $E_{J4}$ (top panel) and (2) a mesoscopic model of tunnelling through a non-uniform oxide barrier (bottom panel). We note that the phenomenological $E_{J4}$ model guarantees agreement for the lowest four transitions (Methods), and while many samples have physically reasonable $E_{Jm}$ coefficients when truncating at $E_{J4}$, a few JJs require terms up to $E_{J6}$ (Supplementary Section IC3).

The mesoscopic model allows us to derive $ρ(T; \tilde{d}, σ)$ based on a Gaussian thickness distribution with average thickness $\tilde{d}$ and standard deviation $σ$ (Supplementary Section IB4). As a consequence, all Josephson harmonics for $m \ge 2$ are parameterized in terms of the two parameters $\tilde{d}$ and $σ$ according to
<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Standard transmon model fails to describe the measured frequency spectra. a, Differences between the frequencies $f_{0j}^{\text{model}}$ predicted by the standard transmon model in equation (3) and the measured transitions $f_{0j}^{\text{experiment}}$. The markers indicate the different experiments at KIT (red stars), ENS (yellow squares), Köln (blue diamonds) and IBM (green circles). For the KIT experiment, we show results for three successive CDs of the same sample (CD1–3, dark red to bright red, respectively). For the Köln experiment, we chose a set of measured transitions at a fixed magnetic field (blue arrow in Fig. 4a). For the IBM experiment, we show results for 20 qubits in the IBM Hanoi device, using different marker sizes and shades of green. Measurement imprecisions are on the order of 1 MHz and not visible in the figure. Note that the scale on the vertical axis is linear between ±100 MHz and logarithmic onward. Dashed and dotted lines are guides to the eye. b, Same as a, with $f_{0j}^{\text{model}}$ given by the Josephson harmonics Hamiltonian in equation (4). Top, model truncated at $E_{14}$. Bottom, mesoscopic</figcaption></figure>

model of tunnelling through an inhomogeneous AlO$_x$ barrier, where $E_{1m}(\vec{d}, \sigma)$ is parameterized in terms of the average barrier thickness $\vec{d}$ and the standard deviation $\sigma$ (equation (5); the fit values are listed in the table inset). c, Ranges of the Josephson harmonics ratios $|E_{1m}/E_{1j}|$ that are consistent with the measured spectra. The ranges are represented by coloured vertical bars using the same colouring as in a. For the IBM Hanoi device, we show the ranges for qubits 0–2 from left to right (ranges for the other qubits are shown in Supplementary Section IC3). The shaded grey area highlights the region between two limiting cases: the fully open quantum point contact with unit transparency and a homogeneous barrier with $T_n = 10^{-6}$ for all n. Turquoise markers on the vertical bars indicate the harmonics ratios calculated from the mesoscopic model, where the average thickness $\vec{d}$ and the standard deviation $\sigma$ are given in b. Inset, an Al-AlO$_x$–Al junction obtained from molecular dynamics simulations (Supplementary Fig. 25) with average barrier thickness $\vec{d} = 1.5$ nm (Fig. 2c).

$$E_{1m}(\vec{d}, \sigma) \propto \int_0^1 c_m(T) \rho(T; \vec{d}, \sigma) dT, \quad (5)$$

where the Fourier coefficients $c_m(T)$ (equation (2) and Fig. 2d) are weighted by the channel-transparency distribution $\rho(T; \vec{d}, \sigma)$. In this model, relatively large ratios $|E_{1m}/E_{1j}|$ originate from higher-transparency contributions from the narrower regions of the barrier (compare the STEM images in Supplementary Fig. 27). The model can describe the samples at KIT, ENS and Köln (Fig. 3b) but not the IBM device (Supplementary Section IB4). The model parameters $\vec{d}$ and $\sigma$ (Fig. 3b) are comparable to results from molecular dynamics simulation and STEM pictures of the oxide barrier (Supplementary Section IV).

In Fig. 3c, we indicate the ranges of $E_{1m}$ coefficients consistent with the measured spectra. The bars represent the lower and upper limits of Josephson harmonics ratios $|E_{1m}/E_{1j}|$. The corresponding $\sin(m\phi)$ contribution to the C$\phi$R is given by $m|E_{1m}/E_{1j}|$ (see Fig. 1 for the KIT sample). The ratios lie between two limiting cases spanning the physical regime (shaded grey area): (1) the upper limit, $|E_{1m}/E_{1j}| = 3/(4m^2 - 1)$, corresponds to an open quantum point contact—that is, one channel with $T=1$—and (2) the lower limit, $|E_{1m}/E_{1j}| \approx (T/4)^{m-1}m^{3/2}$, corresponds to a perfectly homogeneous low-transparency barrier ($T_n = T = 10^{-6}$ for all n). For the scanning routine, we include harmonics up to $E_{110}$ to obtain results within the physical regime and to see when truncation is possible (Methods). Remarkably, for all samples, the $E_{12}$ contribution is in the few percent range even after considering additional corrections such as series inductance or gap asymmetry in the superconducting electrodes (Supplementary Section ID).

The Josephson harmonics ratios computed from the mesoscopic model in equation (5) are shown with turquoise markers. Notice that the barrier evolved between CDs of the KIT sample due to ageing (CD1 to CD2) and thermal annealing (CD2 to CD3) (Supplementary Section IIIA). Even for the most homogeneous barrier (CD3), the second-harmonic contribution is $E_{12}/E_{11} \approx -2.4\%$, implying that there would be at least one conduction channel with a transparency $T \ge 0.29$ (Supplementary Section IA). The methodology presented in Fig. 3 can serve as a tool to

characterize Josephson harmonics and tunnel barrier homogeneity, independent of circuit design.

Since the charge dispersion increases for higher transmon levels (even for the standard transmon Hamiltonian$^{34}$; Fig. 2f) and is exponentially sensitive to the shape of the JJ potential (Fig. 2e), a natural question arises: what are the consequences of the Josephson harmonics on the transmon’s susceptibility to offset charges? In Fig. 4a, we show the measured charge dispersion $\delta f_{0j}$ of the Köln device for states $j=1, 2, 3$ versus the first transition frequency $f_{01}$, which is tuned by an in-plane magnetic field $B_{\parallel}$ of up to 0.4 T (see Supplementary Section IIIC for details). The charge dispersion predicted by the standard model (dashed grey lines) falls short of the measurements by a factor of 2–7 for the three measured transitions. In contrast, when using the Josephson harmonics model, the computed charge dispersion matches the data (blue lines). We emphasize that for both models, we use the same parameters as in the Fig. 3 analysis (that is, the standard model and the $E_{14}$ model) and vary the first Josephson energy to match the qubit frequency $f_{01}$ for different magnetic fields while keeping the $E_{1m}/E_{1j}$ ratios constant.

Interestingly, the presence of large Josephson harmonics, as in the case of the IBM qubits (Fig. 3c), can also reduce the charge dispersion, which directly decreases charge noise decoherence. We show evidence for this in Fig. 4b, on the first three IBM qubits, for which the charge dispersion of the qubit transition can be a factor of 4 lower than expected from the standard transmon model. This observation indicates a possible optimization route in which Josephson harmonics are engineered (for example, by shaping the channel transparencies or adding inductive elements in series) and the spectrum is steered towards regions of reduced charge dispersion and increased anharmonicity (Supplementary Fig. 8). A recent work$^{54}$ proposes a similar approach to engineer arbitrary-shaped C$\phi$Rs using networks of effective high-transparency JJs, each of which is a series of tunnel JJs.

The main reason for the failure of the standard transmon model in describing the charge dispersion (when fitted to $f_{01}$ and $f_{02}$) is that it misjudges the value of $E_j/E_c$. To quantify this effect, in Fig. 4c we plot the values of $E_{1j}/E_c$ from the Josephson harmonics model against the value
<figure><img src="image_7.png" /></figure>

**Fig. 4 | Influence of Josephson harmonics on the charge dispersion.**
a, Measured charge dispersion $\delta f_0/f_0$ (blue diamonds) of states $j=1, 2, 3$ for the experiment in Köln, plotted as a function of the $f_{01}$ frequency. All transition frequencies are tuned, as the Josephson energy is suppressed by up to 35% by means of an in-plane magnetic field $B_1$ swept to 0.4 T. The standard model in equation (3), shown in dashed grey lines, underestimates the charge dispersion by a factor of 2–7 (grey arrows), while the Josephson harmonics model in equation (4) plotted in solid blue overlaps the measured data. Note that both are computed with the same parameters used for Fig. 3; the Josephson energy is reduced with increasing magnetic field, and the other parameters such as the $E_{J1}/E_{J1}$ ratios are kept constant. The blue arrow indicates $f_{01} = 5.079$ GHz, corresponding to the dataset shown in Fig. 3. b, Evidence that Josephson harmonics can reduce the charge dispersion by an order of magnitude (grey arrows). The dashed grey lines represent the standard model predictions.

In contrast, the green bars show results from all Josephson harmonics models. The data correspond to IBM qubits 0–2 (green bars in Fig. 3c) for the levels $j=1, 2, 3, 4$; results for all other samples are shown in Supplementary Fig. 6. c, The values of $E_{J1}/E_C$ change compared to the standard model $E_{J1}/E_C$, which constitutes the main correction to the predicted charge dispersions in a and b. The bars represent the range of suitable ratios $E_{J1}/E_C$ (Fig. 3c) for the successive CDs of the KIT sample (red bars), the ENS sample (yellow bar), the Köln sample (blue diamonds, using the same colour coding as in s) and the IBM Hanoi device (green bars). The dashed diagonal line indicates the case in which the ratios $E_{J1}/E_C$ of the harmonics model and $E_{J1}/E_C$ of the standard model are equal. Inset, correction $(e^{\text{har}} - e^{\text{std}})/e^{\text{std}}$ to the relative charge dispersion $\epsilon = \delta f_{01}/f_{01}$ for fixed $E^{\text{std}}/E_C = E^{\text{har}}/E_C$ for the Köln sample, where $e^{\text{std}}$ is given by the standard charge dispersion$^{34}$ and $e^{\text{har}}$ is computed using the Josephson harmonics model.

of $E_{J1}/E_C$ from the standard model. Indeed, the $E_{J1}/E_C$ ranges for many of our measurements are not compatible with the standard model $E_{J1}/E_C$ ratio (dashed diagonal). We note that when evaluated for the same $E_{J1}/E_C$, the Josephson harmonics correction to the charge dispersion is relatively small (inset of Fig. 4c).

In summary, we have shown that for ubiquitous AlOₓ tunnel junctions, the microscopic structure, currently underappreciated in its complexity, causes level shifts and modifies the charge dispersion in superconducting artificial atoms. In order to fully describe the measured transmon energy spectra, we amend the standard sin φ Josephson CφR for tunnel junctions to include higher-order sin(mφ) harmonics, with the relative amplitude of the m = 2 term in the few percent range. We confirm this finding in various sample geometries from four different laboratories, and we argue that the source of the Josephson harmonics is the presence of relatively higher-transparency channels with $T \gg 10^{-6}$ in the AlOₓ tunnel barrier. The methodology shown here can reveal percent-level deviations from a sinusoidal CφR, which are hard to detect in more standard measurements based on asymmetric direct current superconducting quantum interference devices$^{35}$.

The observation of Josephson harmonics in tunnel junctions highlights the need to revisit established models for superconducting circuits. Our work directly impacts the design and measurement of transmon qubits and processors: As an illustration, we show that by engineering Josephson harmonics, the dephasing due to charge noise can be reduced by an order of magnitude without sacrificing anharmonicity. These results ask for future research studying the implications of Josephson harmonics and associated Andreev bound states in other tunnel-JJ-based circuits, for example fluxonium or generalized flux qubits$^{36}$.

In general, we expect the inclusion of the harmonics will refine the understanding of superconducting artificial atoms and will directly

benefit, among others, quantum gate and computation schemes relying on higher levels$^{57–63}$, quantum-non-demolition readout fidelities$^{64–66}$ and frequency crowding mitigation in quantum processors$^{67}$. Josephson harmonics will probably also have to be accounted for in topological JJ circuits$^{68–70}$, parametric pumping schemes employed in microwave amplifiers and bosonic codes$^{71,72}$, amplification and mixing$^{7,73,74}$, JJ metrological devices$^{8–10}$, Floquet qubits$^{75,76}$, protected Josephson qubits$^{68,70,77}$ and so on, and they can be harnessed to realize Josephson diodes$^{78}$. As devices become increasingly sophisticated with progressively smaller error margins, higher-order Josephson harmonics will need to be either suppressed via the development of highly uniform and low-transparency barriers or engineered and included as an integral part of the device physics.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41567-024-02400-8.

## References

1. Josephson, B. D. Possible new effects in superconductive tunnelling. Phys. Lett. **1**, 251 (1962).

2. Josephson, B. D. The discovery of tunnelling supercurrents. Rev. Mod. Phys. **46**, 251 (1974).

3. Google Quantum AI. Exponential suppression of bit or phase errors with cyclic error correction. Nature **595**, 383–387 (2021).

4. Krinner, S. et al. Realizing repeated quantum error correction in a distance-three surface code. Nature **605**, 669–674 (2022).
5. Sivak, V. V. et al. Real-time quantum error correction beyond break-even. *Nature* **616**, 50–55 (2023).

6. Wendin, G. Quantum information processing with superconducting circuits: a perspective. Preprint at https://arxiv.org/abs/2302.04558 (2023).

7. Roy, A. & Devoret, M. Introduction to parametric amplification of quantum signals with Josephson circuits. *C. R. Phys.* **17**, 740–755 (2016).

8. Gramm, K., Lundgren, L. & Beckman, O. SQUID magnetometer for mangetization measurements. *Phys. Scr.* **13**, 93 (1976).

9. Shapiro, S. Josephson currents in superconducting tunneling: the effect of microwaves and other observations. *Phys. Rev. Lett.* **11**, 80 (1963).

10. Crescini, N. et al. Evidence of dual Shapiro steps in a Josephson junction array. *Nat. Phys.* **19**, 851–856 (2023).

11. Albertinale, E. et al. Detecting spins by their fluorescence with a microwave photon counter. *Nature* **600**, 434–438 (2021).

12. Arute, F. et al. Quantum supremacy using a programmable superconducting processor. *Nature* **574**, 505–510 (2019).

13. King, A. D. et al. Quantum critical dynamics in a 5,000-qubit programmable spin glass. *Nature* **617**, 61–66 (2023).

14. Kim, Y. et al. Evidence for the utility of quantum computing before fault tolerance. *Nature* **618**, 500–505 (2023).

15. Likharev, K. K. *Dynamics of Josephson Junctions and Circuits* (Gordon and Breach, 1986).

16. Golubov, A. A., Kupriyanov, M. Y. & Il'ichev, E. The current-phase relation in Josephson junctions. *Rev. Mod. Phys.* **76**, 411–469 (2004).

17. Goldobin, E., Koelle, D., Kleiner, R. & Buzdin, A. Josephson junctions with second harmonic in the current-phase relation: properties of φ junctions. *Phys. Rev. B* **76**, 224523 (2007).

18. de Lange, G. et al. Realization of microwave quantum circuits using hybrid superconducting-semiconducting nanowire Josephson elements. *Phys. Rev. Lett.* **115**, 127002 (2015).

19. Kringhøj, A. et al. Anharmonicity of a superconducting qubit with a few-mode Josephson junction. *Phys. Rev. B* **97**, 060508 (2018).

20. Stoutimore, M. J. A. et al. Second-harmonic current-phase relation in Josephson junctions with ferromagnetic barriers. *Phys. Rev. Lett.* **121**, 177702 (2018).

21. Bargerbos, A. et al. Observation of vanishing charge dispersion of a nearly open superconducting island. *Phys. Rev. Lett.* **124**, 246802 (2020).

22. Dolan, G. J. Offset masks for lift-off photoprocessing. *Appl. Phys. Lett.* **31**, 337–339 (1977).

23. Beenakker, C. W. J. in *Transport Phenomena in Mesoscopic Systems* (eds Fukuyama, H. & Ando, T.) 235–253 (Springer Berlin Heidelberg, 1992); https://doi.org/10.1007/978-3-642-84818-6_22

24. Glazman, L. & Catelani, G. Bogoliubov quasiparticles in superconducting qubits. *SciPost Phys. Lect. Notes* https://doi.org/10.21468/SciPostPhysLectNotes (2021).

25. Kittel, C. *Introduction to Solid State Physics* (Wiley, 2004).

26. Fritz, S., Radtke, L., Schneider, R., Weides, M. & Gerthsen, D. Optimization of Al/AlOx/Al-layer systems for Josephson junctions from a microstructure point of view. *J. Appl. Phys.* **125**, 165301 (2019).

27. Da Costa, V., Tiusan, C., Dimopoulos, T. & Ounadjela, K. Tunneling phenomena as a probe to investigate atomic scale fluctuations in metal/oxide/metal magnetic tunnel junctions. *Phys. Rev. Lett.* **85**, 876–879 (2000).

28. Parkin, S. S. P. et al. Giant tunnelling magnetoresistance at room temperature with MgO (100) tunnel barriers. *Nat. Mater.* **3**, 862–867 (2004).

29. Yuasa, S., Nagahama, T., Fukushima, A., Suzuki, Y. & Ando, K. Giant room-temperature magnetoresistance in single-crystal fe/MgO/fe magnetic tunnel junctions. *Nat. Mater.* **3**, 868–871 (2004).

30. Cyster, M. J. et al. Effect of atomic structure on the electrical response of aluminum oxide tunnel junctions. *Phys. Rev. Res.* **2**, 013110 (2020).

31. Cyster, M. J. et al. Simulating the fabrication of aluminium oxide tunnel junctions. *npj Quantum Inf.* **7**, 12 (2021).

32. Aref, T. et al. Characterization of aluminum oxide tunnel barriers by combining transport measurements and transmission electron microscopy imaging. *J. Appl. Phys.* **116**, 073702 (2014).

33. Zeng, L. J. et al. Direct observation of the thickness distribution of ultra thin AlOₓ barriers in Al/AlOₓ/Al Josephson junctions. *J. Phys. D: Appl. Phys.* **48**, 395308 (2015).

34. Koch, J. et al. Charge-insensitive qubit design derived from the Cooper pair box. *Phys. Rev. A* **76**, 042319 (2007).

35. Blais, A., Grimsmo, A. L., Girvin, S. M. & Wallraff, A. Circuit quantum electrodynamics. *Rev. Mod. Phys.* **93**, 025005 (2021).

36. Nigg, S. E. et al. Black-box superconducting circuit quantization. *Phys. Rev. Lett.* **108**, 240502 (2012).

37. Ansari, M. H. Superconducting qubits beyond the dispersive regime. *Phys. Rev. B* **100**, 024509 (2019).

38. Riwar, R.-P. & DiVincenzo, D. P. Circuit quantization with time-dependent magnetic fields for realistic geometries. *npj Quantum Inf.* **8**, 36 (2022).

39. Miano, A. et al. Hamiltonian extrema of an arbitrary flux-biased Josephson circuit. *PRX Quantum* **4**, 030324 (2023).

40. Lescanne, R. et al. Escape of a driven quantum Josephson circuit into unconfined states. *Phys. Rev. Applied* **11**, 014030 (2019).

41. Krause, J. et al. Magnetic field resilience of three-dimensional transmons with thin-film Al/AlOₓ/Al Josephson junctions approaching 1 T. *Phys. Rev. Appl.* **17**, 034032 (2022).

42. Blais, A., Huang, R.-S., Wallraff, A., Girvin, S. M. & Schoelkopf, R. J. Cavity quantum electrodynamics for superconducting electrical circuits: an architecture for quantum computation. *Phys. Rev. A* **69**, 062320 (2004).

43. Blais, A. et al. Quantum-information processing with circuit quantum electrodynamics. *Phys. Rev. A* **75**, 032329 (2007).

44. Friedland, S. Inverse eigenvalue problems. *Linear Algebra Appl.* **17**, 15–51 (1977).

45. Friedland, S., Nocedal, J. & Overton, M. L. The formulation and analysis of numerical methods for inverse eigenvalue problems. *SIAM J. Numer. Anal.* **24**, 634–667 (1987).

46. Chu, M. T. Inverse eigenvalue problems. *SIAM Rev.* **40**, 1–39 (1998).

47. Chu, M. & Golub, G. Inverse Eigenvalue Problems: Theory, Algorithms, and Applications (Oxford Univ. Press, 2005); https://doi.org/10.1093/acprof:oso/9780198566649.001.0001

48. Downing, A. C. & Householder, A. S. Some Inverse Characteristic Value Problems Which Arise in the Study of Simple Molecules. Technical report no. CF-55-10-95 (Oak Ridge National Laboratory, 1955); https://doi.org/10.2172/4357461

49. Toman, S. & Plíva, J. Multiplicity of solutions of the inverse secular problem. *J. Mol. Spectrosc.* **21**, 362–371 (1966).

50. Brussaard, P. J., Glaudemans, P. W. M. & Klein, A. Shell Model Applications in Nuclear Spectroscopy (North-Holland, 1977); https://doi.org/10.1002/piuz.19790100516

51. Peterer, M. J. et al. Coherence and decay of higher energy levels of a superconducting transmon qubit. *Phys. Rev. Lett.* **114**, 010501 (2015).

52. Xie, E. et al. Compact 3D quantum memory. *Appl. Phys. Lett.* **112**, 202601 (2018).

53. Xie, E. Scalable 3D Quantum Memory. Ph.D. thesis, Technische Universität München (2019).

54. Bozkurt, A. M., Brookman, J., Fatemi, V. & Akhmerov, A. R. Double-fourier engineering of Josephson energy-phase relationships applied to diodes. *SciPost Phys.* **15**, 204 (2023).
55. Della Rocca, M. L. et al. Measurement of the current-phase relation of superconducting atomic contacts. *Phys. Rev. Lett.* **99**, 127005 (2007).

56. Nguyen, L. B. et al. Blueprint for a high-performance fluxonium quantum processor. *PRX Quantum* **3**, 037001 (2022).

57. DiCarlo, L. et al. Demonstration of two-qubit algorithms with a superconducting quantum processor. *Nature* **460**, 240–244 (2009).

58. Rigetti, C. & Devoret, M. Fully microwave-tunable universal gates in superconducting qubits with linear couplings and fixed transition frequencies. *Phys. Rev. B* **81**, 134507 (2010).

59. Chow, J. M. et al. Microwave-activated conditional-phase gate for superconducting qubits. *New J. Phys.* **15**, 115012 (2013).

60. Sank, D. et al. Measurement-induced state transitions in a superconducting qubit: beyond the rotating wave approximation. *Phys. Rev. Lett.* **117**, 190503 (2016).

61. Negirneac, V. et al. High-fidelity controlled-z gate with maximal intermediate leakage operating at the speed limit in a superconducting quantum processor. *Phys. Rev. Lett.* **126**, 220502 (2021).

62. Rol, M. A. et al. Fast, high-fidelity conditional-phase gate exploiting leakage interference in weakly anharmonic superconducting qubits. *Phys. Rev. Lett.* **123**, 120502 (2019).

63. Roy, T., Li, Z., Kapit, E. & Schuster, D. Two-qutrit quantum algorithms on a programmable superconducting processor. *Phys. Rev. Appl.* **19**, 064024 (2023).

64. Gusenkova, D. et al. Quantum nondemolition dispersive readout of a superconducting artificial atom using large photon numbers. *Phys. Rev. Applied* **15**, 064030 (2021).

65. Shillito, R. et al. Dynamics of transmon ionization. *Phys. Rev. Appl.* **18**, 034031 (2022).

66. Cohen, J., Petrescu, A., Shillito, R. & Blais, A. Reminiscence of classical chaos in driven transmons. *PRX Quantum* **4**, 020312 (2023).

67. Hertzberg, J. B. et al. Laser-annealing Josephson junctions for yielding scaled-up superconducting quantum processors. *npj Quantum Inf.* **7**, 129 (2021).

68. Gyenis, A. et al. Experimental realization of a protected superconducting circuit derived from the $0-\pi$ qubit. *PRX Quantum* **2**, 010339 (2021).

69. Rymarz, M., Bosco, S., Ciani, A. & DiVincenzo, D. P. Hardware-encoding grid states in a nonreciprocal superconducting circuit. *Phys. Rev. X* **11**, 011032 (2021).

70. Smith, W. C. et al. Magnifying quantum phase fluctuations with Cooper-pair pairing. *Phys. Rev. X* **12**, 021002 (2022).

71. Campagne-Ibarcq, P. et al. Quantum error correction of a qubit encoded in grid states of an oscillator. *Nature* **584**, 368–372 (2020).

72. Lescanne, R. et al. Exponential suppression of bit-flips in a qubit encoded in an oscillator. *Nat. Phys.* **16**, 509–513 (2020).

73. Frattini, N. E. et al. 3-wave mixing Josephson dipole element. *Appl. Phys. Lett.* **110**, 222603 (2017).

74. Grimm, A. et al. Stabilization and operation of a Kerr-cat qubit. *Nature* **584**, 205–209 (2020).

75. Huang, Z. et al. Engineering dynamical sweet spots to protect qubits from 1/f noise. *Phys. Rev. Appl.* **15**, 034065 (2021).

76. Gandon, A., Le Calonnec, C., Shillito, R., Petrescu, A. & Blais, A. Engineering, control, and longitudinal readout of floquet qubits. *Phys. Rev. Appl.* **17**, 064006 (2022).

77. Gladchenko, S. et al. Superconducting nanocircuits for topologically protected qubits. *Nat. Phys.* **5**, 48–53 (2008).

78. Fominov, Y. V. & Mikhailov, D. S. Asymmetric higher-harmonic SQUID as a Josephson diode. *Phys. Rev. B* **106**, 134514 (2022).

79. Catelani, G. Parity switching and decoherence by quasiparticles in single-junction transmons. *Phys. Rev. B* **89**, 094522 (2014).

80. Serniak, K. et al. Direct dispersive monitoring of charge parity in offset-charge-sensitive transmons. *Phys. Rev. Appl.* **12**, 014052 (2019).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024, corrected publication 2024

<footnote n="1">¹Jülich Supercomputing Centre, Forschungszentrum Jülich, Jülich, Germany. ²IQMT, Karlsruhe Institute of Technology, Eggenstein-Leopoldshafen, Germany. ³PHI, Karlsruhe Institute of Technology, Karlsruhe, Germany. ⁴Departments of Applied Physics and Physics, Yale University, New Haven, CT, USA. ⁵Yale Quantum Institute, Yale University, New Haven, CT, USA. ⁶AIDAS, Jülich, Germany. ⁷Physics Institute II, University of Cologne, Köln, Germany. ⁸LPENS, Mines Paris-PSL, ENS-PSL, Inria, Université PSL, CNRS, Paris, France. ⁹Alice & Bob, Paris, France. ¹⁰IBM Quantum, IBM T. J. Watson Research Center, Yorktown Heights, NY, USA. ¹¹PGI-9, Forschungszentrum Jülich and JARA Jülich-Aachen Research Alliance, Jülich, Germany. ¹²RWTH Aachen University, Aachen, Germany. ¹³CETATEA, INCDTIM, Cluj-Napoca, Romania. ¹⁴PGI-2, Forschungszentrum Jülich, Jülich, Germany. ¹⁵PGI-11, Forschungszentrum Jülich, Jülich, Germany. ¹⁶Quantum Research Center, Technology Innovation Institute, Abu Dhabi, UAE. ¹⁷These authors contributed equally: Dennis Willsch, Dennis Rieger. ✉e-mail: ioan.pop@kit.edu</footnote>
## Methods

### Diagonalizing the Hamiltonians to obtain model predictions

We construct the matrices of $H_{\text{std}}$ in equation (3) and $H_{\text{har}}$ in equation (4) by first diagonalizing the bare transmon matrix (excluding $H_{\text{res}}$) in the charge basis $|n\rangle$, where $4E_C(n-n_g)^2 = \sum_n 4E_C(n-n_g)^2 |n\rangle\langle n|$ is diagonal and $-E_{lm}\cos(m\phi) = -\sum_n E_{lm}/2(|n\rangle\langle n+m| + |n+m\rangle\langle n|)$ has constant entries $-E_{lm}/2$ on the $m$th subdiagonal (we ensure enough terms by generally verifying that the predictions do not change if more terms are included). This yields the transmon eigenenergies $E_l$ and eigenstates $|j\rangle$. Then we diagonalize the joint transmon-resonator Hamiltonian $H_{\text{std/har}} = \sum_l E_l |j\rangle\langle j| + \Omega a^\dagger a + \sum_{i,j,l} G |j\rangle\langle j| n |j'\rangle\langle j'| (a + a^\dagger)$, where $a = \sum_k \sqrt{k+1}|k\rangle\langle k+1|$. To each resulting eigenenergy $E_l$ and eigenstate $|l\rangle$, we assign a photon label $k$ and a transmon label $j$ based on the largest overlap $\max_{k,j} \langle kj|l\rangle$ (this only works for small $k$; Supplementary Section IIC), which yields the dressed energies $E_{kj}$ and states $|kj\rangle$. This procedure is done for both $n_g=0$ and $n_g=1/2$. From the resulting dressed energies $E_{kj}(n_g)$, we compute the transmon transition frequencies $f_{0j}^{\text{model}}(n_g) = (E_{0j}(n_g) - E_{00}(n_g))/2\pi$ and the resonator frequencies $f_{res,j}^{\text{model}}(n_g) = (E_{lj}(n_g) - E_{0j}(n_g))/2\pi$ (setting $\hbar=1$). The predicted frequencies are then given by $f_{0j}^{\text{model}} = (f_{0j}^{\text{model}}(0) + f_{0j}^{\text{model}}(1/2))/2$, $f_{res,j}^{\text{model}} = (f_{res,j}^{\text{model}}(0) + f_{res,j}^{\text{model}}(1/2))/2$, and the charge dispersion is $\delta f_{0j}^{\text{model}} = |f_{0j}^{\text{model}}(0) - f_{0j}^{\text{model}}(1/2)|$. We consistently use $n = -N, ..., N$ with $N=14$ and thus $2N+1=29$ charge states, $j=0, ..., M-1$ with $M=12$ transmon states and $k=0, ..., K-1$ with $K=9$ resonator states, where $N, M$ and $K$ have been chosen by verifying that the model predictions change by less than a few kHz when adding more states.

### Solving the IEP to obtain model parameters

The inverse problem$^{47,81}$ to obtain the parameters $\mathbf{x}^{\text{std}}$ of the standard model Hamiltonian in equation (3) and $\mathbf{x}^{\text{har}}$ of the harmonics model Hamiltonian in equation (4), such that the linear combinations of eigenvalues $\mathbf{f} = (f_{01}^{\text{model}} \ f_{02}^{\text{model}} \ \dots \ f_{N_1}^{\text{model}} \ f_{\text{res},0}^{\text{model}} \ f_{\text{res},1}^{\text{model}})$ agree with the measured data, is an instance of the Hamiltonian parameterized IEP (HamPIEP; Supplementary Section IIA2). We solve the HamPIEP using the globally convergent Newton method$^{82}$ with cubic line search and backtracking$^{83}$ and the Broyden-Fletcher-Goldfarb-Shanno algorithm$^{84}$ as implemented in TensorFlow Probability$^{85,86}$. The Jacobian $\partial\mathbf{f}/\partial\mathbf{x}$ is obtained by performing automatic differentiation through the diagonalization with TensorFlow. For the $E_{j4}$ model shown in Fig. 3b, the IEP is solved unambiguously for $\mathbf{x} = (E_{j1}, E_{j2}, E_{j3}, E_{j4}, \Omega, G)$ using the lowest four transmon transition frequencies, and we fix the values $E_C^{\text{KIT}}/h = 242$ MHz, $E_C^{\text{ENS}}/h = 180$ MHz and $E_C^{\text{IBM}}/h = 300$ MHz, respectively, to make the models consistent with further available information such as accurate finite-element simulations (Supplementary Section IIIA) or knowledge of the transmon capacitance. For the mesoscopic model (Supplementary Section IB4), the parameters $\mathbf{x} = (\tilde{d}, \sigma, E_C, E_j, \Omega, G)$ are found by minimizing the function $\sum_{j=1}^{N_1} |f_{0j}^{\text{model}}/j - f_{0j}^{\text{experiment}}/j| + \sum_{j=0}^{N_2} |f_{\text{res},j}^{\text{model}} - f_{\text{res},j}^{\text{experiment}}|$ using the Broyden-Fletcher-Goldfarb-Shanno algorithm. The initial values for the minimization are given by $\tilde{d} = 1.64$ nm (taken from the molecular dynamics result in Supplementary Section IV), $\sigma = \tilde{d}/4$ (also Supplementary Table 2) and $(E_C, E_j, \Omega, G)$ from the standard transmon model. For the Köln data, where 288 data points have to be described by the same model parameters $\mathbf{x}$ (Fig. 4a) and only the Josephson energy is varied, we use cubic interpolation as a function of $f_{01}^{\text{model}}$ and include only a few central points for the available frequencies in the solution of the IEP (the residuals are given in Supplementary Fig. 17). All model parameters are provided in the repository$^{87}$ accompanying this manuscript.

### Scanning the Josephson energies

To obtain the range of suitable Josephson energies $\{E_{lm}\}$ (shown in Fig. 3c) that are consistent with a measured spectrum, we use an exhaustive scanning procedure. A spectroscopy dataset of $N_1$ measured transition frequencies $f_{0j,l}=1, ..., N_1$ and two resonator frequencies $f_{\text{res},0}$ and $f_{\text{res},1}$ uniquely determines—via the HamPIEP—the values

$\mathbf{x} = (E_{j1}, ..., E_{jN_1}, \Omega, G)$. We then scan the values of four additional ratios $\mathbf{y} = (E_{jN_1+1}/E_{j1}, ..., E_{jN_1+4}/E_{j1})$ for each of these four $E_{jlm}/E_{j1}$ over 16 geometrically spaced values between the point contact limit $3(-1)^{m+1}/(4m^2-1)$ and $(-1)^{m+1} \min\{10^{-7}, |E_{jlm+1}/E_{j1}|\}$ (always skipping the first to ensure $|E_{jlm}/E_{j1}| > |E_{jlm+1}/E_{j1}|$). Additionally, we include $\mathbf{y} = (0, 0, 0, 0)$ to see if truncation at $E_{jN_1}$ is allowed. For each combination $\mathbf{y}$, we solve the HamPIEP for the spectroscopy data to obtain the unique solution $\mathbf{x}$. We call the ratios $\mathbf{e} = (1, E_{j12}/E_{j1}, ..., E_{jN_1+4}/E_{j1})$ a trajectory that can reproduce the spectrum. However, the trajectory $\mathbf{e}$ may not be physical, since (1) some of the leading ratios $E_{jlm}/E_j$ for $m \le N_1$ might be beyond the quantum point-contact limit, (2) the Josephson energies might not be strictly decreasing in absolute value for increasing order $m$, or (3) the signs might not be alternating. Note that this can also happen when the Josephson harmonics model in equation (4) is truncated at too-low orders (Supplementary Section IC3). For all $E_{jlm}$, the maximum and minimum possible ratios $|E_{jlm}/E_j|$ define the vertical bars in Fig. 3c.

## Data availability

The spectroscopy data and the model parameters that support the findings of this study are available in the Jülich DATA repository at https://doi.org/10.26165/JUELICH-DATA/LGRHUH.

## References

81. Jaynes, E. T. & Bretthorst, G. L. *Probability Theory: The Logic of Science* (Cambridge Univ. Press, 2003); https://doi.org/10.1017/CBO9780511790423

82. Dennis, J. E. & Schnabel, R. B. *Numerical Methods for Unconstrained Optimization and Nonlinear Equations* (Society for Industrial and Applied Mathematics, 1996); https://doi.org/10.1137/1.9781611971200

83. Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. *Numerical Recipes: The Art of Scientific Computing* 3rd edn (Cambridge Univ. Press, 2007).

84. Nocedal, J. & Wright, S. J. *Numerical Optimization* (Springer, 2006); https://doi.org/10.1007/b98874

85. TensorFlow Developers. TensorFlow v.2.8.0 https://doi.org/10.5281/zenodo.5949125 (2022).

86. TensorFlow Developers. TensorFlow Probability v.0.15.0 https://github.com/tensorflow/probability/releases/tag/v0.15.0 (2022).

87. Willsch, D. & Rieger, D. Spectroscopy data and model parameters for 'Observation of Josephson harmonics in tunnel junctions'. Jülich DATA https://doi.org/10.26165/JUELICH-DATA/LGRHUH (2023).

88. Jülich Supercomputing Centre JUWELS cluster and booster: exascale pathfinder with modular supercomputing architecture at Juelich Supercomputing Centre. JLSRF 7, A138 (2021).

## Acknowledgements

D.W., M.W. and K.M. thank H. De Raedt and H. Lagemann for stimulating discussions. C.D., J.K. and Y.A. thank L.M. Janssen for assistance with the measurements and P. Janke for contributions to the data analysis. I.M.P. thanks U. Vool for providing comments on an early version of the manuscript. C.D. and I.M.P. thank L. DiCarlo for discussions that helped improve the manuscript. L.B.-I., C.M. and I.M.P. thank J. Cole for insightful discussions about simulating the $\text{AlO}_x$ structure. D.W., M.W. and K.M. gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer JUWELS$^{88}$ at Jülich Supercomputing Centre. D.W. and M.W. acknowledge support from the project Jülich UNified Infrastructure for Quantum computing, which that has received funding from the German Federal Ministry of Education and Research (BMBF) and the Ministry of Culture and Science of the State of North Rhine-Westphalia. P.W., B.D., T.R. and G.C. acknowledge support from the BMBF within the project GEQCOS (FKZ: 13N15683 and 13N15685).
D.R., S. Geisert, S. Günzler, S.I., W.W., G.C. and I.M.P. acknowledge
support from the BMBF within the project QSolid (FKZ: 13N16151
and 13N16149). C.D., J.K. and Y.A. acknowledge support from the
European Research Council under the European Union's Horizon
2020 research and innovation programme (grant no. 741121) and
from Germany's Excellence Strategy—Cluster of Excellence Matter
and Light for Quantum Computing EXC 2004/1-390534769. P.P.
acknowledges support from the BMBF within the QUANTERA project
SiUCs (FKZ: 13N15209). P.D. was supported by the IBM Quantum
Community Advocate internship programme. R.H., J.H.B., P.S. and D.G.
acknowledge the support of Hitachi High-Technologies. This work
has been supported financially by the BMBF via the TLE4HSQ project
(grant no. 13N15983). P.S. acknowledges financial support by the BMBF
via the Quantum Futur project MajoranaChips (grant no. 13N15264)
within the funding programme Photonic Research Germany. L.B.-I.
and C.M. acknowledge support from UEFISCDI Romania through the
contract ERANET-QUANTERA QuCos 120/16.09.2019 and from ANCS
through Core Program 27N/2023, project no. PN 23 24 01 04.

**Author contributions**

D.W., D.R., P.W., M.W. and I.M.P. conceived of the presented study.
D.W., D.R., P.W., M.W. and I.M.P. wrote the original draft. D.R., P.W.,
B.D., S. Günzler, P.P. and T.R. performed the experiments on the KIT
sample. C.D. and J.K. performed the experiments on the Köln sample.
R.L. and Z.L. performed the experiments on the ENS sample. N.T.B.
and P.D. performed the experiments on the IBM sample. R.H., J.H.B.,
P.S., S. Geisert and S.I. acquired the STEM images and prepared the

corresponding samples. D.W., M.W. and G.C. performed the theoretical
modelling and the numerical simulations. L.B.-I. and C.M. performed
the molecular dynamics simulations. D.P.D., K.M., G.C. and I.M.P.
supervised the work. All authors analysed the data and contributed
to reviewing and editing the manuscript and the Supplementary
Information.

**Funding**

Open access funding provided by Karlsruher Institut für
Technologie (KIT).

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version
contains supplementary material available at
https://doi.org/10.1038/s41567-024-02400-8.

**Correspondence and requests for materials** should be addressed to
loan M. Pop.

**Peer review information** *Nature Physics* thanks the anonymous
reviewers for their contribution to the peer review of this work.

**Reprints and permissions information** is available at
www.nature.com/reprints.