Stringent test of QED with hydrogen-like tin – Nature

Ion production

The ions produced in the EBIT are already rather slow (7 keV × N_q, with N_q as the charge state), which simplifies the capture of ions after ejection. In a pulsed drift tube, the kinetic energy is reduced to roughly 500 eV × N_q, which is low enough that we can capture the ion bunch within the capture section of the trap. A detailed description about the ion capture can be found in ref. ⁹. For the injection from the Heidelberg EBIT, a lot of time was spent optimizing the production efficiency as well as the ion transport through the beamline. At present, the largest limit in the amount of ions that can be trapped in our experiment is the transport efficiency. From the EBIT, several tens of thousand hydrogen-like Sn-118 are ejected. In a single shot, because of poor transport efficiency in the beamline, at most a handful of particles can be trapped in our setup. The particles get lost during transport through the beamline. The phase-space distribution in the EBIT is rather large owing to the high temperature in the plasma. This results in a large spread of the ejected ions, which makes it difficult to efficiently transport and trap the whole ion bunch. Most ions are lost at the pulsed drift tube, in which only a fraction of the bunch is decelerated. As this is a technical limitation, improvement in the transport efficiency would also be feasible in the future. But because we perform the measurement mostly on single particles, we only have to trap a few particles, assuming that the vacuum of the trap is good. In total, the required setup time is about two weeks, with a following measurement time of several months. For future measurements in even heavier hydrogen-like systems, such as hydrogen-like lead, the electron-beam energy has to be increased further. As of now, the Heidelberg EBIT is limited to around 65 kV of acceleration voltage, making the ionization of higher binding energies impossible. When this is overcome, either by improving the Heidelberg EBIT or setting up an EBIT with enough electron-beam energy, the g-factor measurement is expected to be straightforward, as the measurement scheme for tin could also be used for hydrogen-like lead or uranium.

Measurement scheme

After capture, the next step is to perform the high-precision measurement and corresponding systematic checks. As this normally takes a few weeks, the vacuum in our trap has to be extremely good. This is ensured by the cryogenic valve as shown in Fig. 1, which blocks the inflow of gas from the room-temperature beamline⁹. The valve and the cryogenic environment keep the vacuum below 10⁻¹⁶ mbar, allowing ion storage of many weeks.

In the presented Γ₀ measurement, a method similar to ref. ²⁶ was used. The measurement sequence is shown in Extended Data Fig. 1. First, the initial spin state is investigated in the analysis trap. After adiabatic transport to the precision trap and a waiting time of two minutes to allow the voltages to settle, the modified cyclotron frequency is determined by sideband coupling. This is done by measuring and fitting of the noise spectrum that the detection circuit coupled to the ion produces⁵⁶. Afterwards, the axial frequency is measured with the detection circuit as well; a spectrum is shown in Fig. 1d. This is followed by a ‘pulse and amplify’ (PnA) measurement sequence⁵⁷. In PnA, we determine the accumulated phase of the modified cyclotron motion after a fixed evolution time. By initial excitation of the mode to a certain radius, a free evolution time for undisturbed phase accumulation followed by readout of the accumulated phase, we can track the phase of the particle over time, allowing to extract the frequency with high precision. In this measurement, the set of evolution times consists of five reference phase measurements with 0.2-s evolution time, two unwrapping phases and two 5.2-s evolution times, which are used for the determination of the magnetic field. In the reference phases, the magnetic-field jitter is negligible on the phase stability, as this is not yet dominating. For the long phases, the magnetic-field stability is the dominating jitter/drift source. The two unwrapping phases 0.5 s and 2.2 s are measured to allow a consistent phase unwrapping to the final 5.2-s measurement time. All except the two 5.2-s measurements are performed in a random order. The 5.2-s measurement is done twice at the end of the cycle. The first is used to precisely determine the magnetic field, which is used to calculate the expected Larmor frequency for the microwave injection. During the second 5.2-s measurement, a microwave with a random frequency offset to the expected Γ is irradiated. With this, we measure the magnetic field one more time, while trying to perform a spin-flip with the microwave. For the final determination of the cyclotron frequency, we only use the five reference phases and the last 5.2-s phase during the microwave excitation. All others are used for unwrapping and to test certain systematics. After the PnA cycle, the axial frequency is measured a second time. The two axial-frequency measurements are used to identify systematic shifts in the measurement, although mostly the second is used for the calculation of the free-space cyclotron frequency, as it is immediately after the microwave injection. Afterwards, the ion is brought back to the analysis trap, in which we test if the precision trap microwave injection ‘flipped’ the spin. The measurement scheme is shown in Extended Data Fig. 1.

