Effect of correlated hyperfine theory errors in the determination of rotational and vibrational transition frequencies in HD⁺

Abstract

We revisit the problem of separating the purely rovibrational and hyperfine contributions to observed rovibrational transition frequencies. The separation of hyperfine frequency shifts played an essential role in recent works in which the value of the proton-electron mass ratio was determined from Doppler-free rovibrational spectroscopy and high-precision theory of the ${\mathrm{HD}}^{+}$ molecular ion. Here, we present an alternative approach to removing theoretically computed hyperfine shifts from a limited set of observed rovibrational transition frequencies, based on a least-squares adjustment. The adjustment facilitates a quantitative study of the effect of correlations between the uncertainties of lower- and upper-state hyperfine energies on the derived value of the rovibrational transition frequency. We find that correlations may lead to significant shifts of the extracted rovibrational transition frequencies. In addition, we propose and demonstrate a method to quantify the additional shift and uncertainty, which is of importance to other applications of the spectroscopic data, such as the determination of fundamental physical constants.

GRAPHICAL ABSTRACT

Keywords:

• Hydrogen molecular ions
• rotational-vibrational spectroscopy
• molecular hyperfine structure
• least-squares adjustment of physical constants

1. Introduction

Since the development of wave mechanics nearly one century ago, the three-body hydrogen molecular ion (HMI) has been a benchmark system for quantum molecular structure calculations (for a historical overview, see [1]). A number of milestones in the development of the HMI molecular theory – in particular concerning methods for performing nonadiabatic calculations – were achieved by Kołos and Wolniewicz [2], and later by Wolniewicz and Poll [3–5]. With the advent of fully nonadiabatic three-body nonrelativistic calculations around the turn of the millennium, and the subsequent evaluation of high-order relativistic and quantum electrodynamics (QED) corrections, the state of the art has now advanced to the point where transition frequencies of rotational and vibrational transitions in the HMI electronic ground state can be computed with a theoretical uncertainty of the order of 10 parts per trillion (ppt) [6,7]. This level of precision is remarkable because for the first time, the uncertainty of theoretical transition frequencies is no longer limited by the uncertainty contribution from quantum molecular theory, but rather by the uncertainty introduced by the value of the proton-electron mass ratio, $m_p/m_e$, published by the Committee on Data of the International Science Council (CODATA) [8]. In 1976, Wing et al. had already suggested that under such conditions, it may be possible to derive an improved value of $m_p/m_e$ provided that sufficiently accurate experimental transition frequencies are available [9], a suggestion that was more recently extended by Karr et al. to include other physical constants as well [10].

In 2020, measurements performed on the deuterated HMI, ${\mathrm{HD}}^{+}$, were reported, including frequency values of a vibrational overtone [11] and of a rotational transition [12], with relative uncertainties of 3 ppt and 17 ppt, respectively. Both experiments made use of ion traps and sympathetic cooling and detection of ${\mathrm{HD}}^{+}$ via laser-cooled beryllium ions, in combination with techniques for Doppler-free spectroscopy. For the $(v,N)$: $(0,3)\to (9,3)$ vibrational overtone (where v and N stand for the vibrational and rotational quantum numbers, respectively), Doppler-free laser excitation in the Lamb-Dicke regime was achieved using counter-propagating lasers with similar wavelengths [11,13]. For the $(v,N)$: $(0,0)\to (0,1)$ transition, it was the relatively long wavelength of the terahertz radiation used for rotational spectroscopy that enabled Doppler-free excitation. The results of each of the experiments were combined with theoretical calculations to derive a value of $m_p/m_e$ with an unprecedented relative uncertainty of $2\times {10}^{-11}$. More recently, a measurement of the $(v,N)$: $(0,0)\to (1,1)$ rovibrational transition in ${\mathrm{HD}}^{+}$ was reported, also with sufficient accuracy for a meaningful derivation of $m_p/m_e$ [14]. The three high-resolution measurements [11,12,14] were later combined with other spectroscopic data from three-body systems by Germann et al. to set bounds on the validity of QED and possible manifestations of physics beyond the standard model [15]. In this particular analysis, Germann et al. used values from independent determinations of $m_p/m_e$, and of the deuteron-proton mass ratio, $m_d/m_p$, to evaluate the theory [16–20].

Whether it is the goal to test QED or to determine $m_p/m_e$, rotational or rovibrational transition frequencies need to be measured. In practice, however, experimentally accessible transition frequencies also contain large shifts (of the order of 100 MHz) due to hyperfine structure in the lower and upper rovibrational levels. The measured frequencies therefore need to be corrected for these hyperfine shifts, in order to obtain the spin-averaged rovibrational transition frequency. In the absence of experimental data on the hyperfine structure for the transitions considered here, various approaches have been developed that all make use of theoretical ab initio predictions of the ${\mathrm{HD}}^{+}$ hyperfine structure. The theoretical hyperfine structure has been computed with an uncertainty of the order of 1 kHz [21,22], which is comparable to the experimental uncertainty achieved in the $(v,N)$: $(0,3)\to (9,3)$ vibrational spectroscopy [11]. However, rotational spectroscopy of hyperfine components of the $(v,N)$: $(0,0)\to (0,1)$ transition has been carried out with a much smaller uncertainty of 17 Hz [12]. Correcting these transitions using the theoretical hyperfine shifts would reduce the accuracy by nearly two orders of magnitude, to about 1 kHz. This issue is circumvented through an approach that exploits the tracelessness of the effective spin interaction Hamiltonian that is used to compute the energies of hyperfine levels: by measuring a suitable selection of hyperfine components, and constructing a composite transition frequency by taking an appropriately weighted sum of the measured frequencies, the uncertainty of the hyperfine structure can be made to cancel out to a large extent [12,23]. It was pointed out by Alighanbari et al. that in this process, one has to make assumptions about (possibly) correlated uncertainties of the lower- and upper-level hyperfine energies [12]. The experimental frequency value of the $(v,N)$: $(0,0)\to (0,1)$ spin-averaged transition was eventually derived for one specific correlation scenario, in which such correlations were assumed to be absent. However, the possible effect (in terms of frequency shift and frequency ucertainty) of other conceivable correlation scenarios on this result was not quantitatively discussed [12].

Correlations between theory errors have been considered for atomic hydrogen and deuterium, which have traditionally played a dominant role in CODATA adjustment of the Rydberg constant, and proton and deuteron electric charge radii. For these atomic systems, the theory permits identifying state-independent and state-dependent uncertainties for the various uncalculated terms, resulting into unambiguous correlation coefficients [8]. Such an approach could to some extent be employed for the HMIs as well. For rovibrational states of ${\mathrm{HD}}^{+}$, however, the nature of higher-order terms can differ substantially for lower and upper vibrational states [24]; see also Section 2.2 below. In addition, both recent ${\mathrm{HD}}^{+}$ measurements revealed deviations between the theoretically predicted and experimentally determined hyperfine structure [11,12,22]. Such deviations might point to a missed term in the theory, the nature of which is unknown, such that correlation coefficients cannot be reliably estimated.

Here we introduce a different approach, based on a least-squares adjustment of the spin-averaged transition frequency to an arbitrary selection of experimentally observed hyperfine components and the corresponding theoretical expressions. The above-mentioned correlations can be taken into account through the covariance matrix of the input data. The least-squares adjustment automatically takes advantage of the tracelessness of the effective spin Hamiltonian (which enters the adjustment through the theoretical expressions for the hyperfine shifts), and produces results with the same precision as the composite frequency method put forward in Refs. [12,23]. However, our analysis also reveals that, in general, the adjusted spin-averaged frequency may vary significantly depending on the assumed correlation coefficients. As such it lays bare a hitherto unrecognised systematic effect (which is of relevance to the both the results published by Alighanbari et al. and Patra et al.) that should be taken into consideration when determining spin-averaged transition frequencies, and derived quantities such as the value and uncertainty of $m_p/m_e$ from ${\mathrm{HD}}^{+}$. We also find that the ${\chi }^2$ value produced by the least-squares adjustment can be helpful in establishing which correlation scenarios are consistent with the observed data, and which other correlation scenarios are less likely. To obtain a quantitative measure of the newly identified 'correlation shift' and the additional uncertainty associated with it, we carry out a Monte-Carlo simulation that takes advantage of relative weights that can be derived from the ${\chi }^2$ values.

