Testing Standard Model extensions with few-electron ions

When collecting spectroscopic data on at least four isotopes, nonlinearities in the King plot are a possible sign of Physics beyond the Standard Model. In this work, an improved approach to the search for hypothetical new interactions with isotope shift spectroscopy of few-electron ions is presented. Very careful account is taken of the small nuclear corrections to the energy levels and the gyromagnetic factors, which cause deviations from King linearity within the Standard Model and are hence a possible source of confounds. In this new approach, the experimental King nonlinearity is not compared to the vanishing prediction of the Standard Model at the leading order, but to the calculated full Standard Model contribution to King nonlinearity. This makes searching for beyond-the-Standard-Model physics with King linearity analysis possible in a very-high-precision experimental regime, avoiding confounds. The bounds which can be set on beyond-the-Standard-Model parameters remain limited by the uncertainties on the small Standard Model nuclear corrections which cause King nonlinearity. Direct comparison between theory and experiment on a single pair of isotopes is advocated as a more suitable approach for few-electron ions.


I. INTRODUCTION
Recent years have seen the accelerating development of tests of Standard Model (SM) extensions at the precision frontier [1][2][3][4][5][6][7][8][9][10][11]. In a context where the SM of fundamental physics is known to be incomplete, but where direct signals of new particles or interactions have so far not been reported, the precision frontier constitutes an interesting complement to the search for physics beyond the SM (also known as New Physics (NP)) at particle colliders, which probe the so-called energy and intensity frontiers. At the precision frontier, rigorous searches for NP signatures can be carried out with relatively smallsized experiments. In particular, in atomic physics, very precise experiments are often supported by calculations. Such high-precision calculations are possible in particular for few electron-ions, and the development of bound-state quantum electrodynamics (QED) has allowed predictions of fundamental spectroscopic quantities with relative precisions as good as 3 × 10 −15 in the case of an optical transition [12], and 2.8 × 10 −11 in the case of the boundelectron g (gyromagnetic) factor [13,14].
Comparisons between experimental and theoretical results on fundamental atomic spectroscopic quantities, such as transition frequencies and g factors, constitute stringent tests for QED in strong external fields (here, the electrostatic field of the nucleus). Such comparisons have also been used to determine fundamental constants, such as the electron rest mass [14], or the proton charge radius [15], with unprecedented precision. There also have been proposals for an improved determination of the finestructure constant α of electrodynamics [16][17][18][19][20]).
As has been shown in recent works, precision atomic spectroscopy can be used to set bounds on proposed SM * vincent.debierre@mpi-hd.mpg.de † natalia.oreshkina@mpi-hd.mpg.de ‡ zoltan.harman@mpi-hd.mpg.de extensions (conversely, it could be used to find signals of such beyond-SM phenomena). In particular, studies of the isotope shift of transition frequencies in singlycharged ions (and neutral atoms) have been shown [4-6, 8, 9] to allow for the setting of stringent bounds on NP, with "minimal theory inputs". Isotope shift data can be checked for deviations from King planarity which, with some care, could be interpreted as signatures of NP. Care is important, because departures from King planarity are also caused by small nuclear corrections to spectroscopic quantities within the SM [6][7][8]21]. A recent study on the transition frequencies of argon ions has shown that such effects had been underestimated in the past [7,21]. For a more systematically reliable interpretation of future experimental data on isotope shifts, it is important to estimate these small nuclear corrections accurately.
In the case of few-electron ions, calculations are generally more tractable than in the many-electron case, and hence can reach higher precisions. As such, these systems, and in particular their g factor, could become prime testing grounds for proposed SM extensions, with joint efforts from the theory and the experiment sides. On the theory side, more accurate calculations of the energy levels, the g factors, and their respective isotope shifts are made possible by the progress of bound-state QED. Both radiative QED contributions (two-loop radiative corrections) [22][23][24][25][26] and nuclear contributions [27][28][29][30][31][32] (we will discuss them in more detail throughout this work) are under active consideration, and are being determined with increasing precision. On the experiment side, the AL-PHATRAP group has reported the measurement of the g factor of B-like Ar 13+ as its first result [33], with a precision of 1.4 parts per billion, which represents an improvement on the previous value by seven orders of magnitude. It is projected [33] that the precision of future such measurements in Penning traps will outmatch the currently most precise reported results, which concern carbon and silicon (electron mass). Moreover, it is thought that not only g factors, but differences of g factors between two trapped ions, can be measured with a precision of one part in 10 11 and better [11]. This includes isotope shifts, which makes the present investigation very timely. Highprecision measurements of the g factor are also planned within the HITRAP project [34][35][36].
In this work, we study the setting of bounds on NP with few-electron ions, and show a simple approach based on direct comparison between theory and experiment to be best suited. This is in contrast to the case of many-electron ions, where a data-driven approach based on King representations and their generalisations is favoured [37]. We present in detail an advanced Kingtype formalism to set bounds on NP through isotope shift spectroscopy of few-electron ions by taking careful account of the aforementioned small nuclear corrections and their contribution to King nonplanarity, and show that this complex formalism is suboptimal for few-electron ions, compared to a simpler approach. The rest of this manuscript is organized as follows. In Sec. II, we briefly present proposed extensions of the SM which can be addressed through atomic spectroscopy, and through isotope shifts in particular. Radiative and inter-electronic-interaction QED corrections to the NP contributions to the g factor are discussed in Sec. III. In Sec. IV, the central part of this work, we study in detail the isotope shift of bound-electron energy levels and g factors, and build a modified King formalism to search for NP with spectroscopic data on three isotope pairs. It is shown that this indirect approach, while well-suited for many-electron systems, is unnecessary and even suboptimal for few-electron systems. In Sec. V, we argue that direct comparisons between experimental and theoretical results on few-electron ions can be used to set more competitive bounds on NP with a single isotope pair. Sec. VI is reserved for concluding remarks.