Resonance analysis

The resulting resonance was analysed with a maximum-likelihood fit. Of the 387 spin-flip tries in the precision trap, 54 have been successful. The precision trap is at a position at which the second-order magnetic-field inhomogeneity is smaller than 10 mT m⁻²; the influence on the field inhomogeneity can be safely neglected as the resulting error is less than 1 part per trillion (ppt)⁵⁸. Odd-order inhomogeneity effects led by the first-order B₁ ≈ 2.64(3) mT m⁻¹ (ref. ⁹) are even further suppressed, as the odd orders are cancelled in a harmonic trap. The remaining influences in the line shape are given by magnetic-field jitter and by the power of the microwave injection.

Because the microwave injection time (5 s) is longer than the time the spin stays coherent to the (weak) drive, the spin-flip probability can be at most 50%. If high microwave powers are used, the chance to drive a spin-flip increases accordingly and broadens the resonance. A microwave-power-dominated resonance would follow a Lorentzian line shape. On the other hand, if the power is small enough, the magnetic-field jitter dominates, the spin-flip probability drops below 50% and the line shape changes from a Lorentzian to a Gaussian distribution. For the described resonance, a Gaussian function is used:

with A being the amplitude, Γ the irradiated Larmor frequency divided by the measured cyclotron frequency, Γ_stat the centre of the resonance and σ the standard deviation of the normally distributed data. Using a maximum-likelihood fitting method, which uses the unbinned dataset of Γ ratios, we extract the parameters of the resonance. It has an amplitude of A = 29(4)% and a full width at half maximum of 5.6(3) × 10⁻¹⁰. Owing to the low statistics, a Lorentzian fit results in a rather similar likelihood, but because the resonance width is consistent with the expected jitter based on the PnA phase stability, this is unlikely to be the case. Further, the Lorentzian fit would give a smaller uncertainty on the resonance centre. Using the Gaussian shape is thus the conservative approach in the extraction of Γ_stat and its uncertainty. In Extended Data Fig. 2, the likelihood surfaces of the maximum-likelihood fit are shown. They show the maximum-likelihood planes in the three-dimensional parameter space spanned by the free parameters Γ_stat, σ and A.

Systematic effects

Of all the systematic effects, three dominate the uncertainty/correction. The largest shift is the image-charge shift. It is calculated according to ref. ⁵⁵ and amounts to a relative shift of 1.50(8) × 10⁻¹⁰. Here the error is 5% of the total value, following ref. ⁵⁵.

The systematic effect with the largest uncertainty on Γ₀ arises from the axial frequency fit. The axial resonator causes a frequency pushing depending on the relative frequency difference of the particle and the resonator. This is accounted for in the fit function, but because the frequency of the resonator cannot be measured more accurately than a couple of hertz, the resulting uncertainty in the line shape of the axial dip cannot be neglected. With a resonator-frequency uncertainty of around 2 Hz, the fitted axial frequency changes by about 20 mHz. This translates to a 2.0 × 10⁻¹¹ relative uncertainty on Γ₀.

The relativistic correction owing to the mass increase on a cyclotron radius of 12.8(13) μm is 2.4(5) × 10⁻¹¹ (ref. ⁵⁹). Here we assume an error of 10% on the cyclotron radius during PnA. This was also confirmed by measuring a resonance with 161 cycles, resulting in 18 spin-flips at a much higher cyclotron radius (102(10) μm). Owing to the stronger mass increase, this is shifted by roughly 1.6 ppb and is in agreement with the expected shift (see Extended Data Fig. 3). This resonance also tests for hypothetical systematic effects from the unwrapping, as—in this dataset—the longest PnA evolution time is set to 10.2 s, nearly twice as long as in the main resonance.

Owing to the exceptionally harmonic trapping field in our precision trap and a good magnetic homogeneity, the residual shift of the result from higher-order coefficients in the electromagnetic potential can be limited to less than 1 ppt.

Another uncertainty arises from the less frequent/less precise measurement of the magnetron frequency. This is only measured every several tens of cycles, as its influence on the free-space cyclotron frequency is small. Conservatively, we assume a maximum error of 0.3 Hz, which results in a relative uncertainty of 3.8 × 10⁻¹². An extensive table of the different systematic effects is given in the Extended Data Table 1.