This article is organised as follows. The least-squares adjustment used to separate rovibrational and hyperfine contributions is introduced in Section 2, along with the theoretical model of the ${\mathrm{HD}}^{+}$ hyperfine structure. In the same chapter, the covariance between different hyperfine components is derived, and the various types of correlation involved are discussed. In Section 3, the results of the adjustment (considering various correlation scenarios) and of the Monte-Carlo runs are given for both the $(v,N)$: $(0,3)\to (9,3)$ and $(v,N)$: $(0,0)\to (0,1)$ transitions described in Refs. [11,12]. The impact of correlations on the spin-averaged transition frequencies and the interpretation of the Monte-Carlo results are discussed in Section 4, followed by concluding remarks in Section 5.

2. Separating hyperfine shifts from rovibrational transition frequencies

Let $f_{\mathrm{ex}}^{ij}$ denote an experimentally determined frequency of a given rotational or rovibrational transition within the $1s\sigma$ electronic ground state of ${\mathrm{HD}}^{+}$. The superscripts i and j are used to label the hyperfine state of the lower and upper level, respectively. Below it will become clear that these superscripts actually represent a set of five quantum numbers that uniquely identify each state. Furthermore, we will loosely use the term 'rovibrational' to indicate vibrational as well as rotational transitions. From molecular theory it is known that $f_{\mathrm{ex}}^{ij}$ can be written as the sum of the pure rovibrational transition frequency (also known as the spin-averaged rovibrational transition frequency), $f_{\mathrm{sa}}$, and the hyperfine frequency shift, $f_{\mathrm{hf}}^{ij}$. This separation is useful for various reasons. First, it can be used to determine an 'experimental' spin-averaged transition frequency, $f_{\mathrm{ex},\mathrm{sa}}$, which may be readily compared with the theoretical prediction of the same quantity, $f_{\mathrm{th},\mathrm{sa}}$, in order to test the theory. Second, $f_{\mathrm{th},\mathrm{sa}}=f_{\mathrm{th},\mathrm{sa}}(m_p/m_e)$ is a function of $m_p/m_e$ (as well as other physical constants; see Section 4). Hence, $m_p/m_e$ can be determined from the comparison between theory and experiment by solving the equation $f_{\mathrm{ex},\mathrm{sa}}=f_{\mathrm{th},\mathrm{sa}}(m_p/m_e)$.

To determine $f_{\mathrm{ex},\mathrm{sa}}$ from a single measured hyperfine component of a rovibrational transition, we can solve the equation (1) $f_{\mathrm{ex}}^{ij}=f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{hf}}^{ij},$(1) provided that the hyperfine shift $f_{\mathrm{hf}}^{ij}$ is known. For the rovibrational states of ${\mathrm{HD}}^{+}$ considered here, however, no independent experimental data exists on the hyperfine structure. Instead, we can use the theoretical hyperfine structure, which has been calculated with an uncertainty of the order of 1 kHz, or ${10}^{-6}$ relative to the hyperfine splitting ($\sim 1$ GHz) [21,22].

In each of the recent experiments [11,12], a number $M\ge 2$ of hyperfine components were measured. Identifying each component by the pair of quantum labels ij, we can form a set of M equations of the form (2) $f_{\mathrm{ex}}^{ij}\doteq f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{th},\mathrm{hf}}^{ij}.$(2) The dotted equal sign indicates that the two sides of each equation are in general found to be only approximately equal, because the set of equations is overdetermined. If we again use the hyperfine shifts $f_{\mathrm{th},\mathrm{hf}}^{ij}$ from theory, then the value of $f_{\mathrm{ex},\mathrm{sa}}$ can be found by a least-squares adjustment based on Equation (2(2) $f_{\mathrm{ex}}^{ij}\doteq f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{th},\mathrm{hf}}^{ij}.$(2) ). In what follows, we will use the terminology of the CODATA least-squares adjustments [25], and refer to Equation (2(2) $f_{\mathrm{ex}}^{ij}\doteq f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{th},\mathrm{hf}}^{ij}.$(2) ) as 'observational equations'.

2.1. Hyperfine structure of ${\mathrm{HD}}^{+}$

In ${\mathrm{HD}}^{+}$, the hyperfine structure arises from magnetic interactions between the proton, deuteron, electron, and molecular rotation, as well as the interaction of the deuteron electric quadrupole moment with the molecular internal electric field gradient. These interactions are described (to lowest order) by the Breit-Pauli Hamiltonian, but in practice it is more convenient to capture all interactions within an effective spin Hamiltonian that is specific to a given rovibrational state v, N. The currently most elaborate version of the effective spin Hamiltonian contains nine interaction terms [26]. The strength of each of the nine terms is characterised by a spin coefficient, $E_k^{vN}$ ($k=1,\dots ,9$), which are obtained numerically by averaging the corresponding terms of the Breit-Pauli Hamiltonian, as well as possible high-order terms derived using the nonrelativistic quantum electrodynamics (NRQED) approach, over the rovibrational wavefunction [21,22,26].

To find the eigenvalues of the effective spin Hamiltonian, we choose a basis set of pure spin states $|G_1,G,F⟩$. The quantum numbers $G_1$, G, and F follow from the angular momentum coupling scheme ${\mathbf{G}}_1={\mathbf{I}}_p+{\mathbf{s}}_e$, $\mathbf{G}={\mathbf{G}}_1+{\mathbf{I}}_d$, $\mathbf{F}=\mathbf{G}+\mathbf{N}$. Here, ${\mathbf{I}}_p$, ${\mathbf{I}}_d$, and ${\mathbf{s}}_e$ stand for the proton, deuteron and electron spin, respectively, with corresponding quantum numbers $I_p=1/2$, $I_d=1$, and $s_e=1/2$. The theoretical hyperfine level energies, $W_{\mathrm{hf}}^i$, are obtained by diagonalisation of the effective spin Hamiltonian matrix. The diagonalisation introduces some hyperfine state mixing, so that the quantum numbers $G_1$, G are only approximately good.

Obviously, the theoretical hyperfine energies are functions of the coefficients $E_k^{vN}$, and we express them (in frequency units) as $W_{\mathrm{hf}}^i\equiv W_{\mathrm{hf}}^{vN{G}_{1}GF}(E_1^{vN},\dots ,E_9^{vN})$. Using double and single primes to indicate quantum numbers of lower and upper states, respectively, Equation (2(2) $f_{\mathrm{ex}}^{ij}\doteq f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{th},\mathrm{hf}}^{ij}.$(2) ) then takes the form (3) $\begin{array}{rl}& f_{\mathrm{ex}}^{v^{″}{v}^{\prime }{N}^{″}{N}^{\prime }{G}_1^{″}{G}_1^{\prime }{G}^{″}{G}^{\prime }{F}^{″}{F}^{\prime }}\\ & \doteq f_{\mathrm{ex},\mathrm{sa}}^{v^{″}{v}^{\prime }{N}^{″}{N}^{\prime }}+W_{\mathrm{hf}}^{v^{\prime }{N}^{\prime }{G}_1^{\prime }{G}^{\prime }{F}^{\prime }}(E_1^{v^{\prime }{N}^{\prime }},\dots ,E_9^{v^{\prime }{N}^{\prime }})\\ & -W_{\mathrm{hf}}^{v^{″}{N}^{″}{G}_1^{″}{G}^{″}{F}^{″}}(E_1^{v^{″}{N}^{″}},\dots ,E_9^{v^{″}{N}^{″}}).\end{array}$(3) To simplify notation, we will use ρ to denote the rovibrational quantum numbers v, N, and β to indicate the hyperfine quantum numbers $G_1,G,F$ of a given rovibrational-hyperfine state. With this notation, the superscripts i, j in Equation (1(1) $f_{\mathrm{ex}}^{ij}=f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{hf}}^{ij},$(1) ) become $i={\rho }^{″}{\beta }^{″}$ and $j={\rho }^{\prime }{\beta }^{\prime }$. The observational equations can thus be written as (4) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}}),\end{array}$(4)