II. SOME EXTENSIONS OF THE STANDARD MODEL
Several proposed extensions of the SM can be probed with atomic physics experiments. We briefly discuss these extensions in the present section.

A. Higgs portal relaxion
The long-standing electroweak hierarchy problem [38], which is connected to the lack of an explanation of the weakness of gravity with respect to the other forces, concerns the mass of the Higgs boson, which could be expected to be much larger than it is due to coupling of the Higgs boson with any beyond-SM phenomena. Indeed, radiative corrections to the Higgs mass vary as the square of the energy scale at which beyond-SM phenomena might arise. The Higgs portal provides a solution to this problem. It involves the mixing of a new massive scalar boson, the relaxion, with the Higgs boson. These massive scalar bosons would mediate a new fundamental force, resulting, as far as atomic physics is concerned, in an interaction between nucleons and electrons [5,6,39]. The spin-independent potential exerted on electrons by this hypothetical force is of the Yukawa type [4]: where m φ is the mass of the relaxion, α HR = y e y n /4π is the Higgs portal coupling constant, with y e and y n the coupling of the massive scalar boson to the electrons and the nucleons, respectively, and c are Planck's reduced constant and the vacuum velocity of light, and A is the nuclear mass number of the considered ion. The first-order correction to the energy level of a bound electron in the quantum state a due to this potential corresponds to the diagram in Fig. 1. The contribution is Here g a and f a are the radial wave functions (large and small component, respectively) of the bound electron in state a [40]. For the H-like ground state a = 1s, we obtain Here and below, we use the point-like-nucleus approximation for wavefunctions and energies to deliver analytical formulas. The first-order correction to the g factor of a bound electron in the quantum state a due to this potential corresponds to the diagram in Fig. 1, together with the one in which the order of the two interactions is swapped. The contribution is FIG. 1. Feynman diagram corresponding to the leading contribution to the hypothetical fifth force correction to the energy level (E) and the g factor (g) of a bound electron. The double line represents the bound electron, the dashed line terminated by a square denotes the New Physics potential, and the wavy line terminated by a triangle denotes a photon from the external magnetic field. Diagram (g) has an equivalent diagram, as such, its contributions should be counted twice.
where µ 0 = e /2m e is the Bohr magneton, m a the magnetic projection quantum number of state a, B the magnitude of the external, static, homogeneous magnetic field, and X a and Y a the corrections to the large and small components of the bound electron radial wave function, due to the interaction with the magnetic field, given in Ref. [41]. With these wave function corrections and after performing the angular integration, Eq. (4) can be rewritten as with j a and κ a the total and Dirac angular momentum numbers, respectively, and λ e = /2πm e c the Compton wavelength of the electron. For the H-like ground state a = 1s, we obtain [11] g HR(1s It can be checked that the correction to the energy level and the g factor obey the relation which is derived in Ref. [42] for arbitrary central potentials. The full exact results for the a = 2s ground state of Li-like ions, and for the a = 2p 1/2 ground state of B-like ions, are given in the Supplemental Material to Ref. [7]. The difference B − L between the baryon and lepton numbers is conserved in the SM, and one proposed extension of the SM is built around the introduction of the U (1) B−L gauge symmetry. Unlike most group symmetries proposed to supplement the SM, this symmetry would be unbroken in the SM extension, giving rise to a new vector boson Z ′ , a massive hidden photon. Quantum anomalies in this model are canceled by the introduction of a right-handed neutrino for each lepton family, which introduces neutrino masses. In this SM extension, the boson Z ′ couples electrons to nucleons with a Yukawa potential [4,43], meaning that the results of Ref. [7], recapitulated in Sec. II A, are directly applicable, by replacing the Higgs portal coupling y e y n with the B − L symmetry coupling g 2 B−L . In this case, the couplings of the boson to the neutron and electrons is identical, therefore, atomic physics is sufficient to resolve the individual couplings of massive hidden photons to SM particles, which was not the case for the relaxions studied in Sec. II A.

C. Chameleon models
Chameleon particles have been proposed as dark energy candidates [44,45]. The most remarkable property of these scalar particles is that, due to their non-linear self-interaction, their mass is an increasing function of the energy density of the local environment. As such, in dense environments, these hypothetical particles would mediate a new force with a very limited range. Such models can be tested from atomic data. For practical purposes, the influence of this hypothetical particle on atomic spectra is described [3,46] by the potential Here M A is the mass of the nucleus of the considered ion, while M m and M γ are mass scales which are inversely proportional to the coupling strengths of the chameleon to matter and to electromagnetic energy density, respectively. Indeed, writing φ the chameleon field, ρ the matter energy density, and F µν the Maxwell-Faray tensor, an effective field theory treatment yields the potential V = φ ρ/M m for the chameleon-matter coupling potential, and V = φ F µν F µν /2M m for the chameleon-photon coupling potential. In the following, we will use the case of Higgs relaxion portal for the explicit expressions, and refer to it as NP; however our conclusions are also valid for the cases of gauged B-L symmetry and chameleon models, and can be analogously extended there.

III. QED CORRECTIONS TO THE NP CORRECTIONS
NP contributions to spectroscopic quantities can be expected to be very small, and QED corrections thereto should be even smaller (typically by a factor of α for radiative corrections at the one-loop level). But as was shown in Refs. [47,48], the magnitude of such radiative QED corrections to a contribution from a potential can be comparable to, or even larger than the leading contribution from that potential, when the latter is generated by a highly localised potential. In our case, the leading contribution is the one-electron hypothetical NP correction, shown in Fig. 1. For the ground-state energy of Hlike, Li-like and B-like ions, we have shown that radiative corrections to the NP contribution are much smaller than the leading NP contribution [49]. However, for heavy new bosons, photon-exchange corrections to the NP contribution to the energy levels of Li-like and B-like ions can be comparable or even much larger than the one-electron NP contribution [49]. We anticipate that the same could happen for the g factor of Li-like and B-like ions.

IV. INDIRECT TESTS WITH THE KING APPROACH
We now turn to the detailed study of isotope shifts of the energy levels and g factors of highly charged ions. In Sec. IV A, we first present the standard, simple formulation of the isotope shift in spectroscopic data, and introduce the King representation of isotope shifts. In Sec. IV B, we present in detail the various subleading nuclear corrections to the energy levels and to the g factor, which contribute to the isotope shift, and complicate the analysis of data. Methods to distinguish SM and NP contributions to the isotope shift are discussed in Sec. IV C. Projected bounds are derived in Sec. IV D.

A. Isotope shifts and the King representation
Isotope shift data is a promising avenue [4][5][6][7] to obtain strong bounds on NP parameters. This is easily understood when recalling that several proposed SM extensions would result in new forces between nucleons and electrons. Isotope shift data can thus carry information on potential neutron-electron interactions due to NP. Isotope shift data can conveniently be handled through the King representation [50]. The bulk of the following explanation will be given for the explicit case of the g factor but, as will be further clarified, the same construction can be elaborated for energy levels. Let us consider two levels 1 and 2 of a bound electron in an ion. For both levels, considering two isotopes A and A ′ of the same ion, we write [50] the isotope shift in the g factor as Within the SM, the leading-order (LO) contributions to the isotope shift are The first summand on the r.h.s. of (10) is the leadingorder contribution to the mass shift, featuring the reduced mass The second summand is the leading-order contribution to the field shift, in the nuclear squared charge radii between the two isotopes. The reader should note that both summands are expressed as the product of an electronic and a nuclear factor. Indeed, the coefficients K i and F i only depend on the electronic level considered, while the quantities µ AA ′ and δR 2 AA ′ are purely nuclear properties. As was done in Ref. [7], we introduce, for any pair of isotopes, the notation  10), it can easily be seen that This relation between isotope shifts for two different levels is hence linear within the LO treatment, and the offset does not depend on the specific isotope pair considered. Hence, King plots will always be linear at the LO: points corresponding to any pair of isotopes will fall on the same line. Hence, deviations from King linearity in experimental data can either be explained by smaller, subleading nuclear contributions to the g factor, or by possible NP contributions. We will explore this in detail in the following. But first, we note that, in some cases, the relativistic correction to the field shift can be absorbed within this leading-order, linear framework. Indeed, in a relativistic treatment, the second summand on the r.h.s. of Eq. (10) can be rewritten as Here A ′ , and RLO stands for 'relativistic leading order'. Obviously, the second summand on the r.h.s. of Eq. (12) is no longer the product of an electronic factor with a nuclear one. Nevertheless, if one considers only electronic levels that share the same |κ i |, as we will do in this work, then γ i is in effect fixed, and, substituting H i/j for F i/j , the linear relation (11) still holds between G AA ′ i(RLO) and G AA ′ j(RLO) . As we will see in the upcoming Sec. IV B, however, King linearity is broken by various subleading contributions to the isotope shift. Before we turn to these subleading corrections, we give, for reference, the expressions for the coefficients K i and H i for the levels 1s, 2s and 2p 1/2 which we consider in this work.
The mass shift coefficients are These three results are derived, respectively, in Refs. [51][52][53]. Here E 2s refers to the Dirac bound energy of the 2s level. The functions P (1s) , P (2s) are tabulated in Refs. [51,52], and the function F is tabulated in Ref. [53] for Z > 20.
The field shift coefficients are where λ e ≡ /m e c is the Compton wavelength of the electron. The expressions for the functions f 1s , f 2s , f 2p 1/2 are given in Ref. [54].
As was mentioned earlier, the foregoing considerations are also valid for the (dimensionless) energy levels E i / m e c 2 of bound electrons and their isotope shifts. In particular, for the ground-state energy of H-like ions, the coefficients are with P 0 tabulated in Ref. [55], and (16) Again, in this case, the relativistic correction to the field shift can be included directly, without breaking King linearity, as long as all the levels considered share the same γ i . Note that (also see Ref. [42]) which inspired the introduction in Ref. [19] of the reduced g factor of H-like ions, in which the leading finite nuclear size correction cancels to a great degree. As we will see, for several subleading nuclear corrections, the same relation obtains at least approximately, with the same proportionality factor x. This motivates the introduction of a King representation where the isotope shifts of the energy level and g factor of the ground state of H-like ions are expressed as a function of each other.