Mass measurement

To exclude the AME value of the neutral ¹¹⁸Sn isotope as an error source to the g factor, we performed a direct mass measurement. We trapped a hydrogen-like ¹¹⁸Sn⁴⁹⁺ ion along with a hydrogen-like ¹²C⁵⁺ ion. Both have a relatively similar q/m value, which results in suppressed systematic effects in the measurement. The axial frequencies are separated by roughly 820 Hz for the same electric potential in the precision trap. To determine the mass, we measure the cyclotron frequency ratio (CFR), which allows to extract the mass ratio without precise knowledge of the magnetic field:

To determine the mass, we used two measurement methods with partially different systematic effects to cross-check the result. A sketch of the measurement schemes in shown in Extended Data Fig. 4. In both methods, we measure the ions in an interleaved manner, transporting each particle subsequently into the precision trap, while storing the other in a neighbouring section. In the main method, we use the PnA scheme to determine the modified cyclotron frequency without any line-shape uncertainty. Atypically, we perform this measurement with one particle on resonance and the other 820 Hz apart. That way, we can apply the same voltages to the electrodes and effects resulting from a shifted potential/ion position are avoided. For this, we put the carbon ion in resonance; as a result of lower charge, this couples less to the detector system. It would be virtually impossible to detect its image currents when detuned by that amount. But because of the stronger coupling in hydrogen-like tin, this is still detectable even 820 Hz detuned, allowing to measure the motional frequencies.

In the second method, we change the electrostatic potential to tune each axial frequency into resonance with the detector. We use sideband coupling, also known as a ‘double-dip’ measurement, to measure the modified cyclotron frequencies of each particle. To move each particle on resonance, the axial potential is changed by roughly 150 mV out of the initial V ≈ −59 V. As possible patch potentials on the electrodes are not necessarily symmetric, the trap centre might shift slightly owing to the different voltages. Combined with a magnetic-field gradient B₁ = 2.64(3) mT m⁻¹, the measured frequency ratio could be systematically shifted because of this possible misalignment.

Another possible systematic shift arises from the different line shapes of the ions on the detector. The dip width of hydrogen-like tin and hydrogen-like carbon differs by about a factor of ten. Thus, effects owing to an erroneous input parameter to the fitting model of the ‘double dip’ could result in a systematic effect, as each ion would be affected differently. In the PnA method, the modified cyclotron frequency is determined independent of the detector line shape, therefore removing the connected systematic uncertainties. Owing to the better understanding of the systematic effects in the PnA method, we use the ‘double-dip’ measurement solely as a cross-check of the measured mass. When considering all systematic effects and uncertainties, this agrees with the PnA measurement.

Note that, here, only five frequency ratios have been collected in the phase-sensitive measurement. Therefore, we conservatively use the standard deviation as our statistical uncertainty and not the standard deviation of the mean.

Correcting the measurement for relativistic effects and the image-charge shift, we get a final value for the mass ratio:

The two brackets represent the statistical and systematic uncertainties, respectively. The mass of the ¹²C⁵⁺ ion can be expressed in relation to the neutral carbon atom as unit of mass after accounting for the missing electrons and binding energies. This gives a mass value for the ¹²C⁵⁺ ion of 11.99725768029217(43)(8) u (refs. ^28,31). The brackets are the uncertainty of the binding energies and the electron-mass uncertainty, respectively. As the total uncertainty is less than 0.01 ppt, its influence on the resulting hydrogen-like tin mass can be safely neglected. From this, we can infer a value for the atomic mass of the hydrogen-like tin-118 ion:

The first bracket is the statistical uncertainty and the second bracket is the systematic uncertainty, which is dominated by the relativistic effect.

Binding energies and neutral mass

The AME2020 value for the mass of neutral ¹¹⁸Sn is 117.90160663(54) u, which has an uncertainty of 466 eV/c² (u = 9.3149410242(28) × 10⁸ eV/c²). With the measured mass of ¹¹⁸Sn⁴⁹⁺, we can improve the accuracy of the mass of the neutral atom through

Here m_e = 0.000548579909070(16) u is the electron rest mass³¹. ΔE is the energy required to ionize the 49 electrons from a neutral Sn atom and is theoretically calculated to be 132,746(5) eV. The final result is m(¹¹⁸Sn) = 117.901606974(56)(5) u, which is a factor of 9.5 more accurate than the previous best value.