2.2. Hyperfine structure: theoretical uncertainty and correlated errors

We now make the following observations regarding Equation (4(4) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}}),\end{array}$(4) ). First, each spin coefficient has an associated uncertainty, $u(E_k^{\rho })$. This uncertainty reflects the estimated size of the error that is made by various approximations in the theory. Although this error is not a stochastic quantity, we will follow the common practice to treat it as a normally distributed variable with zero mean and standard deviation $u(E_k^{\rho })$. Second, the observational equations are inherently correlated, since each of them depends on the same set of 18 hyperfine spin coefficients. This inherent correlation also provides the mechanism that converts the tracelessness of the effective spin Hamiltonian [23] into a reduced uncertainty of the adjusted value of $f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}$. Third (and crucially), the errors of different spin coefficients may be correlated, and we can distinguish three types of correlations:

• Correlated uncertainties of different spin coefficients in the same (lower or upper) rovibrational state, e.g. $u(E_k^{{\rho }^{″}})$ and $u(E_l^{{\rho }^{″}})$ with $k\ne l$ (type A)

• Correlated uncertainties of pairs of lower-state and upper-state spin coefficients, e.g. $u(E_k^{{\rho }^{″}})$ and $u(E_k^{{\rho }^{\prime }})$ (type B)

• Correlated uncertainties of pairs of different spin coefficients, each in a different (lower or upper) state, e.g. $u(E_k^{{\rho }^{″}})$ and $u(E_l^{{\rho }^{\prime }})$ (type C)

Such correlations may arise because the coefficients are derived using the same underlying theoretical framework (NRQED), and uncertainties may stem from a common source of error. This might be the case for, for example, the coefficients $E_4^{\rho }$ and $E_5^{\rho }$, associated with the proton-electron and deuteron-electron Fermi-contact interaction, respectively [22]. However, if the nature of the errors is different, then the uncertainties might be less correlated. This could occur, for instance, for a pair of uncertainties $u(E_k^{{\rho }^{″}})$ and $u(E_k^{{\rho }^{\prime }})$, where ${\rho }^{\prime }$ corresponds to a highly excited vibrational state. For the excited state, an theoretical error might then occur in the value of $E_k^{{\rho }^{\prime }}$ due to neglected nonadiabatic coupling to electronically excited states, which would be nearly absent for the lower-state coefficient $E_k^{{\rho }^{″}}$ [24]. The fundamental problem here is that the errors and therefore the correlation coefficients are a priori unknown, so that one has to consider different correlation scenarios (i.e. perform the adjustment with different sets of correlation coefficients).

The uncertainties and various correlations are calculated with the help of the known functions $W_{\mathrm{hf}}^{\rho \beta }(E_1^{\rho },\dots ,E_9^{\rho })$, and included in the covariance matrix for the input data. The fractional uncertainties $u(E_k^{\rho })/E_k^{\rho }$ are between ${10}^{-4}$ and ${10}^{-6}$ [11,12], which justifies using linearised expressions for $W_{\mathrm{hf}}^{\rho \beta }(E_1^{\rho },\dots ,E_9^{\rho })$. Let $f_{\mathrm{te}}$ denote the transition frequency that is represented by the theoretical expression, i.e. the right-hand side of Equation (4(4) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}}),\end{array}$(4) ): $\begin{array}{rl}{f}_{\mathrm{te}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \equiv f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}}).\end{array}$Then, for any pair of transitions with theoretical frequency values $f_1\equiv f_{\mathrm{te}}^{{\rho }_1^{″}{\rho }_1^{\prime }{\beta }_1^{″}{\beta }_1^{\prime }}$ and $f_2\equiv f_{\mathrm{te}}^{{\rho }_2^{″}{\rho }_2^{\prime }{\beta }_2^{″}{\beta }_2^{\prime }}$, the covariance matrix elements [25] can be expressed as follows: (5) $\begin{array}{rl}& \mathrm{cov}(f_1,f_2)\\ & =\sum _{k,l=1}^9[\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}r(E_k^{{\rho }^{″}},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{″}})u(E_l^{{\rho }^{″}})\\ & +\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{\prime }}}r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{\prime }})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{\prime }})\\ & +(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}+\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}})\\ & \times r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{″}})\phantom{\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}}].\end{array}$(5) The first two terms in the summation of Equation (5(5) $\begin{array}{rl}& \mathrm{cov}(f_1,f_2)\\ & =\sum _{k,l=1}^9[\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}r(E_k^{{\rho }^{″}},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{″}})u(E_l^{{\rho }^{″}})\\ & +\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{\prime }}}r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{\prime }})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{\prime }})\\ & +(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}+\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}})\\ & \times r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{″}})\phantom{\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}}].\end{array}$(5) ) contain two types of correlations. First, there are the terms with diagonal correlation coefficients (i.e. k = l) for which by definition $r(E_k^{{\rho }^{″}},E_k^{{\rho }^{″}})=r(E_k^{{\rho }^{\prime }},E_k^{{\rho }^{\prime }})=1$. These terms represent the correlation between $f_1$ and $f_2$ due to the shared dependence of the observational equations on the same spin coefficients, and they mediate the effect of the tracelessness of the effective spin Hamiltonian. Second, they contain terms with off-diagonal correlation coefficients ($k\ne l$) within the same rovibrational level (i.e. type-A correlation, for example between $E_4^{\rho }$ and $E_5^{\rho }$). The third term in Equation (5(5) $\begin{array}{rl}& \mathrm{cov}(f_1,f_2)\\ & =\sum _{k,l=1}^9[\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}r(E_k^{{\rho }^{″}},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{″}})u(E_l^{{\rho }^{″}})\\ & +\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{\prime }}}r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{\prime }})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{\prime }})\\ & +(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}+\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}})\\ & \times r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{″}})\phantom{\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}}].\end{array}$(5) ) describes the covariance due to correlations $r(E_k^{{\rho }^{″}},E_l^{{\rho }^{\prime }})$ between similar spin coefficients (k = l, type B) and different ($k\ne l$, type C) spin coefficients in the lower- and upper-state theoretical expressions. As pointed out above, the values of the various correlation coefficients r are a priori unknown (except for the diagonal coefficients), and we will assign values $r\in [-1,1]$ to them depending on the scenario under consideration.

We will only consider correlations for the five largest spin coefficients. These are $E_1^{\rho }$ (for the electron spin-rotation interaction), $E_4^{\rho }$ and $E_5^{\rho }$ (related to the aforementioned Fermi-contact terms), and $E_6^{\rho }$, $E_7^{\rho }$, which represent the tensor part of the proton-electron and deuteron-electron dipolar magnetic interactions, respectively. The spin coefficients for the proton spin-rotation ($E_2^{\rho }$), deuteron spin-rotation ($E_3^{\rho }$), proton-deuteron spin-spin tensor ($E_8^{\rho }$), and deuteron electric quadrupole ($E_9^{\rho }$) interaction have uncertainties that are negligible compared to any other uncertainty involved in the observational equations (including the experimental uncertainty). Therefore, correlation coefficients involving these spin coefficients are set to zero. Furthermore, we will assume that the uncertainties of the various experimental transition frequencies are uncorrelated. It should be noted that for the v = 0, N = 0 state, only $E_4^{00}$ and $E_5^{00}$ have nonzero values [26]. As a consequence, for the $(v,N)$: $(0,0)\to (0,1)$ transition, type-B and type-C correlations only occur for these two spin coefficients.