B. Subleading nuclear corrections and Standard-Model King nonplanarities
Let us introduce an extra term in the isotope shift. As was done in the previous Sec. IV A, we write the general abstract expressions for the g factor, but an identical construction can be, and indeed is, made for the energy levels. The isotope shift now reads Here SM stands for 'Standard Model', as s i AA ′ encompasses all the SM other contributions to the isotope shift, besides the leading-order mass shift, and the relativistic leading-order field shift. The third summand s i AA ′ on the r.h.s. is thus generated by subleading nuclear corrections to the g factor. We will examine these various corrections in detail in the following. For now, we can handle this term in a more abstract way, and start by noting that it is not the product of an electronic and a nuclear factor. Hence, the relation between isotope shifts for two electronic levels is There is now an extra offset, in the relation between the isotope shifts of two electronic levels. This offset depends on the specific isotope pair considered, which means that the King plot will no longer be linear. Typically, though, the subleading terms s i AA ′ are much smaller than the other two summands on the r.h.s. of (19), so that the deviations from King linearity due to these terms are only detected on isotope shift data if that data is sufficiently precise.
We now turn to the specific determination of the term s i AA ′ . Our analysis includes the seven largest subleading nuclear corrections to the energy level and the g factor: 1. the higher-order nuclear recoil and 2. nuclear size corrections, 3. the mixed finite-size recoil correction, 4. the radiative recoil and 5. radiative nuclear size corrections, as well as 6. the nuclear deformation and 7. nuclear polarization corrections. For reasons that will become clear, we will devote specific attention to the theoretical uncertainties on these various contributions. We will focus, from this point onwards, on the g factor and the (dimensionless) ground-state energy level of H-like ions. A representation of the contributions from all seven considered nuclear corrections, to a specific isotope shift in H-like calcium, is given in Fig. 2. In the case of ions with more than one electron, nuclear recoil and finite nuclear size corrections to the inter-electronic interaction would also contribute to King nonlinearities, but we do not consider this case in further detail here.