The value ΔE is derived from the electron-binding-energy difference between neutral Sn and hydrogen-like Sn⁴⁹⁺. As the electron-binding energy of Sn⁴⁹⁺ is known to be 35,192.501(11) eV (ref. ²²), we only need to calculate the electron-binding energy of neutral tin. This is performed with an ab initio, fully relativistic, multiconfiguration Dirac–Hartree–Fock (MCDHF) together with a relativistic configuration interaction (RCI) method^60,61,62 implemented in the GRASP2018 code⁶². However, because the binding energy of the four outermost electrons of Sn has been experimentally determined to be 93.22(4) eV (ref. ²⁸), the binding energy of the ground state of Pd-like Sn⁴⁺ ([Kr]4d^10 1S₀) is calculated instead. We note that the ionization potential for the fifth electron is also known experimentally. However, the ground state of Sn⁵⁺ is an open-shell configuration with a total angular momentum of 5/2. Thus, it requires a much larger basis set and amount of computational power to achieve the same accuracy as for the closed-shell ion Sn⁴⁺.

Within the MCDHF scheme, the many-electron atomic-state function is constructed as a linear combination of configuration state functions (CSFs) with common total angular momentum (J), magnetic (M) and parity (P) quantum numbers: \(\left|\Gamma PJM\right\rangle ={\sum }_{k}{c}_{k}\left|{\gamma }_{k}PJM\right\rangle \). The CSFs |γ_kPJM⟩ are given as jj-coupled Slater determinants of one-electron orbitals and γ_k summarizes all the parameters needed to fully define the CSF, that is, the orbital occupation and coupling of single-electron angular momenta. Γ collectively denotes all the γ_k included in the representation of the atomic-state function. The mixing coefficients c_k and the radial orbital wavefunctions are obtained by solving the MCDHF equations self-consistently^60,61, including the Dirac–Coulomb Hamiltonian. After that, the RCI method is used to calculate the contributions from mass shift, transverse photon interactions and QED effects.

We start with a Dirac–Hartree–Fock (DHF) calculation, in which only the ground-state configuration is considered. This gives a binding energy of 167,973.14 eV, with a −4.42-eV correction from the finite nuclear size effect. Further calculation with RCI adds contributions of −0.57 eV from mass shift, −120.34 eV from frequency-independent transverse photon interactions (or Breit interactions), 1.16 eV from frequency-dependent transverse photon interactions and −79.08 eV from QED terms. To derive the electron correlation energy, the size of the CSF basis set is gradually expanded through single and double (SD) excitation of electrons from the ground-state configuration to high-lying virtual orbitals. This allows us to monitor the convergence of the correlation energy by adding and optimizing virtual orbitals layer by layer up to n = 10 (n is the principal quantum number), with all orbital angular momenta ranging from 0 to n − 1 being included.

Six terms may lead the error bar: the uncertainty in the nuclear parameters, the finite basis set, the uncounted higher-order electron correlations, the insufficient basis functions, the inaccurate estimations for QED corrections and the uncalculated QED effect to the mass shift. The uncertainty in the nuclear radius gives a 0.06-eV uncertainty to the corrections in the finite nuclear size effect. To exclude the error caused by the finite basis set, we extrapolate our calculated value to n = ∞. This results in a SD correlation energy of 67.26(23) eV for the ground state of Sn⁴⁺. As this analysis is performed under the RCI calculation, the contribution from Breit interactions is fully accounted for. The frequency-dependent transverse photon interaction cannot be included in the multiconfiguration calculations of GRASP2018, thus its uncertainty will be accounted for later together with other untreated minor effects when deriving the systematic errors. The relativistic mass-shift operator implemented in the GRASP2018 code is accurate to the order of (m_e/M)(αZ)⁴m_ec² (M is the mass of the nucleus). Therefore, it bears an uncertainty of 0.02 eV for the mass-shift correction.