We here also introduce the correlation coefficients $r_A$, $r_B$, and $r_C$. These are associated with the coefficients r in Equation (5(5) $\begin{array}{rl}& \mathrm{cov}(f_1,f_2)\\ & =\sum _{k,l=1}^9[\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}r(E_k^{{\rho }^{″}},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{″}})u(E_l^{{\rho }^{″}})\\ & +\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{\prime }}}r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{\prime }})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{\prime }})\\ & +(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}+\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}})\\ & \times r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{″}})\phantom{\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}}].\end{array}$(5) ) as follows (grouped by type A,B or C): $\begin{array}{rl}{r}_A:& r(E_k^{{\rho }^{″}},E_l^{{\rho }^{″}}),r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{\prime }})(k\ne l)\\ r_B:& r(E_k^{{\rho }^{″}},E_k^{{\rho }^{\prime }})\\ r_C:& r(E_k^{{\rho }^{″}},E_l^{{\rho }^{\prime }})(k\ne l).\end{array}$Type-A and type-C correlations are only considered for the large spin coefficients which describe similar interactions, i.e. $E_4^{\rho },E_5^{\rho }$ and $E_6^{\rho },E_7^{\rho }$. This implies that we always set (for example) $0=r(E_1^{{\rho }^{″}},E_5^{{\rho }^{″}})=r(E_1^{{\rho }^{″}},E_6^{{\rho }^{″}})=r(E_4^{{\rho }^{″}},E_6^{{\rho }^{″}})$. Note also that the definitions of $r_B$ and $r_C$ above are somewhat restrictive as they imply $r_B=r(E_k^{{\rho }^{″}},E_k^{{\rho }^{\prime }})$ for all k, and likewise for $r_C=r(E_k^{{\rho }^{″}},E_l^{{\rho }^{\prime }})$, i.e. they assign the same value to correlation coefficients which in reality may not be necessarily equal. This restriction does not appreciably affect the results and insights obtained here, as was verified using computations with a greater flexibility in the different correlation coefficients. We nevertheless choose to work with the simplified (yet more restrictive) scenarios represented by the three coefficients $r_A,r_B,r_C$, as it makes the analysis more comprehensible. To aid the reader, a graphical overview of all nonzero correlation coefficients is provided in Figure A1 of Appendix 1.

2.3. Least-squares adjustment

The least-squares adjustment is implemented as done for the CODATA adjustments (see Appendix E of Mohr et al. [25]). Of special interest is the covariance matrix, V, of the input data, through which the various types of correlation are incorporated. Specifically, we will consider scenarios with various amounts of correlation of type A, B, and C, which are denoted by the tuple $(r_A,r_B,r_C)$. We use the theoretical spin-averaged transition frequency, $f_{\mathrm{th},\mathrm{sa}}$, as starting value for the adjusted parameter $f_{\mathrm{ex},\mathrm{sa}}$. Since the observational equations [Equation (2(2) $f_{\mathrm{ex}}^{ij}\doteq f_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{th},\mathrm{hf}}^{ij}.$(2) )] are linear in $f_{\mathrm{ex},\mathrm{sa}}$, a single adjustment suffices and no iterations are needed. The adjustment produces the adjusted value, ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$, and its uncertainty, $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}})$. In addition, it produces several other useful quantities, such as the normalised residuals (6) $\frac{f_{\mathrm{ex}}^{ij}-({\hat{f}}_{\mathrm{ex},\mathrm{sa}}+f_{\mathrm{th},\mathrm{hf}}^{ij})}{u(f_{\mathrm{ex}}^{ij})},$(6) as well as the ${\chi }^2$ value [25]. The associated number of degrees of freedom equals the number of observational equations (M) minus the number of adjusted constants (which equals 1 in this case, so that we have M−1 degrees of freedom).

In Appendix 2, a different approach to incorporating correlation coefficients in the least-squares adjustment is discussed. This method explicitly considers the errors of the theoretical spin coefficients, $\delta E_k^{\rho }$, and includes these as adjustable parameters. This approach is similar to that employed by CODATA, and it produces similar results for the adjusted values of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$ and $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}})$, as well as the ${\chi }^2$ value.

2.4. Cases considered in this work

Here, we will consider two recent cases, namely the terahertz rotational spectroscopy of ${\mathrm{HD}}^{+}$ as reported by Alighanbari et al. [12], and the laser vibrational spectroscopy of ${\mathrm{HD}}^{+}$ by Patra et al. [11]. Alighanbari et al. measured the transition frequencies of M = 6 hyperfine components of the $(v,N)$: $(0,0)\to (0,1)$ transition. We will contract the previously used superscript ij into a single superscript $n=1\cdots 6$ to indicate the six transition frequencies $f_{\mathrm{ex}}^1\cdots f_{\mathrm{ex}}^6$. In order to determine $f_{\mathrm{ex},\mathrm{sa}}^{v^{″}{v}^{\prime }{N}^{″}{N}^{\prime }}=f_{\mathrm{ex},\mathrm{sa}}^{0001}$, Alighanbari et al. combined the six frequencies into a composite frequency, $f_{\mathrm{c}}$, given by $f_{\mathrm{c}}=\sum _n^6{b}_n{f}_{\mathrm{ex}}^n.$With the condition $\sum b_n=1$, this can be written as [12] (7) $f_{\mathrm{c}}=f_{\mathrm{ex},\mathrm{sa}}^{0001}+\sum _n^6{b}_n{f}_{\mathrm{th},\mathrm{hf}}^n.$(7) Owing to the tracelessness of the effective spin Hamiltonian, the terms in the sum in Equation (7(7) $f_{\mathrm{c}}=f_{\mathrm{ex},\mathrm{sa}}^{0001}+\sum _n^6{b}_n{f}_{\mathrm{th},\mathrm{hf}}^n.$(7) ) largely cancel. Moreover, the uncertainty introduced through the sum is minimised by an appropriate choice of the coefficients $b_n$, which are determined numerically using a procedure that is conceptually similar to a least-squares minimisation. From Equation (7(7) $f_{\mathrm{c}}=f_{\mathrm{ex},\mathrm{sa}}^{0001}+\sum _n^6{b}_n{f}_{\mathrm{th},\mathrm{hf}}^n.$(7) ), the value of $f_{\mathrm{ex},\mathrm{sa}}^{0001}$ is subsequently obtained [12]. Alighanbari et al. also argued that from a theoretical viewpoint, strong correlations between some spin coefficients may be suspected. However, they did not quantify the effect of such possible correlations on the found value and uncertainty of $f_{\mathrm{ex},\mathrm{sa}}^{0001}$. Below, we will determine $f_{\mathrm{ex},\mathrm{sa}}^{0001}$ in an alternative way, through a direct least-squares adjustment to $f_{\mathrm{ex}}^1\cdots f_{\mathrm{ex}}^6$, thus enabling a quantitative treatment of the effect of correlations.

The other work considered here, that of Patra et al., involves measurements of M = 2 hyperfine components of the $(v,N)$: $(0,3)\to (9,3)$ transition. Rather than performing an adjustment, they determined $f_{\mathrm{ex},\mathrm{sa}}^{0933}=f_{\mathrm{ex}}^n-f_{\mathrm{th},\mathrm{hf}}^n$ for each value of n = 1, 2, and then computed the weighted mean of the two resulting values of $f_{\mathrm{ex},\mathrm{sa}}^{0933}$. However, contrary to expectations, they found the two values $f_{\mathrm{ex},\mathrm{sa}}^{0933}$ to be in disagreement by about 7 kHz or $4\sigma$. Patra et al. subsequently expanded the theoretical uncertainties $u(E_k^{\rho })$ by a uniform factor of $\eta \approx 2.2$ to reduce the discrepancy to $2\sigma$, which is commonly considered an 'acceptable' deviation. In this way, the $4\sigma$ discrepancy became absorbed in the larger resulting uncertainty of $f_{\mathrm{ex},\mathrm{sa}}^{0933}$. Patra et al. did take into account correlations due to the shared dependence of the theoretical hyperfine shifts $f_{\mathrm{th},\mathrm{hf}}^n$ on the same set of spin coefficients $E_k^{\rho }$, but they did not consider correlations between the spin coefficients themselves [i.e. they effectively considered a scenario ($r_A,r_B,r_C)=(0,0,0)$]. In this work, we will re-analyse the results of Patra et al. using a least-squares adjustment to find $f_{\mathrm{ex},\mathrm{sa}}$, this time also considering the effect of (different types) of correlation.

3. Results