Higher-order nuclear recoil correction
The leading-order nuclear recoil correction expressed by Eqs. (13) and (15) gives a contribution to the boundelectron g factor proportional to m e /M A , where m e is the electron mass, and M A the nuclear mass (for a specific isotope). A result valid to all orders in the electron-nucleus mass ratio can be obtained, which does not take account of radiative or electron-electron interaction contributions. The result is given by Eqs. (29)-(30) of Ref. [56] (also see Ref. [57]); one should, to obtain the higher-order nuclear recoil correction, substract from it the leading-order nuclear recoil term (m e /M A ) (Zα) 2 /n 2 . The result of Ref. [56] is valid, as we said, at all orders in m e /M A , but is given as an expansion in (Zα). Inspection indicates that the next, unacccounted-for terms should feature an extra (Zα) 2 , which is how we estimate the uncertainty on this correction.
For the energy level, the higher-order nuclear recoil correction is obtained from Eqs. (48) and (49) of Ref. [58]. The unaccounted term should be of order (Zα) 6 (m e /M A ) 2 , which yields the estimate of the uncertainty on this correction.

Higher-order finite nuclear size correction
The leading-order finite-nuclear size correction expressed by Eqs. (14) and (16) should be supplemented by smaller, subleading contributions. Calculation of these contributions can be achieved by numerically calculating the full finite-nuclear size correction to the g factor and energy level, e.g. following the method of Ref. [32], and subtracting the already-included leading-order contribution. The uncertainty on these contributions is dominated by that on the nuclear radius.