In the GRASP2018 code, the QED corrections are estimated by means of a screened-hydrogenic approximation⁶². This correction is dominated by inner-shell electrons. With an absolute value of 79.08 eV for Sn⁴⁺, one already has a QED correction of 76.64 eV for Be-like Sn⁴⁶⁺. Fortunately, the QED effect for a nearby element, Be-like Xe⁵⁰⁺, has been known to a sub-eV accuracy through ab initio QED calculations^63,64. This allows us to infer the accuracy of our QED calculations: with a calculated value of 99.03 eV for Xe⁵⁰⁺, it is 0.85 eV larger than its ab initio result. Assuming the same relative deviation, we derive a QED correction of 78.40(68) eV for Sn⁴⁺. However, for such an ion, the many-electron QED effects are difficult to evaluate accurately. In the following, we will effectively include these contributions into the systematic errors.

The systematic errors caused by the uncalculated effects can be estimated from the ionization potentials of the outermost electrons of Sn. For example, the ionization potential of the 5s electron in Sn³⁺ ([Kr]4d¹⁰5s ²S_1/2) is determined to be 40.74(4) eV by experiment²⁸. With the calculated binding energy of Sn³⁺ under a similar SD excitation scheme, we derive an ionization potential of 40.07 eV, which is 0.67 eV smaller than the experimental value. This deviation mainly originated from the high-order correlation effects, many-electron QED effects, insufficient basis functions and frequency-dependent transverse photon interactions. Nevertheless, this deviation becomes smaller for highly charged ions⁴⁵. Therefore, one could conservatively assume that the corrections to the ionization potentials decrease linearly to 0.10 eV for Cu-like Sn²¹⁺ (because we have found that the deviation is already below 0.10 eV for the ionization potential of Cu-like Kr⁶⁺). For the ionization potentials of ions throughout Sn²¹⁺ to Sn⁴⁸⁺, we conservatively assume that they all have a correction of 0.10 eV. In total, the contribution from all unaccounted for terms is in the range 0–9.73 eV. To cover this whole range, we can add a correction of 4.86 eV to the total binding energy of Sn⁴⁺ and simultaneously assume a systematic error of 4.86 eV.

Finally, we arrive at a total binding energy of 167,847(5) eV for the ground state of Sn⁴⁺. The different contributions and their uncertainties are summarized in Extended Data Table 2. This gives a binding-energy difference of 132,748(5) eV between neutral Sn and hydrogen-like Sn⁴⁹⁺. Combining the measured mass for Sn⁴⁹⁺, we obtain

for the neutral tin-118 atom. The first bracket is the measurement uncertainty and the second is the uncertainty of the electron-binding energies.

Theory of the bound-electron g factor

The leading g-factor contribution was first calculated in ref. ³⁴. It is based on the approximation of an infinitely small and infinitely heavy nucleus. Therefore, corrections owing to the finite nuclear size and mass need to be taken into account. Furthermore, QED corrections contribute to the total g factor, just as in the case of the free electron. However, QED corrections in the case of the bound electron differ from the free-electron case. In Fig. 4, we highlight especially QED binding corrections, that is, the difference between QED corrections for bound and free electrons. In the following, we discuss all relevant contributions.

Free-electron contributions

Contributions to the free-electron g factor were taken from ref. ³¹, namely, the one-loop to five-loop QED contributions. Hadronic as well as electroweak contributions as given in ref. ³¹ are too small to be relevant in this work.

Nuclear corrections

The finite size (FS) correction to the g factor as given in Extended Data Table 3 was calculated for the two-parameter Fermi distribution using formulas and tabulated parameters from ref. ⁶⁵. The first uncertainty corresponds to the nuclear root mean square charge radius as given in ref. ⁶⁶. Note that the uncertainty corresponding to the number of digits specified for relevant parameters in ref. ⁶⁵ is much smaller than the radius uncertainty, 7 × 10⁻¹⁰ ≪ 2 × 10⁻⁸. The second uncertainty of the FS correction from the table is a combination of the uncertainty owing to the nuclear polarization, deformation and susceptibility, together with a conservative estimate of the nuclear model dependence. The model dependence is the leading uncertainty and expresses the difference of the FS corrections for the two-parameter Fermi and the homogeneously charged sphere model, again following ref. ⁶⁵. A direct calculation of the FS correction for the sphere model using semianalytic wavefunctions was consistent with the result from ref. ⁶⁵.

Calculating the FS correction analytically using formulas from refs. ^67,68, we find a disagreement with numerical results corresponding to three times the nuclear radius uncertainty for both the sphere and Fermi models. This suggests that further, higher-order contributions need to be determined for the analytic approach to be accurate at high Z, such as Z = 50.