We have analysed the two cases in the following ways. First, we have performed the adjustment for $2^3=8$ possible correlation scenarios $(r_A,r_B,r_C)$, with each of the coefficients $r_A$, $r_B$, $r_C$ set to either 0 or 1. These scenarios include one scenario described by Alighanbari et al., who pointed out that correlation coefficients of type B are expected to be close to 1 for spin coefficients $E_k^{{\rho }^{″}},E_k^{{\rho }^{\prime }}$ [12]. At the same time, we emphasise here that estimated theoretical errors are by definition uncertain, as are their correlation coefficients. For example, past research has shown that estimated errors sometimes do not include all the physics that plays a significant role, including unknown state-dependent error contributions that vary with v, N [24]. Such errors could reduce the type-B correlation identified by Alighanbari et al. To fully appreciate the possible effect of correlations, one should therefore consider other scenarios for which $r_A$, $r_B$, and $r_C$ have different independent values, potentially in a range as wide as $[-1,1]$.

We address this issue in the second step of our analysis with the help of Monte-Carlo simulations, where each of the coefficients $r_A,r_B,r_C$ is selected pseudo-randomly from the interval $[-1,1]$, after which the adjustment is performed. Repeating this 50,000 times provides insight in the distribution of adjusted values ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$ and, through the ${\chi }^2$ value of each adjustment, the relative probability of finding the value ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$ for the corresponding coefficients $r_A,r_B,r_C$, and given the set of observed transition frequencies. This probability is obtained by inserting the found value of ${\chi }^2$ in the ${\chi }^2$ probability density function for $\nu =M-1$ degrees of freedom, i.e. we compute the value of $\frac{1}{\mathrm{\Gamma }(\nu /2)}{\left(\frac{x}{2}\right)}^{\frac{\nu }{2}}\frac{e^{-x/2}}{x},$with Γ the gamma function, by inserting $x={\chi }^2$. The thus found probabilities will later be used as (normalised) weights, w. Note that as an alternative to the Monte-Carlo runs, we could have carried out the same calculations based on a three-dimensional grid of evenly spaced values $r_A,r_B,r_C$. However, such a grid leads to an undesired rasterization of collections of data points in two-dimensional charts, which reduces the resolution of such charts (as compared to collections of data points that are quasi-continuously distributed).

To facilitate the comparison with earlier published results, the theoretical hyperfine structure and corresponding covariances were evaluated using the spin coefficients published by Alighanbari et al. [12] and Patra et al. [11], and not with the updated spin coefficients that were published more recently [21,22].

3.1. Results $(v,N)$: $(0,3)\to (9,3)$ transition

Figure 1 shows the results of the adjustment for eight scenarios involving discrete values (0 or 1) of the correlation coefficients $r_A,r_B,r_C$. For the $(v,N)$: $(0,3)\to (9,3)$ transition [Figure 1(a,b)] a new expansion factor η was determined before each adjustment. The procedure to determine η used here differs from that deployed by Patra et al. In keeping with CODATA practice, we now uniformly expand the entire covariance matrix (expanding not only the theoretical but also the experimental uncertainty) by a factor ${\eta }^2$ such that none of absolute values of the normalised residuals exceeds 2 [8]. Most expansion factors η are between 1.5 and 2.7, depending on the correlation scenario, with one exception of $\eta =5.5$ for the $(0,1,1)$ case, which corresponds to a physically unlikely combination of correlation coefficients.

Figure 1. Adjustment of $f_{\mathrm{ex},\mathrm{sa}}^{0933}$ and $f_{\mathrm{ex},\mathrm{sa}}^{0001}$ for various correlation scenarios [labelled by the tuples $(r_A,r_B,r_C)$]. The corresponding theoretical predictions $f_{\mathrm{th},\mathrm{sa}}$ are those from Refs. [11,12] and are indicated by the vertical dashed lines, while the shaded area represents the corresponding theoretical uncertainty $u(f_{\mathrm{th},\mathrm{sa}})$. For completeness, a more recent value of $f_{\mathrm{th},\mathrm{sa}}$ by Korobov and Karr [7] is shown as a separate data point. (a) Results for the $(v,N)$: $(0,3)\to (9,3)$ transition, considering correlations only for $E_4^{\rho }$ and $E_5^{\rho }$. The value of $f_{\mathrm{ex},\mathrm{sa}}^{0933}$ as determined by Patra et al. [11] is also indicated. (b) Same as (a), but for the $(v,N)$: $(0,0)\to (0,1)$ transition using the data of Alighanbari et al. [12]. (c) Same as (a), but now considering $E_1^{\rho }$ (type-B correlation only), and the pairs $E_4^{\rho },E_5^{\rho }$ and $E_6^{\rho },E_7^{\rho }$. (d) Same as (b), but now including also type-A correlations between $E_6^{\rho },E_7^{\rho }$.

In Figure 1(a), only correlations involving $E_4^{\rho }$ and $E_5^{\rho }$ are considered. The $(0,0,0)$ scenario in Figure 1(a) corresponds to the correlation assumptions made by Patra et al. The shift of this value relative to that of Patra et al. [11] shows the (small) effect of our modified approach to determining the expansion factor η. Clearly, different correlation scenarios $(r_A,r_B,r_C)$ may lead to additional shifts of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}$ of about one standard uncertainty. If we in addition consider correlations between $E_1^{\rho }$, $E_6^{\rho }$ and $E_7^{\rho }$, then the spread in values ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}$ becomes slightly smaller (if we exclude the $(0,1,1)$ point); see Figure 1(c).

The results of the Monte-Carlo run are shown in Figure 2(a). For all adjustments during this run, we necessarily assumed the same expansion factor ($\eta =1.5$); if we were to determine η for each adjustment individually, then the ${\chi }^2$ values of different adjustments would become incomparable. Out of 50,000 realisations, 299 adjustments were discarded because they produced negative ${\chi }^2$ values or complex-valued uncertainties $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933})$. The cause of these problematic cases can be traced back to the covariance matrix (for the correlation scenario in question), for example when correlations coincidentally lead to a much larger variance for one observational equation than for the other. For each of the remaining 49,701 realisations, the resulting value ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}$ is assigned a weight w as explained above. The distribution of weights is shown in Figure 2(a). Here its worth pointing out the difference between the weights (or probabilities) w, which can be regarded as a measure of the compatibility of the measured frequencies with the particular correlation scenario, and the probabilities that are expressed by the bin heights in the three-dimensional histogram, which do not take the compatibility between measurements and correlation scenario into account. Figure 2(a) also shows the distribution of frequency shifts $\mathrm{\Delta }{\hat{f}}^{0933}$, which for each given combination $(r_A,r_B,r_C)$ are computed as follows: $\mathrm{\Delta }{\hat{f}}^{0933}\equiv {\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}(r_A,r_B,r_C)-{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}(0,0,0).$The minimum and maximum values of $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}$ found during the run are $-5.4$ kHz and 2.5 kHz, respectively, which corresponds to $-5.1$ and 2.4 standard uncertainties. However, as will be discussed in Section 4, the probabilities of such strongly deviating values are small, as is their impact on the weighted mean and standard uncertainty derived from the entire set of 49,701 adjustments. The relationship between weight w and the different types of correlation $r_A,r_B,r_C$ is illustrated in Figure 3. It can be seen that the cases with larger weight occur for positive type-A correlation, and negative type-B and type-C correlation. This forms an indication that the combination of the measured transition frequencies and the theoretical model favours such correlation scenarios.

Figure 2. Results of Monte-Carlo runs. (a) $(v,N)$: $(0,3)\to (9,3)$ transition, 49,701 realisations. (b) $(v,N)$: $(0,0)\to (0,1)$ transition, 49,550 realisations. (c) Zoomed-in portion of the histogram shown in (b).

Figure 3. Weight (or probability) w versus r for the Monte-Carlo run shown in Figure 2(a). Although not visible from the graph, the different coloured dots occur in vertically aligned groups of three, each group corresponding to a single adjustment with correlation coefficients $r_A$ (blue), $r_B$ (green), and $r_C$ (red) (colour online only).

3.2. Results $(v,N)$: $(0,0)\to (0,1)$ transition