Finite-size recoil correction
As mentioned in Ref. [21], the nonrelativistic approximation of the finite-size recoil correction to the energy levels can be obtained through the substitution H Hence, following Ref. [42], we have, in the nonrelativistic approximation, H i → H i 1 − 4 me MA for the finite-size recoil correction to the g factors of the levels under study here. The uncertainties in this estimate are given by multiplying this nonrelativistic estimate by (Zα) 2 , however, finite-size recoil corrections to the H-like energy level were computed with three significant digits in Ref. [55], so that we can retain this latter level of uncertainty as very realistic.

Radiative corrections to the nuclear recoil correction
Radiative (QED) corrections to the nuclear recoil contribution are mixed QED+nuclear corrections. For the g factor of H-like ions, they are given by [59] ∆g rad. rec.
(21) Similarly, the radiative recoil correction to the energy level is given by [58] ∆E rad. rec. 1s Clearly, when considering the isotope shift, the first summand on the r.h.s. of Eqs. (21) and (22), linear in m e /M A , will bring a small correction to the coefficient K i , and hence no King nonlinearity. The terms quadratic in (m e /M A ) 2 are extremely small, but taken into account in our treatment. Their uncertainties are estimated by multiplying the quadratic mass ratio term by (Zα) 2 in the case of the g factor (with Eq. (13) in mind), and by (Zα) in the case of the energy level (with Eq. (15) in mind).

Radiative corrections to the finite size
The most detailed study of radiative nuclear size corrections to the g factor and energy of H-like ions are found in Refs. [28,60], respectively. The correction to the g factor is cast as (23) The correction to the energy level can be recast as (24) The coefficients where G NS and G NSQED , as well as D NS and D NSQED are functions of Z and of the nuclear radius R, and are numerically calculated coefficients, with G NS ≃ 1 (within 20%) and −2 < G NSQED < −0.4 for Z ≤ 92. The weighted difference only cancels the radiative nuclear size correction quite weakly (less than onedigit cancellation). When considering the isotope shift, it should be kept in mind that all of R, G NS and G NSQED are isotope-specific. In principle, the uncertainties on these radiative nuclear size corrections is only limited by the knowledge of the nuclear radii. This indicates that the relative uncertainties on the radiative nuclear size corrections can be as small as roughly one part pert thousand.
Finally, the two-loop radiative corrections to the finite nuclear size, and to the nuclear recoil contributions, are extremely small, of course, and can be ignored altogether.