Finally, we also calculated the FS correction using the GRASP code⁶². In the absence of results for ¹¹⁸Sn⁴⁹⁺ in refs. ^69,70, the nuclear polarization correction was estimated as zero following ref. ⁶⁵, with an uncertainty estimated as 50% of the nuclear model uncertainty of the FS correction. Our estimates for the extra uncertainty owing to nuclear deformation⁷¹ and nuclear susceptibility⁷² corrections are negligible.

For recoil calculations, we used the mass value presented in this paper. The leading recoil term of the first order in the mass ratio was calculated to all orders in Zα using formulas and tabulated parameters in ref. ⁷³. Recoil corrections of higher order in \(\frac{m}{M}\) were calculated to the leading order in Zα but exactly in the mass ratio^74,75,76. We find the result up to \({\mathcal{O}}\left({\left(\frac{m}{M}\right)}^{2}\right)\) (refs. ^77,78) to deviate by less than 1 × 10⁻¹⁴ from the all-order result.

Radiative recoil corrections were calculated using formulas from refs. ^{1,74,76,77,78}. So far, recoil corrections to the g factor have been derived only for the model of a point-like nucleus.

One-loop QED

Binding corrections at the one-loop level have been calculated to all orders in Zα. For the one-loop self-energy correction, we used the result from ref. ³⁶, which is based on the model of a point-like nucleus. We calculated the Uehling part of the electric-loop vacuum polarization diagram for the model of a point-like nucleus both numerically, using the Uehling potential from ref. ⁷⁹ and the bound-electron wavefunction perturbed by a constant external magnetic field as derived in ref. ⁸⁰, as well as analytically using formulas from ref. ⁸¹, with both results in excellent agreement. Values for the Wichmann–Kroll electric-loop vacuum polarization correction as well as the magnetic-loop vacuum polarization correction were taken from ref. ¹ and were calculated for extended nuclei in that work.

Combined QED-FS corrections

Bound-electron QED corrections, when carried out for the model of an extended nucleus, give slightly different results compared with point-nucleus calculations. Therefore, for an accurate theoretical description of the bound-electron g factor, we have to take into account the correction to QED contributions owing to the finite nuclear size (QED-FS corrections). As mentioned in the previous section, the Wichmann–Kroll part of the electric-loop vacuum polarization as well as the magnetic-loop vacuum polarization corrections already include the QED-FS corrections. We calculated the FS corrections to the one-loop self-energy correction and the Uehling contribution to the electric-loop vacuum polarization correction for the two-parameter Fermi distribution of the nuclear charge using formulas, and tabulated parameters for Z = 50, as given in ref. ⁶⁵. An older calculation of the FS correction to the one-loop self-energy contribution based on the spherical shell model of the nucleus from ref. ³⁶ differs by 2.1 × 10⁻⁹ from the result for the Fermi distribution. Using semianalytic wavefunctions for the homogeneous sphere model of the nucleus, we calculated the FS correction to the Uehling part of the electric-loop vacuum polarization correction⁸². Our result differs by 5 × 10⁻¹⁰ from the result for the Fermi distribution. We therefore assign (nuclear model) uncertainties of 2.1 × 10⁻⁹ to the one-loop self-energy FS correction, as well as 5 × 10⁻¹⁰ to the FS correction to the Uehling part of the electric-loop vacuum polarization correction.

Two-loop and higher QED diagrams

We calculated binding corrections to two-to-five-loop QED diagrams of order (Zα)² following ref. ⁸³. See also refs. ^74,77,78 for earlier derivations of these binding corrections. (Results given in the lines ‘2-loop QED, (Zα)²’ and ‘≥3-loop QED, binding’ in Extended Data Table 3.)

Two-loop binding corrections of order (Zα)⁴ were derived in refs. ^40,84. All-order calculations in Zα were carried out in ref. ³⁷ for a subset of two-loop diagrams, namely those diagrams with at least one vacuum polarization loop. However, magnetic-loop vacuum polarization diagrams were not considered in that work. In ref. ³⁷, results are given explicitly for the contribution of orders (Zα)⁵ and higher. We give this result in the line ‘(Zα)⁵⁺ S(VP)E, SEVP, VPVP’. For the remaining two-loop Feynman diagrams, QED corrections of order (Zα)⁵ were calculated in ref. ⁴⁰. In ref. ⁴⁰, a relative uncertainty of 13% is mentioned, owing to uncalculated Feynman diagrams contributing to order (Zα)⁵. This corresponds to our uncertainty in the line (Zα)⁵ in Extended Data Table 3. The uncertainty in the line ‘(Zα)⁵⁺ S(VP)E, SEVP, VPVP’ corresponds to higher-order corrections of order (Zα)⁶⁺ of two-loop Feynman diagrams with one vacuum polarization magnetic loop, by interpolating between tabulated results from ref. ³⁹ for nearby Z.