Effects of correlation are also visible in the adjustments pertaining to the $(v,N)$: $(0,0)\to (0,1)$ transitions; see Figure 1(b,d). If we consider correlation coefficients with discrete values (0 or 1), then, depending on the correlation scenario, we find shifts by more than three standard uncertainties relative to the uncorrelated case $(0,0,0)$. The value and uncertainty of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$ found for the latter scenario corresponds closely to the value of $f_{\mathrm{ex},\mathrm{sa}}^{0001}$ reported by Alighanbari et al. Figure 1(b,d). This is not unexpected because these authors also assumed uncorrelated spin coefficients in their analysis.

The spread of the adjusted frequencies in Figure 1 appears to be considerably larger for the $(v,N)$: $(0,0)\to (0,1)$ transition than for the $(v,N)$: $(0,3)\to (9,3)$ transition. However, closer inspection of the dependence of the shift on the coefficients r reveals that the shift may drop rapidly already for values only slightly below r = 1. This behaviour is illustrated in Figure 4(b) for the scenario $(1,r_B,1)$ and the $(v,N)$: $(0,0)\to (0,1)$ transition. For values $r_B$ approaching 1, the sudden rise in frequency is accompanied by ${\chi }^2$ values that imply a drop in probability by orders of magnitude (not shown). For the $(0,r_B,1)$ scenario concerning the $(v,N)$: $(0,3)\to (9,3)$ transition, the frequency dependence on $r_B$ is shown in Figure 4(a). In this case, the frequency dependence close to $r_B=1$ is considerably less steep than for the $(v,N)$: $(0,0)\to (0,1)$ case.

Figure 4. Dependence of the adjusted frequency on the correlation coefficient $r_B$. (a) $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}$ versus $r_B$ for the scenario $(0,r_B,1)$, and considering correlations between $E_4^{\rho }$ and $E_5^{\rho }$ (upper curve), and between $E_1^{\rho }$, $E_4^{\rho }$, $E_5^{\rho }$, $E_6^{\rho }$, and $E_7^{\rho }$ (lower curve). For the upper and lower curves, the limiting cases $r_B=0$ and $r_B=1$ are also shown in Figure 1(a) and Figure 1(c), respectively. A uniform expansion factor $\eta =1.5$ was applied. (b) $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$ versus $r_B$ for the scenario $(1,r_B,1)$, and considering correlations between $E_4^{\rho }$ and $E_5^{\rho }$ (lower curve), and between $E_4^{\rho }$, $E_5^{\rho }$, $E_6^{\rho }$, and $E_7^{\rho }$ (upper curve). In this case, for the lower and upper curves, the limiting cases $r_B=0$ and $r_B=1$ are also shown in Figure 1(b) and Figure 1(d), respectively.

The results of the Monte-Carlo run are shown in Figure 2(b) and (c). As before, we define $\mathrm{\Delta }{\hat{f}}^{0001}={\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}(r_A,r_B,r_C)-{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}(0,0,0).$Out of 50,000 realisations, 450 adjustments were discarded because they produced negative ${\chi }^2$ values or complex-valued uncertainties $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001})$. The minimum and maximum values of $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$ found during the run are $-63$ Hz and 320 Hz, respectively, corresponding to $-3.6$ and 18 standard uncertainties. However, as will be discussed in Section 4, the probabilities of such strongly deviating values are small, as is their impact on the weighted mean and standard uncertainty derived from the set of 49,550 adjustments.

A salient feature of the histogram in Figure 2(b) is the ${\chi }^2$-distribution-like 'curve' visible in the $(\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001},w)$-plane, along which most results are found. As can be seen in Figure 5(a), the realisations in this part of the histogram correspond to adjustments with values of $r_A$ spread around zero for values $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}\approx 0$, but with $r_A$-values becoming more concentrated and increasing smoothly to asymptotic values of $\pm 1$ if the 'curve' is followed towards ever larger values of $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$. By contrast, no such smooth behaviour is seen in the values of $r_B,r_C$, which appear scattered over the range $[-1,1]$ independent of the value of $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$.

Figure 5. Shift of $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$ versus values of $r_A$ (blue), $r_B$ (green), and $r_C$ (red), obtained from Monte-Carlo runs (colour online only). Although not visible from the graphs, the different coloured dots occur in vertically aligned groups of three, each group corresponding to a single adjustment with correlation coefficients $(r_A,r_B,r_C)$. (a) Results of a ∼50,000-point Monte-Carlo run involving $E_4^{\rho }$, $E_5^{\rho }$, $E_6^{\rho }$, and $E_7^{\rho }$. (b) Same as in (a), but now involving only $E_4^{\rho }$ and $E_5^{\rho }$. (c) Same as (a), but now involving only $E_6^{\rho }$ and $E_7^{\rho }$. (d) Weight (or probability) w versus r for the run involving $E_4^{\rho }$, $E_5^{\rho }$, $E_6^{\rho }$, and $E_7^{\rho }$. This run was also used to produce Figure 2(b).

To investigate which spin coefficients are responsible for the observed relationship between ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$ and $r_A$, the Monte-Carlo run was repeated once considering only correlations involving $E_4^{\rho }$ and $E_5^{\rho }$, and once involving only $E_6^{\rho }$ and $E_7^{\rho }$. The results in the $(\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001},w)$-plane are presented in Figure 5(b,c), respectively As is visible in Figure 5(c), the latter pair of coefficients is responsible for the behaviour seen in Figure 2(b), with their type-A correlation coefficient displaying a nearly functional dependence on ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$. Plotting the outcomes of the full Monte-Carlo run (i.e. involving $E_4^{\rho }$, $E_5^{\rho }$, $E_6^{\rho }$, and $E_7^{\rho }$) in the $(w,r)$-plane reveals that the weights peak around $r_A=0.76$ and $r_A=-0.87$ [Figure 5(d)].

4. Discussion

The results in Section 3 show that if one uses a theoretical hyperfine structure model to derive the spin-averaged transition frequency from a set of measured hyperfine components of a rovibrational transition, assumptions regarding correlated errors in the theoretical spin coefficients influence the result (in our case the adjusted value of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$). As correlation coefficients between various theory errors cannot be exactly known, some liberty (or rather uncertainty) exists in the choice of correlation coefficients. The consequence of this liberty is twofold: the value of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$ could be significantly shifted away from the value that is obtained when assuming zero correlation (as was done in Refs. [11,12]), and the concomitant spread of values ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$ introduces additional uncertainty to its final value. For example, computing the mean and standard deviation of all values ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}$ in the 49,701-point Monte-Carlo run yields a shift of $\mathrm{\Delta }{\hat{f}}^{0933}=-0.07$ kHz and a standard deviation of 0.58 kHz. While both the shift and the standard deviation are still smaller than $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933})(0,0,0)=1.0$ kHz, the effect cannot be ignored, in particular in view of some of the very large observed shifts $\mathrm{\Delta }{\hat{f}}^{0933}$. Similar considerations apply to the $(v,N)$: $(0,0)\to (0,1)$ transition also studied here: taking the mean and standard deviation of all values ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$ in the 49,550-point Monte-Carlo run yields a shift of $\mathrm{\Delta }{\hat{f}}^{0001}=1.5$ Hz and a standard deviation of 3.6 Hz. Again, both the shift and the standard deviation are smaller than $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001})(0,0,0)=17$ Hz, but also for this case the effect cannot be ignored in view of some of the very large shifts $\mathrm{\Delta }{\hat{f}}^{0001}$. Our analysis furthermore shows that one should be particularly careful with scenarios where correlation coefficients r = 1 are assumed, as evidenced by the very steep frequency shifts shown in Figure 4(b).

4.1. Interpretation and treatment of correlation-induced shifts

We may interpret this effect in the context of observational equations, which generally contain a measured quantity in the left-hand side of the equation, and a corresponding theoretical expression in the right-hand side of the equation. Usually, all parameters of the theoretical model are included in the theoretical expression, and their values are either adjusted, or estimated using external means. In the latter case, the uncertainties of the estimated parameters (and their covariances with other parameters) are included in the covariance matrix of the input data. However, in cases where correlation coefficients are possibly nonzero but a priori unknown, the correlation coefficients themselves become parameters of the theoretical model too.