Nuclear deformation correction
For nonspherical nuclei, which represent most of the nuclei, there exists a nuclear deformation correction (also called nuclear shape correction) to the g factor. Typically, it is considered that nuclei possess a quadrupole and a hexadecapole deformation, although they may also exhibit an octupole deformation [61]. A nonperturbative calculation of the nuclear deformation correction to the g factor has also been carried out, and presented in Ref. [31]. In this approach, the wave functions for the bound electrons are obtained by numerically solving the Dirac equation for the potential generated by a deformed nucleus. It is found that, taking into account nuclear deformation at all orders modifies the results obtained through the perturbative method of Ref. [30], to a relatively large extent. Nevertheless, the nuclear shape corrections to the g factor and energy obey ∆g NS 1s = x∆E NS 1s / m e c 2 to a very good level of approximation, leading to a cancellation of about two digits in the reduced g factor, and in the uncertainties, which is very favourable for our purposes.

Nuclear polarization correction
The bound electron and the nucleus can exchange photons, which excite virtual nuclear transitions from and to the nuclear ground state [29,62], giving rise to the nuclear polarization correction. The detailed expressions for the corresponding correction to the boundelectron g factor can be found in Refs. [19,62]. It has been shown through explicit numerical calculations that ∆g NPOL 1s ≃ x∆E NPOL 1s / m e c 2 to a good level of approximation [19], leading to a cancellation of at least one digit in the reduced g factor difference.

C. A modified King representation
As we just discussed in detail, small nuclear corrections to the g factor and the energy levels cause violations of King planarity within the SM. This limits the range of the test of NP by inspection of King plot data used in Ref. [5]. In Ref. [7], we used the method of Ref. [5] to study tests of NP with the bound-electron g factor, while also making sure that the sensitivity to NP at the projected level of precision for potential Penning trap experiments would not be confounded by these subleading nuclear corrections. In the present work, we propose an extension of this method, which circumvents its fundamental limitation: namely, we propose to test the King nonlinearity of experimental data against its predicted theoretical value within the SM, instead of testing it against the zero nonlinearity predicted by the leading-order SM contributions to the isotope shift (see Sec. IV A). In all approaches, we need to consider the hypothetical NP contribution to King nonlinearity.
In the presence of NP, the isotope shift reads where the fourth summand on the r.h.s. corresponds to the hypothetical NP contribution to the g factor. The same general expression holds for the (dimensionless) energy levels. The minimal number of isotopes to consider to study King plots is four. As such, it is interesting to retain the general formalism developed in Ref. [5], even though it is not directly applicable to data sets with more isotopes. We note that, in any case, little high-precision data on the isotope shift of g factors is available, with the notable exception of Ref. [63], which motivates the use of this method, adapted to the case with the smallest number of isotope shift data points, that is sufficient to assess King planarity. We recall the definition We now consider four isotopes A, A ′ 1 , A ′ 2 and A ′ 3 and define the vectors Note, for future purposes, that for the Higgs portal, In the rest of this work, we thus give expressions for the Higgs portal, to be explicit, but these expressions are straightforwardly applied to the case of the unbroken B−L symmetry when substituting the coupling constant 4πg 2 B−L for α HP , and, as can be seen from Eq. (8), by substituting m e m p/n /4πM 2 m . King nonplanarity (or nonlinearity) is measured [5] by the parameter which can be calculated to be equal to Note that if n i AA ′ is expressed as the product of A − A ′ with a purely electronic factor, as is the case for the Higgs portal and the gauged B − L symmetry, then N 1 × N 2 = 0. If we neglect the subleading nuclear contributions to the isotope shift (S i = 0), and consider that there are no NP contributions (N i = 0), then it is seen from Eq. (29) that King nonplanarity will vanish. Hence, nonplanar King data might be interpreted as a sign of NP, provided that this nonplanarity is not expected to arise from the subleading nuclear contributions. Experimentally, the presence or absence of nonplanarity can only be assessed at a certain level of accuracy, since, due to experimental uncertainties, we could not expect King data to verify N = 0 exactly, even in a Gedankenexperiment in which S i = N i = 0. We propose here to introduce a modified King nonplanarity, that accounts for (and hence, in a way, sets aside) the deviation from the nonplanarity expected within the SM. This nonplanarity is obtained by a combination of experimental data and theoretical calculations of subleading nuclear corrections, and reads: The first summand can be obtained experimentally, by collecting isotope shift data. The second summand can be calculated theoretically with the methods summarised in Sec. IV B. The new nonplanarity test is performed as follows. The modified nonplanarity parameter (30) is first compared to its first-order propagated error σ N ′ . If |N ′ | < σ N ′ , the data is considered planar. Then, the NP parameter which is to be constrained, is considered to be bound by its own first-order propagated error, which can be computed from Eq. (29). Note that the error σ N ′ features not only experimental (as in Refs. [5,7,39]) but also theoretical contributions. Indeed, the subleading nuclear contributions S i to the isotope shift have uncertainties. The explicit expression of the error is: where G i(j) (resp. S i(j) ) is the j component of G i (resp. of S i ), which is an experimental (resp. theoretical) quantity. This modified King planarity test can detect nonplanarities caused by NP contributions to the g factor, even when these nonplanarities are much smaller than those caused by the subleading nuclear corrections to the g factor. This can be seen as follows. Let us write E i the (signed) difference between the measured isotope shifts G i and the corresponding SM predictions, obtained through Eq. (19). In the very high precision regime of experiments, where E i is expected to be smaller than the subleading nuclear corrections S i , we can approximate σ N ′ by its second summand in Eq. (31), that is, by where δ (i) is the i component of δ; andī = 3 − i. This first-order propagated error is to be compared to where E i(j) is the j component of E i and ǫ (i) = (−1) 1+i . If there are no NP contributions, the deviations E i between the experimental results and the SM predictions can be expected to become smaller with improving experimental precision. In this case, |N ′ | decreases too, and we can expect to reliably have |N ′ | < σ N ′ , which allows one to set bounds on the NP parameters. If the theory of the nuclear corrections to the g factor improves, then the uncertainties S i(j) will decrease, causing a similar decrease in σ N ′ . Nevertheless, if there are no NP contributions to the g factor, E i will also decrease with improving theoretical precision, as the theory will bet-ter match the experiments. This means that, within this programme, there is in principle no limit of applicability of the King planarity test, contrasting with the test used in Refs. [5,7].