The uncertainties owing to uncalculated Feynman diagrams with two self-energy loops of order (Zα)⁶⁺ were estimated using the methods from refs. ^40,84, with the larger of the two estimates chosen as the uncertainty for these contributions in Extended Data Table 3. Furthermore, uncalculated binding corrections of \({\mathcal{O}}\left({(Z\alpha )}^{4+}\right)\) to Feynman diagrams with three and more loops were estimated by adapting the method in ref. ⁸⁴ (uncertainty in the line ‘≥3-loop QED, binding’).

As can be seen from Extended Data Table 3, the total uncertainty of the theoretical g-factor value is dominated by uncalculated higher order in Zα two-loop QED contributions. Calculations to improve the accuracy of two-loop corrections are underway^38,39. To improve the theoretical accuracy, an all-order (in Zα) calculation of the two-loop QED correction is required. Such calculations are underway^38,39. The most difficult part of the calculation is the two-loop self-energy, which is split into several parts according to the degree of their ultraviolet divergence, which are known as the loop-after-loop (LAL) correction and the F, M and P terms⁸⁵.

The F term is the part with overlapping ultraviolet divergences. It can be represented by Feynman diagrams with only free-electron propagators inside the self-energy loops and is evaluated in momentum space, thus avoiding any partial-wave expansion.

The M term is the ultraviolet finite part of the two-loop self-energy correction with Coulomb Dirac propagators inside the self-energy loops. Because the Coulomb Dirac propagator is best known in coordinate space, M-term calculations need to be carried out using a coordinate representation. Typically, M-term contributions are a double infinite sum of partial waves over angular momentum quantum numbers, which requires a very large number of partial waves to be calculated in practice. The calculation of every partial wave requires a multidimensional integration to be carried out numerically.

The P term is the part of the two-loop self-energy correction that contains Coulomb Dirac propagators in one part of the Feynman diagrams, as well as an ultraviolet divergent subdiagram. This requires P-term contributions to be calculated in a mixed coordinate–momentum representation, which involves the numerical Fourier transforms of the Coulomb Dirac propagators over one of the radial arguments. Details can be found in our earlier work³⁸.

Results for the so-called LAL and F-term contributions to two-loop self-energy corrections have been obtained, with the uncertainty of the LAL correction given in ref. ³⁸ for Z = 50 being 6.5 × 10⁻⁹ and the uncertainty of the F-term being orders of magnitude better. The calculation of the remaining parts of the two-loop self-energy correction, that is, the M and P terms, is continuing. A notable improvement of the total theoretical uncertainty compared with the value given in Extended Data Table 3 can be anticipated once these calculations are complete.

Muonic and hadronic vacuum polarization

The result for the muonic vacuum polarization correction given in Extended Data Table 3 corresponds to the Uehling part of the electric-loop contribution, calculated for the model of a point-like nucleus. Comparing results for point-like nuclei with results for extended nuclei from ref. ⁸⁶, we find good agreement between both for low Z. For Z = 70, the result for the extended nucleus is only about 50% of the point-nucleus result. For even higher Z, the extended nucleus result is much smaller than 50% of the point-nucleus result. For Z = 50, we therefore expect the muonic vacuum polarization correction for an extended nucleus to be larger than 50% of the point-nucleus result and assign a 50% uncertainty to the muonic vacuum polarization correction. We also calculated the Uehling part of the muonic vacuum polarization correction for the sphere model of the nucleus⁸². Our result of −2.0 × 10⁻⁹ is, within the specified uncertainty, in agreement with the point-nucleus result from Extended Data Table 3.

The hadronic vacuum polarization correction (‘Hadronic Uehling’) was estimated following refs. ^86,87 as 0.671 and 0.664 times the muonic vacuum polarization correction, with both estimates being identical to all digits given. We assigned the same uncertainty to the hadronic vacuum polarization correction that was used for the muonic vacuum polarization correction.