This line of reasoning implies that correlation coefficients, like all other parameters, should meet two requirements: first, they need to be assigned a value (either by some external means, or by making them adjustable parameters); second, they need be assigned an uncertainty. However, assigning a value to correlation coefficients between estimated theory errors that are intrinsically uncertain (and that may not include all relevant physics) can be considered a 'best guess' at best. The second requirement, assigning an uncertainty range to r values, might be achievable, but would most likely necessitate subjective measures such as Bayesian 'degree of belief' likelihood functions. The other option, adjusting the values of correlation coefficients themselves might be challenging too, as it may not always be possible to find a single-valued functional correspondence between r and the adjusted parameter, and because the dependence on r may exhibit strong nonlinear behaviour. This may preclude linearisation of the theoretical expression and – ultimately – the (linear) least-squares adjustment itself. On the other hand, if correlation coefficients can be adjusted, the adjustment would automatically produce their uncertainties as well.

Here we propose a pragmatic alternative approach to deal with the issues outlined above. The Monte-Carlo runs developed here provide a means to determine the spread in a way that makes no assumptions regarding the values of $r_A,r_B,r_C$ whatsoever. As shown above, we were able to quantify the effect of unknown correlation coefficients by computing the mean and standard deviation of all ∼50,000 values of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$. In addition, we can take advantage of the capability of any least-squares adjustment to provide a (relative) measure – the ${\chi }^2$ value – of the compatibility between the measurements and the theoretical model (which now also encompasses the correlation coefficients), in order to get an improved estimate of ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}$ and its uncertainty. To this end, we use the weights w to calculate means and uncertainties of the weighted results of the Monte-Carlo runs. For example, for the $(v,N)$: $(0,3)\to (9,3)$ transition, we earlier found an unweighted mean (standard deviations within parentheses) of $-0.07(0.58)$ kHz relative to ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}(0,0,0)$, whereas consideration of the weights leads to a compatible but more precise value of $0.03(0.27)$ kHz. This should be compared with the uncertainty $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933})=1.0$ kHz for the uncorrelated $(0,0,0)$ scenario (with expansion factor $\eta =1.5$). Also for the $(v,N)$: $(0,0)\to (0,1)$ transition, the weighting reduces the standard deviation significantly, improving the unweighted mean and uncertainty of all values $\mathrm{\Delta }{\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}$, $1.5(3.6)$ Hz, to the weighted result $1.1(1.5)$ Hz. The standard deviation of the weighted result is negligibly small compared to $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001})(0,0,0)$, which equals 17 Hz. The results are summarised in Table 1 for both the $(v,N)$: $(0,3)\to (9,3)$ and $(v,N)$: $(0,0)\to (0,1)$ transitions.

Table 1. Correlation-induced shift and standard uncertainty, as obtained from unweighted and weighted results of the Monte-Carlo runs, to the $(v,N)$: $(0,3)\to (9,3)$ and $(v,N)$: $(0,0)\to (0,1)$ spin-averaged transition frequencies.

	$(v,N)$: $(0,3)\to (9,3)$		$(v,N)$: $(0,0)\to (0,1)$
Method	Shift${}^{\mathrm{a}}$ (kHz)	Uncertainty (kHz)	Shift${}^{\mathrm{b}}$ (Hz)	Uncertainty (Hz)
Monte-Carlo unweighted	−0.07	0.58	1.5	3.6
Monte-Carlo weighted	0.03	0.27	1.1	1.5

${}^{\mathrm{a}}$Shift with respect to ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0933}(0,0,0)=415,264,925,500.7(1.0)$ kHz.

${}^{\mathrm{b}}$Shift with respect to ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{0001}(0,0,0)=1,314,925,752,911(17)$ Hz.

Furthermore, the three-dimensional histogram in Figure 2(b) indicates that most values $\mathrm{\Delta }{\hat{f}}^{0001}$ are located along the ${\chi }^2$-like curve with frequency offsets between 0 Hz and 2.5 Hz. Two peaks are visible: a primary one around $(\mathrm{\Delta }{\hat{f}}^{0001},w\times {10}^3)=(0\mathrm{Hz},0.016)$, and a secondary one around $(2\mathrm{Hz},0.028)$. From Figure 5(a) it can be seen that for the secondary peak (i.e. for $\mathrm{\Delta }{\hat{f}}^{0001}\approx 2$ Hz), $r_A$ attains values close to 0.76 or $-0.87$. These values are the ones with the largest weight w in Figure 5(d), a behaviour that is due to the coefficients $E_6^{\rho }$ and $E_7^{\rho }$. We therefore conclude that the transition frequencies measured by Alighanbari et al. are compatible with a scenario where the coefficients $E_6^{{\rho }^{\prime }}$ and $E_7^{{\rho }^{\prime }}$ are rather strongly correlated. Given the very similar theoretical models that are used to compute these physically related coefficients, one would indeed expect a strong correlation between their errors.

Yet, Figure 2(b) and Figure 5(a,d) show that the measured transition frequencies are also compatible with $\mathrm{\Delta }{\hat{f}}^{0001}\approx 0$ Hz for relatively small values of $r_A$, and for a wide range of values $r_B,r_C$. For these scenarios, which constitute the primary peak in Figure 2(b), the physical interpretation is less obvious. The occurrence of this peak is made possible by the fact that the Monte-Carlo run does not make a priori assumptions regarding the values of $r_A,r_B,r_C$, even if many combinations $(r_A,r_B,r_C)$ might be less likely from a physical point of view. The $\sim 2$ Hz distance between the primary and the secondary peak is therefore indicative of the uncertainty introduced through the unknown correlation coefficients.

It may furthermore be possible to develop a refined approach which includes theoretical input on the possible values of correlation coefficients. For example, anticorrelations appear very unlikely, as it has never been observed that a given NRQED correction term changes sign between two different rovibrational states (at least within the range of investigated rovibrational states studied so far). Such auxiliary information may in the future be used to truncate the range of possible values of the correlation coefficients.

4.2. Adjustment of fundamental physical constants

The correlation-induced shift and uncertainty have an effect on the value of $m_p/m_e$ that can be derived from the $(v,N)$: $(0,3)\to (9,3)$ and $(v,N)$: $(0,0)\to (0,1)$ transition frequencies. We calculate their impact using sensitivity coefficients computed by Karr et al. [10]. To find the fractional effect (in ppt) on the value of $m_p/m_e$, we multiply the entries for the $(v,N)$: $(0,3)\to (9,3)$ transition in Table 1 with $-6.8$ kHz${}^{-1}$, and those for the $(v,N)$: $(0,0)\to (0,1)$ with $-0.77$ Hz${}^{-1}$. The results are shown in Table 2, and we find that the correlation-induced shift and uncertainty affect the value of $m_p/m_e$ at the few-ppt level. This is still one order of magnitude smaller than the $\sim 20$-ppt uncertainties that were reported by Patra et al. and Alighanbari et al. [11,12], which are dominated by the uncertainty of the theoretical values of the spin-averaged transition frequencies, $f_{\mathrm{th},\mathrm{sa}}^{0933}$ and $f_{\mathrm{th},\mathrm{sa}}^{0001}$, respectively.

Table 2. Correlation-induced shift and standard uncertainty (in ppt), to the value of $m_p/m_e$, as obtained from unweighted and weighted results of the Monte-Carlo runs.

	$(v,N)$: $(0,3)\to (9,3)$		$(v,N)$: $(0,0)\to (0,1)$
Method	Relative shift	Relative uncertainty	Relative shift	Relative uncertainty
Monte-Carlo unweighted	0.5	3.9	$-1.2$	2.8
Monte-Carlo weighted	$-0.2$	1.9	$-0.8$	1.1