D. Projected bounds
To derive bounds, we recast the (standard) King nonplanarity as [7] In the following, quantities with index i = 1 will refer to the H-like g factor, and those with index i = 2 will refer to the dimensionless H-like ground-state energy. Our King representation will be explicitly constructed, to compare the isotope shifts of these two quantities. This expression can be inverted to obtain Now, note that while studying the various subleading nuclear corrections to the g factor and the energy levels, we have mentioned that, for several of them, the equality ∆g 1s ≃ x∆E 1s / m e c 2 ≃ 0 holds in good approximation. Hence, we can write with T i coming from the subleading nuclear corrections for which ∆g 1s ≃ x∆E 1s / m e c 2 ≃ 0, and U i from the other corrections. Fortunately, H 2 /H 1 X 1 and X 2 , which appears in the denominator on the r.h.s. of Eq. (35), does not cancel to any appreciable degree except in the very large boson mass limit (as can be checked by Eqs. (3), (6), (14) and (16)). The bounds on α HP are given by the first-order propagated error where the derivatives can be obtained from Eq. (35). We repeat that here, G 1 and G 2 are to be understood as purely experimental quantities, and so too is N , which is given in terms of these two vectors by Eq. (28). The ∆G i(j) are hence experimental uncertainties. On the other hand, S 1 and S 2 are to be calculated, and the ∆S i(j) are theoretical uncertainties. They are the uncertainties of the subleading nuclear corrections. It is important to remember that H 1 = x H 2 . As can then be seen by inspecting the r.h.s. of Eq. (35), this means that, for all occurrences of S i in Eqs. (35) and (37), the corresponding contribution will be dominated by U i , the contribution from the subleading corrections that do not cancel appreciably through the reduced g factor. On Fig. 3, we show the contributions to the predicted SM King nonplanarity in the vanishing boson mass limit from all seven considered subleading nuclear corrections. We also show the contributions, from these same corrections, to the bound on the NP coupling constant. It is seen that the largest contributors to the theoretical uncertainty budget are the higher-order nuclear recoil and nuclear size contributions and, to a lesser extent, the nuclear polarization correction. The latter is currently being investigated with unprecedentedly accurate nuclear models [32,65], and it is anticipated that the calculation uncertainties will be significantly reduced in the near future. It is more difficult to substantially improve the calculation of the higher-order recoil correction.
For that reason, we introduce the square-mass King representation. Let us rewrite the isotope shift (25) as Here, we have written the quadratic recoil contribution explicitly. We now introduce the square-mass normalised quantities Now, with vectors constructed as in Eq. (26), we can proceed to similar manipulations, to obtain the square-mass King nonplanarity Of course, this expression is not equal to that given in Eq. (34), but if one compares the terms one by one, it is seen that the first two summands on the r.h.s. of Eq. (40) have an equivalent in Eq. (34), while the third one does not. That third term is the contribution to N β from the leading-order mass shift. In the absence of NP, and if one neglects the subleading nuclear corrections, this third summand survives as a contribution to the square-mass King nonplanarity (39), which should hence not be expected to vanish when NP is absent and the subleading nuclear corrections are neglected. We can then use Eq. (40) to solve for α HP , and set bounds on the NP coupling constant through The uncertainty on the NP coupling constant, due to the leading nuclear recoil correction, is reduced by four orders of magnitude compared to that generated by the higher-order nuclear recoil correction. The contribution from the six other subleading nuclear corrections remain essentially unchanged from those shown in Fig. 3, so that this new square-mass King representation is a clear improvement.