We intentionally do not provide absolute values of $m_p/m_e$ here: since the publication of the three recent measurements in ${\mathrm{HD}}^{+}$ [11,12,14], no less than nine different values of $m_p/m_e$ derived from these measurements have been reported. These were computed using various methods and using different choices of fundamental constants used in the ${\mathrm{HD}}^{+}$ theory, and the values also depend somewhat on the state of the art of the theory at the time of publication [7,11,12,14]. Instead of extending the list of $m_p/m_e$ determinations from ${\mathrm{HD}}^{+}$ here by four new values, each with its own specific interpretation, we refer the reader to the overview of ${\mathrm{HD}}^{+}$-derived values of $m_p/m_e$ by Korobov and Karr, who based their analysis on the most recent theory [7]. It is conceivable that these values will eventually be superseded by the anticipated CODATA-2022 adjustment, which would produce a single and thoroughly justified value of $m_p/m_e$, which would be based on all relevant measurements and theory (including not only ${\mathrm{HD}}^{+}$ and Penning-trap measurements, but also measurements and theory of atomic hydrogen and hydrogen-like systems for the Rydberg constant and nuclear charge radii, and measurements relevant to the fine-structure constant).

Inclusion of ${\mathrm{HD}}^{+}$ data in a global adjustment of fundamental physical constants such as the CODATA adjustment could proceed in the following way, by augmenting the global set of observational equations used for the adjustment with the following set of equations: (8) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{th},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}(m_p/m_e,m_d/m_p,R_{\mathrm{\infty }},R_p,R_d,\alpha )\\ & +W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}})\end{array}$(8) To arrive at Equation (8(8) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{th},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}(m_p/m_e,m_d/m_p,R_{\mathrm{\infty }},R_p,R_d,\alpha )\\ & +W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}})\end{array}$(8) ), we replaced $f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}$ in Equation (4(4) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}}),\end{array}$(4) ) with the known function $f_{\mathrm{th},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}(m_p/m_e,m_d/m_p,R_{\mathrm{\infty }},R_p,R_d,\alpha )$ [10]. Here, $R_{\mathrm{\infty }}$ stands for the Rydberg constant, α for the fine-structure constant, and $R_p$ and $R_d$ for the proton and deuteron electric charge radii, respectively. Equation (8(8) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{th},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}(m_p/m_e,m_d/m_p,R_{\mathrm{\infty }},R_p,R_d,\alpha )\\ & +W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}})\end{array}$(8) ) would include the two $(v,N)$: $(0,3)\to (9,3)$ hyperfine components as well as the six $(v,N)$: $(0,0)\to (0,1)$ hyperfine components, and may furthermore include the two hyperfine components of the $(v,N)$: $(0,0)\to (1,1)$ transition measured by Kortunov et al. [14]. The arguments of the functions $f_{\mathrm{th},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}$ are parameters that are globally adjusted (i.e. they are adjusted to many other observational equations besides those from ${\mathrm{HD}}^{+}$, derived from measurements of other physical systems). Uncertainties and correlations can be taken into account through a covariance matrix as before, and the impact of the various correlation coefficients could be quantitatively assessed using the approach presented in this work. Alternatively, theory errors and their mutual correlations can be included by introducing additional adjustable parameters that explicitly represent the error terms [8,25]; see also Appendix 2.

5. Conclusion

In conclusion, we have shown that if hyperfine shifts are to be removed from measured rovibrational transtion frequencies by using theoretically calculated hyperfine structure, one needs to consider the effect of (potentially) correlated parameters in the theoretical model. To our knowledge, this constitutes a hitherto unrecognised systematic effect in molecular spectroscopy, which we here analysed using a least-squares approach that can readily take such correlations into account. We analysed the impact of this systematic effect on two recently published results obtained from high-precision spectroscopy of HD${}^{+}$ and calculated hyperfine shifts derived from an effective spin Hamiltonian. We showed that certain correlation scenarios may lead to strongly deviating values of spin-averaged rovibrational transition frequencies. To quantify the effect in terms of the corresponding frequency shift and uncertainty, we presented a Monte-Carlo approach that enables statistical averaging over many possible correlation scenarios. With the help of weights, derived from ${\chi }^2$ values from the least-squared adjustment, a robust estimate of the size and uncertainty of this new systematic shift was obtained, resulting in values much smaller than the uncertainty reported earlier for the spin-averaged transition frequencies [11,12]. This work may be of interest to future spectroscopic studies which make use of theoretical calculations to remove hyperfine shifts, as well as potential future applications of the spectroscopic data analysed here, such as the adjustment of fundamental physical constants.

Acknowledgments

The author wishes to thank Jean-Philippe Karr for fruitful discussions and comments on the manuscript.

Disclosure statement

The author declares no conflicts of interest.

Appendices

Appendix 1: Overview of nonzero correlation coefficients

Figure A1 provides a graphical overview of the correlation coefficients that are assigned nonzero values, and how they relate to the different types of correlation identified in this work.

Figure A1. Graphical overview of correlation coefficients. Empty cells correspond to correlation coefficients which are always set to zero, either because the corresponding covariances are negligibly small on an absolute frequency scale, or because the spin coefficients involved are unlikely to be uncorrelated as they represent physically distinct interactions. Nonzero correlations of type A, B and C are indicated. Cells with a shaded background correspond to correlations that are nonzero only in the $(v,N)$: $(0,3)\to (9,3)$ analysis.

Appendix 2: Including theory errors as adjustable parameters

It is possible to formulate the least-squares problem in a slightly different form, which is due to CODATA. In this case, we define 18 additional adjustable parameters, $\delta E_k^{\rho }$, which represent the errors of the 18 spin coefficients $E_k^{\rho }$ for the lower and upper states. We propagate and add these errors to the observational equations [Equation (4(4) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}}),\end{array}$(4) )] as follows: (A1) $\begin{array}{rl}{f}_{\mathrm{ex}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}& \doteq f_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }}+W_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}(E_1^{{\rho }^{\prime }},\dots ,E_9^{{\rho }^{\prime }})\\ & -W_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}(E_1^{{\rho }^{″}},\dots ,E_9^{{\rho }^{″}})\\ & +\sum _{k=1}^9(\frac{\mathrm{\partial }{W}_{\mathrm{hf}}^{{\rho }^{\prime }{\beta }^{\prime }}}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\delta E_k^{{\rho }^{\prime }}-\frac{\mathrm{\partial }{W}_{\mathrm{hf}}^{{\rho }^{″}{\beta }^{″}}}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}\delta E_k^{{\rho }^{″}}).\end{array}$(A1) Furthermore, we add 18 observational equations of the form $0\doteq \delta E_k^{\rho }$. Note that for the $(v,N)$: $(0,0)\to (0,1)$ transition, the number of spin coefficients is only 11, and the number of parameters $\delta E_k^{\rho }$ (as well as the number of corresponding observational equations) is reduced accordingly. The covariance matrix of input data is now composed of the individual terms from the summations in Equation (5(5) $\begin{array}{rl}& \mathrm{cov}(f_1,f_2)\\ & =\sum _{k,l=1}^9[\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}r(E_k^{{\rho }^{″}},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{″}})u(E_l^{{\rho }^{″}})\\ & +\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{\prime }}}r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{\prime }})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{\prime }})\\ & +(\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}+\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_l^{{\rho }^{″}}}\frac{\mathrm{\partial }{f}_2}{\mathrm{\partial }{E}_k^{{\rho }^{\prime }}})\\ & \times r(E_k^{{\rho }^{\prime }},E_l^{{\rho }^{″}})u(E_k^{{\rho }^{\prime }})u(E_l^{{\rho }^{″}})\phantom{\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }{E}_k^{{\rho }^{″}}}}].\end{array}$(5) ). Performing this alternative least-squares adjustment results into values ${\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }}$, $u({\hat{f}}_{\mathrm{ex},\mathrm{sa}}^{{\rho }^{″}{\rho }^{\prime }{\beta }^{″}{\beta }^{\prime }})$, and ${\chi }^2$ that are identical to the results obtained using the approach of Section 2. In addition, values of the errors $\delta E_k^{\rho }$ are produced. This alternative method, however, produces singular covariance matrices V for some scenarios where $r_A$, $r_B$, and/or $r_C$ equal 1. These matrices cannot be inverted, and therefore in such cases the adjustment cannot be carried out. For this reason, the least-squares approach of Section 2 was preferred in this work.