V. DIRECT TESTS
Despite all the refinements introduced in Sec. IV, searching for NP through King planarity analysis turns out not to be the most efficient approach for few-electron ions. The King approach is needed for many-electron systems, in particular because for these systems, the electron-interaction contributions to the nuclear recoil, the so-called specific mass shift, cannot be computed accurately. The King approach, in which the nuclear recoil is eliminated, is thus useful for such systems. With few-electron ions, on the other hand, calculations can be carried out more comprehensively and more accurately, which indicates the interest of a direct approach to the search for NP.
In this approach, we use the small discrepancy between theory and experiment to derive bounds on hypothetical New Physics (NP) contributions. In some cases, theory and experiment may agree within their respective error bars, while, in other cases, the error bars may not overlap. However, which case is irrelevant to our specific concern here, since the maximum discrepancy between theory and experiment allowed by the error bars may be set as the maximum hypothetical NP contribution in either case (see Fig. 4). This approach was explored further in Ref. [7] in the case of the g factor of few-electron ions, but, implemented on a single g factor, it necessitates challenging improvements in the calculation of QED corrections, and also interelectronic-interaction contributions, to spectroscopic quantities. Many such corrections are under active consideration [23][24][25][26]66].
A better approach consists of the same direct comparison between experiment and theory, but implemented on a single isotope shift [67]. Thus, radiative contributions from QED largely drop out, reducing the amount of necessary theory input. As is clear from the foregoing discussion, the bound that can be set on NP is limited by the uncertainty δs AA ′ on the subleading nuclear corrections. In the Higgs portal scenario, the NP contribution to the isotope shift of the g factor or energy level is n AA ′ = α HP (A − A ′ ) X. The difference in the number of neutrons will typically be of the order of unity, and as can be seen from Eqs. (3) and (6), the electronic coefficient is roughly Zα in the light boson limit, so that for Z > 10, a given uncertainty on the subleading nuclear corrections will yield a 10 times larger contribution, very roughly speaking, to the bound on α HP . This is to be contrasted with the corresponding contribution in the King approach, which, as can be seen in Figs. 2 and 3, is several orders of magnitude larger than the corresponding uncertainties on the subleading nuclear corrections.

VI. CONCLUSION
In this work, the search for New Physics with fewelectron ions through King linearity analysis has been pushed to what is arguably its farthest extent. We have modified the King formalism to explicitly account for subleading Standard Model nuclear corrections, which contribute to isotope shifts and to King nonlinearity. This yielded a King linearity violation test which can be carried out in the very-high-precision experimental regime. Indeed, this test enables detection of New Physics contributions to spectroscopic quantities which would be smaller than the subleading nuclear corrections. We also introduced the square-mass King representation to successfully suppress uncertainties from the higherorder recoil. Even with all these refinements, the bounds which can be set with this approach remain less com- to the g factor from one of the New Physics candidates (Higgs portal). This directly allows the setting of bounds on the New Physics parameters.
petitive than bounds obtained with the spectroscopy of many-electron ions [5,8,9] putting forward a simpler approach with a single isotope pair and no King analysis can yield more stringent bounds with the same amount of theory input [67]. Finally, our results stress the importance of a detailed analysis of the subleading corrections to the King's plot, and would hopefully stimulate for same investigation for the other systems [4-6, 8, 9] being used for the search of the New Physics beyond the Standard Model.