# Higher-order corrections to the spin-orbit and spin-spin tensor interactions
in HD+

Mohammad Haidar1111Present affiliations: 1. Sorbonne Université, Laboratoire
de Chimie Théorique (LCT), 4 place Jussieu, F-75005 Paris, France 2. Sorbonne
Université, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions
(LJLL), 4 place Jussieu, F-75005 Paris, France 3. TotalEnergies, Tour Coupole,
2 Pl. Jean Millier, F-92078 Paris la Défense, France, Vladimir I. Korobov2,
Laurent Hilico1,3, and Jean-Philippe Karr1,3 1Laboratoire Kastler Brossel,
Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, 4 place
Jussieu, F-75005 Paris, France  2Bogoliubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Russia  3Université
d’Evry-Val d’Essonne, Université Paris-Saclay, Boulevard François Mitterrand,
F-91000 Evry, France

###### Abstract

Improved values of hyperfine coefficients related to the electronic spin-orbit
and electron-nucleus spin-spin tensor interactions in the HD+ molecular ion
are obtained through numerical calculation of relativistic corrections at the
$m\alpha^{6}$ order and radiative corrections at the $m\alpha^{7}\ln(\alpha)$
order. The theoretical accuracy is improved by more than one order of
magnitude. Some deviations with recent high-precision ro-vibrational
spectroscopy experiments are observed, in contrast with the good agreement
obtained in H${}_{2}^{+}$.

##  I Introduction

Laser spectroscopy of HD+ has recently made a significant leap in precision,
as a rotational transition Alighanbari20  and two vibrational transitions
Patra20 ; Kortunov21  have been probed in the Lamb-Dicke regime on ensembles
of trapped and sympathetically cooled molecular ions. Transition frequencies
were measured with relative uncertainties in the 10-11-10-12 range, and
comparison with theoretical predictions Korobov17 ; Korobov21  was used to get
improved determinations of the proton-electron mass ratio and constraints on a
“fifth force” between the nuclei Germann21 .

The hyperfine splitting of HD+ rovibrational levels gives rise to 4, 10 or 12
sublevels for $L=0$, $L=1$ and $L\geq 2$ states respectively, where $L$ is the
rotational quantum number Bakalov06 . In the above-mentioned experiments, a
few hyperfine components of the rovibrational transitions (between 2 and 6)
were measured with sub-kHz uncertainties, from which a “spin-averaged”
transition frequency was extracted using theoretical predictions of the
hyperfine structure. A discrepancy between theory and experiment was evidenced
in the hyperfine slitting of the $(v=0,L=3)\to(v^{\prime}=9,L^{\prime}=3)$
transition Karr20 , which affects the uncertainty of the extracted spin-
averaged transition frequency Patra20 ; Koelemeij22 . Beyond that, constant
progress in experimental accuracy opens the way to highly precise tests of the
hyperfine structure theory. It is also worth noting that the HD+ hyperfine
structure is sensitive to the deuteron’s electric quadrupole moment $Q_{d}$. A
value of $Q_{d}$ was extracted from the rotational transition measurement
Alighanbari20 , although it is significantly less precise than that obtained
from magnetic resonance measurements in $D_{2}$ Code71 ; Puchalski20 . For all
these reasons, it is highly desirable to improve further the hyperfine
structure theory in HD+.

On the theoretical side, all the coefficients of the effective spin
Hamiltonian introduced in Bakalov06  have been calculated in the framework of
the Breit-Pauli Hamiltonian with account of the electron’s anomalous magnetic
moment, which yielded a relative accuracy on the order of $\alpha^{2}$. We
recently improved the precision of the largest hyperfine coefficients, i.e.
the electron-proton and electron-deuteron spin-spin scalar interactions, to
slightly below 1 ppm Karr20 . The next largest coefficients are those related
to the electron spin-orbit (denoted by $E_{1}$ in Bakalov06 ) and the
electron-proton ($E_{6}$) and electron-deuteron ($E_{7}$) spin-spin tensor
interactions. The effective Hamiltonian describing $m\alpha^{6}$-order
relativistic corrections to these interactions was derived in Korobov20  in
the framework of nonrelativistic quantum electrodynamics (NRQED), and
preliminary numerical results were given. In Haidar22 , the theory was further
refined by including radiative corrections at the following order
($m\alpha^{7}\ln(\alpha)$). The spin-orbit and spin-spin interaction
coefficients were calculated numerically in the H${}_{2}^{+}$ ion for a range
of rovibrational states, and good agreement with available rf spectroscopy
data was observed. In the present work, we report on the numerical calculation
of these corrections in HD+, which improves the theoretical accuracy of the
$E_{1}$, $E_{6}$ and $E_{7}$ coefficients by almost one order of magnitude
with respect to Ref. Bakalov06 . We then compare our theoretical predictions
of the hyperfine splitting with available experimental data Alighanbari20 ;
Patra20 ; Kortunov21 . Significant deviations are evidenced, possible causes
of which are discussed in the last section.

##  II Numerical results

Higher-order corrections to the spin-orbit and spin-spin tensor interaction
coefficients in the hydrogen molecular ions were studied in our previous works
Korobov20 ; Haidar22 . The expressions of relativistic corrections at orders
$m\alpha^{6}$ can be found in Haidar22 , Eqs. (17)-(24), and radiative
corrections at the $m\alpha^{7}\ln(\alpha)$ order are given in Eqs. (27) and
(29-35). In Tables 1 , 2, and 3, we report the results of numerical
calculations of these correction terms for a few rovibrational states of HD+
of direct interest for experiments Alighanbari20 ; Patra20 ; Kortunov21 ;
Zhong15 . Complete results for the rovibrational states $(1\leq L\leq 4,0\leq
v\leq 9)$ are given in the Appendix F. These numerical calculations are
performed following the methods presented in Korobov20 ; Haidar22  relying on
a variational expansion of the three-body wavefunction in terms of
exponentials of interparticle distances, with pseudo-randomly chosen exponents
Thakkar77 ; Korobov00 .

The numerical and theoretical uncertainties were discussed in Haidar22  in the
context of the H${}_{2}^{+}$ ion, and the main conclusions remain valid for
the results presented here. Briefly, the numerical uncertainty of $E_{1}$
(resp. $E_{6},E_{7}$) is estimated to be smaller than 10 Hz (resp. 1 Hz), see
Haidar22  from a study of convergence of a few terms. The overall uncertainty
is dominated by yet unevaluated nonlogarithmic $m\alpha^{7}$-order
contributions, and is estimated to about one third of the total contribution
of order $m\alpha^{7}\ln(\alpha)$. The resulting uncertainties amount to 3-4
ppm for the electronic spin-orbit interaction and 2 ppm for tensor
interactions respectively. This represents an improvement by factors of about
15 for $E_{1}$, and 25 for $E_{6}$ and $E_{7}$, over the results of Ref.
Bakalov06 .

$(L,v)$ | $\displaystyle E_{1}^{(BP)}$ | $\displaystyle\mathcal{U}_{Y_{1}}$ | $\displaystyle\mathcal{U}_{W}$ | $\displaystyle\mathcal{U}_{CM}$ | $\displaystyle\Delta E_{so\hbox{-}H_{B}}$ | $\displaystyle\Delta E_{so\hbox{-}so}^{(1)}$ | $\displaystyle\Delta E_{so\hbox{-}ret}$ | $\displaystyle\Delta E_{1}^{(6)}$ | $\displaystyle\mathcal{U}_{Y_{1}}$ | $\displaystyle\mathcal{U}_{q^{2}}$ | $\displaystyle\Delta E_{so\hbox{-}H_{(5\ln)}}$ | $\displaystyle\Delta E_{1}^{(7\ln)}$ |  $\displaystyle E_{1}$ (this work)
---|---|---|---|---|---|---|---|---|---|---|---|---|---
(1,0) | 31984.645 | 1.170 | -2.736 | 0.021 | 2.088 | 0.312 | 0.257 | 1.112 | -0.026 | 0.045 | -0.367 | -0.348 | 31985.41(12)
(1,1) | 30280.027 | 1.110 | -2.631 | 0.027 | 1.986 | 0.299 | 0.243 | 1.035 | -0.025 | 0.044 | -0.345 | -0.326 | 30280.74(11)
(1,5) | 24076.867 | 0.893 | -2.205 | 0.042 | 1.592 | 0.233 | 0.193 | 0.747 | -0.020 | 0.036 | -0.268 | -0.252 | 24077.36(8)
(1,6) | 22643.449 | 0.834 | -2.097 | 0.044 | 1.498 | 0.219 | 0.181 | 0.679 | -0.019 | 0.035 | -0.251 | -0.235 | 22643.89(8)
(3,0) | 31627.353 | 1.156 | -2.694 | 0.019 | 2.042 | 0.308 | 0.255 | 1.086 | -0.026 | 0.045 | -0.360 | -0.341 | 31628.10(11)
(3,9) | 18270.560 | 0.680 | -1.732 | 0.043 | 1.158 | 0.182 | 0.146 | 0.478 | -0.015 | 0.028 | -0.198 | -0.185 | 18270.85(6)
Table 1: Corrections to the spin-orbit interaction coefficient $E_{1}$ for a few rovibrational states of HD+ (in kHz). The leading-order (Breit-Pauli) value $E_{1}^{(BP)}$ (Ref. Bakalov06 ) is given in column 2. Columns 3-5 and 6-8 are respectively the first-order and second-order contributions [Eqs. (18) and (19) of Haidar22 ] at the m$\alpha^{6}$ order, and the total correction at this order, $\Delta E_{1}^{(6)}$, is given in column 9. Columns 10-11 are the first-order contributions at the $m\alpha^{7}\ln(\alpha)$ order [Eqs. (27) and (29) of Haidar22 ] whereas column 12 is the second-order term [Eq. (32)] of Haidar22 ]. The total correction at this order, $\Delta E_{1}^{(7\ln)}$, is given in column 13. The last column is our value of $E_{1}$. Its estimated uncertainty (equal to one third of $\Delta E_{1}^{(7\ln)}$) is indicated between parentheses. $(L,v)$ | $\displaystyle E_{6}^{(BP)}$ | $\displaystyle\mathcal{U}_{W}^{(2)}$ | $\displaystyle\mathcal{U}_{CM}^{(2)}$ | $\displaystyle\Delta E_{ss\hbox{-}H_{B}}^{(2)}$ | $\displaystyle\Delta E_{6}^{(6)}$ | $\displaystyle\Delta E_{6}^{(7\ln)}$ |  $\displaystyle E_{6}$ (this work)
---|---|---|---|---|---|---|---
(1,0) | 8611.112 | -0.806 | 0.093 | 0.955 | 0.242 | -0.055 | 8611.299(18)
(1,1) | 8136.686 | -0.770 | 0.090 | 0.905 | 0.225 | -0.052 | 8136.859(17)
(1,5) | 6421.232 | -0.632 | 0.079 | 0.726 | 0.172 | -0.041 | 6421.364(14)
(1,6) | 6027.810 | -0.599 | 0.076 | 0.678 | 0.154 | -0.039 | 6027.925(13)
(3,0) | 948.5222 | -0.0883 | 0.0101 | 0.1042 | 0.0260 | -0.0060 | 948.5421(20)
(3,9) | 538.9906 | -0.0544 | 0.0070 | 0.0593 | 0.0119 | -0.0035 | 539.9991(12)
Table 2: Corrections to the electron-proton spin-spin tensor interaction coefficient $E_{6}$ for a few rovibrational states of HD+ (in kHz). The leading-order (Breit-Pauli) value $E_{6}^{(BP)}$ (Ref. Bakalov06 ) is given in column 2. Columns 3-4 and 5 are respectively the first-order and second-order contributions [Eqs. (22) and (23) of Haidar22 ] at the $m\alpha^{6}$ order. The total correction at this order, $\Delta E_{6}^{(6)}$, is given in column 6. Column 7 is the second-order contribution at the m$\alpha^{7}\ln(\alpha)$ order [Eq. (35) of Haidar22 ]. The last column is our value for $E_{6}$. Its estimated uncertainty (equal to one third of $\Delta E_{6}^{(7\ln)}$) is indicated between parentheses.  $(L,v)$ | $\displaystyle E_{7}^{(BP)}$ | $\displaystyle\mathcal{U}_{W}^{(2)}$ | $\displaystyle\mathcal{U}_{CM}^{(2)}$ | $\displaystyle\Delta E_{ss\hbox{-}H_{B}}^{(2)}$ | $\displaystyle\Delta E_{7}^{(6)}$ | $\displaystyle\Delta E_{7}^{(7\ln)}$ |  $\displaystyle E_{7}$ (this work)
---|---|---|---|---|---|---|---
(1,0) | 1321.7673 | -0.1237 | 0.0142 | 0.1467 | 0.0371 | -0.0085 | 1321.7960(28)
(1,1) | 1248.9359 | -0.1182 | 0.0138 | 0.1390 | 0.0346 | -0.0080 | 1248.9624(27)
(1,5) | 985.5868 | -0.0971 | 0.0121 | 0.1115 | 0.0264 | -0.0063 | 985.6069(21)
(1,6) | 925.1895 | -0.0919 | 0.0115 | 0.1040 | 0.0237 | -0.0060 | 925.2072(20)
(3,0) | 145.5938 | -0.0136 | 0.0016 | 0.0160 | 0.0040 | -0.0009 | 145.5969(3)
(3,9) | 82.7237 | -0.0083 | 0.0010 | 0.0091 | 0.0018 | -0.0005 | 82.7250(2)
Table 3: Same as Table 2, for the electron-deuteron spin-spin tensor
coefficient $E_{7}$.

##  III Comparison with experiments and discussion

###  III.1 Theoretical hyperfine splitting of rovibrational transitions

These results allow us to improve theoretical predictions of the hyperfine
splitting of rovibrational transitions that were recently measured with high
accuracy. These are obtained by diagonalizing the effective spin Hamiltonian
of Ref. Bakalov06 , which comprises nine coefficients $E_{1}\ldots E_{9}$. The
values of $E_{1}$, $E_{6}$ and $E_{7}$ are taken from the present work. For
the largest coefficients, i.e the electron-proton and electron-deuteron spin-
spin scalar interaction coefficients (“Fermi” contact interaction), $E_{4}$
and $E_{5}$, we use the very precise values computed in Karr20 . Finally, the
smaller coefficients $E_{2}$, $E_{3}$, $E_{8}$ and $E_{9}$ are calculated in
the framework of the Breit-Pauli Hamiltonian Bakalov06 . We have updated the
$E_{9}$ coefficient, which corresponds to the effect of the deuteron’s
quadrupole moment $Q_{d}$, with the latest and most precise determination of
$Q_{d}$ Puchalski20 . The values of all the coefficients used here can be
found in the Appendix A. The fractional uncertainties of $E_{4}$ and $E_{5}$
are estimated to 0.93 and 0.57 ppm, respectively Karr20 , and those of
$E_{2}$, $E_{3}$ and $E_{8}$ are taken as equal to $\alpha^{2}\simeq 53$ ppm.
Finally, for $E_{9}$ we add quadratically the $\alpha^{2}$ theoretical
uncertainty and the 82 ppm uncertainty of $Q_{d}$ Puchalski20 , yielding a
total fractional uncertainty of 98 ppm.

The theoretical uncertainties of hyperfine intervals $f_{ij}=f_{j}-f_{i}$,
where $i$ ($j$) denotes the lower (upper) hyperfine states, and $f_{i}$,
$f_{j}$ their respective hyperfine shifts, are then calculated by propagating
the uncertainties of the hyperfine coefficients, using the derivatives

$\gamma_{i,k}=\frac{\partial f_{i}}{\partial E_{k}}\,,1\leq k\leq 9.$ (1)

The values of these derivatives for the hyperfine states involved in the
$v=0\to 9$ transition Patra20  are given in the Appendix B. For the states
relevant to the rotational transition and to the $v=0\to 1$ transition, they
can be found in Alighanbari20  and the Supplementary information of Kortunov21
, respectively.

It is important to take correlations into account in order to get reliable
uncertainty estimates. There are strong correlations between theoretical
errors of a given coefficient $E_{k}$ in different rovibrational states, which
can be understood as follows. First, these errors are dominated by yet
uncalculated QED contributions (numerical uncertainties being negligibly
small), which are given by expectation values of the same effective operators,
evaluated with different wavefunctions. Second, the uncalculated terms are
corrections to the bound electron and essentially depend on the electronic
part of the wavefunction, which varies only slowly with the rovibrational
state. Correlations may also exist between errors of different coefficients
$E_{k}$, $E_{l}$ for the same rovibrational state if they originate from
similar QED terms. The hypotheses we have adopted for our uncertainty
estimates are presented in the Appendix C.

(a) $(v=0,L=0)\to(v^{\prime}=1,L^{\prime}=1)$
---
$i$ | $FSJ\rightarrow F^{\prime}S^{\prime}J^{\prime}$ | $j$ | $FSJ\rightarrow F^{\prime}S^{\prime}J^{\prime}$ |  $f_{ij}^{exp}$ Kortunov21  |  $f_{ij}^{theor}$ Kortunov21  |  $f_{ij}^{theor}$ (this work) | $\Delta_{ij}$ | $\Delta_{ij}/\sigma_{c}$
12 | $122\rightarrow 121$ | 16 | $122\rightarrow 123$ | 41 294.06(32) | 41 293.81(44) | 41 293.66(12) | 0.40 | 1.2
(b) $(v=0,L=3)\to(v^{\prime}=9,L^{\prime}=3)$
$i$ | $FSJ\rightarrow F^{\prime}S^{\prime}J^{\prime}$ | $j$ | $FSJ\rightarrow F^{\prime}S^{\prime}J^{\prime}$ |  $f_{ij}^{exp}$ Patra20  |  $f_{ij}^{theor}$ Karr20  |  $f_{ij}^{theor}$ (this work) | $\Delta_{ij}$ | $\Delta_{ij}/\sigma_{c}$
$F=0$ | $014\rightarrow 014$ | $F=1$ | $125\rightarrow 125$ | 178 254.4(9) | 178 246.2(1.8) | 178 245.89(28) | 8.5 | 9.0

(c) $(v=0,L=0)\to(v^{\prime}=0,L^{\prime}=1)$
---
coefficient | $E_{k}^{exp}$ | $E_{k}^{theor}$ |  $E_{k}^{theor}$ (this work) | $\Delta E_{k}$ | $\Delta E_{k}/\sigma_{c}$
$E_{1}$ | 31 984.9(1) |  31 985.76(35) Korobov20  | 31 985.41(12) | -0.5 | -3.3
$E_{6}$ | 8 611.17(5) |  8 611.1(5) Bakalov06  | 8 611.299(18) | -0.13 | -2.4
$E_{7}$ | 1 321.72(4) |  1 321.77(7) Bakalov06  | 1 321.7960(28) | -0.08 | -2.0

Table 4: Comparison between experimental and theoretical hyperfine intervals
(unit: kHz). In parts (a) and (b), hyperfine components are identified by
their label (following the notations of the original publications Patra20 ;
Kortunov21 ) and by the hyperfine states in the lower and upper rovibrational
states (columns 1-4). Experimental values and previous theoretical predictions
of the hyperfine intervals are given in column 5 and 6. Our predictions are
shown in column 7. The deviations
$\Delta_{ij}=f_{ij}^{exp}-\Delta_{ij}^{theor}$ are given in kHz in column 8,
and in units of the combined uncertainty
$\sigma_{c}=([u(f_{ij}^{exp})]^{2}+[u(f_{ij}^{theor})]^{2})^{1/2}$ in the last
column. In part (c), column 2 contains the values of the $E_{1}$, $E_{6}$ and
$E_{7}$ hyperfine coefficients extracted from the experimental data of Ref.
Alighanbari20  using a least-squares adjustment described in Appendix E.
Previous theoretical predictions are given in column 3. Our predictions are
shown in column 4. Deviations $\Delta E_{k}=E_{k}^{exp}-E_{k}^{theor}$ are
given in kHz in column 5, and in units of the combined uncertainty in the last
column.

###  III.2 Comparison with experimental data

A comparison between experiment and theory is shown in Table 4. Regarding the
hyperfine splitting of the $(v=0,L=0)\to(v^{\prime}=1,L^{\prime}=1)$
transition Kortunov21  and of the $(v=0,L^{\prime}=3)\to(v=9,L^{\prime}=3)$
transition Patra20 , the theoretical precision is significantly improved, and
previous conclusions remain essentially valid: reasonable agreement between
theory and experiment is observed in the first case, whereas a clear
discrepancy appears in the second case, becoming even more significant (9
combined standard deviations) due to the reduced theoretical error bar. The
cause of this discrepancy is currently unknown, but it is clear that it has no
relationship with the corrections calculated in the present work, which shift
the theoretical prediction by only -0.3 kHz with respect to the previous
evaluation Karr20 .

We now discuss the rotational transition Alighanbari20 , which represents the
most stringent test of our improved values of the $E_{1}$, $E_{6}$, and
$E_{7}$ coefficients, in view of the higher absolute precision of the
measurements and of their much lower sensitivity to the Fermi coefficients
$E_{4}$ and $E_{5}$ (due to strong cancellation between the lower and upper
states). Comparison between theory and experiment is more involved in this
case, as six hyperfine components of this transition have been measured, from
which fifteen hyperfine intervals can be deduced (the values of which are, of
course, partially redundant). A detailed comparison for these fifteen
intervals is shown in the Appendix D, revealing significant discrepancies for
several lines. These data are however difficult to interpret, and a useful
overview of the results can be obtained by extracting experimental values of
$E_{1}$, $E_{6}$, and $E_{7}$ from a least-squares adjustment of the
experimental data, whose results are shown in the last line of Table 4. The
difference between experimental and theoretical values is above 3 $\sigma_{c}$
for the spin-orbit coefficient, and 2 $\sigma_{c}$ for spin-spin tensor
coefficients. Details of the least-squares adjustment procedure can be found
in the Appendix E.

###  III.3 Discussion

The observed deviations could be due to a problem in the theory, or in the
experiment, or both; let us discuss here the first of these possibilities.

The largest discrepancy in terms of absolute magnitude (8.5 kHz, or a
fractional difference of $4.8\\!\times\\!10^{-5}$) is in the hyperfine
splitting of the $v=0\to 9$ transition. As discussed in Karr20  it points
towards the largest hyperfine coefficients $E_{4}$ and/or $E_{5}$, because the
error that would be required in the other coefficients to explain the
discrepancy is much larger than the expected order of magnitude of any
corrections beyond the Breit-Pauli Hamiltonian. However, a mistake in the
calculation of the electron-proton spin-spin coefficient $E_{4}$ would also
affect at a similar level the corresponding coefficient (denoted by $b_{F}$)
in H${}_{2}^{+}$, where agreement with experiments at a sub-kHz level is
observed Karr20 . This suggests that the problem might come from the electron-
deuteron interaction coefficient $E_{5}$ (note that this hypothesis is not
contradictory with the reasonably good agreement obtained for the $v=0\to 1$
transition frequency, because the latter only weakly depends on $E_{5}$).

What could then be missing in the theory of the electron-deuteron spin-spin
interaction ? The difference between proton and deuteron cases resides in the
nuclear structure and recoil corrections; these have been determined
phenomenologically from the ground-state hyperfine splitting of the hydrogen
and deuterium atoms Karr20  under the assumptions that they are state-
independent, i.e., they can be described by a contact (delta-function) term.
Possible errors induced by this approximation are linked to the existence of
state-dependent corrections. The largest such term is a recoil correction of
order $(Z\alpha)^{2}(m/M)E_{F}$ Bodwin88 , which contributes to the ground-
state hyperfine splitting of the deuterium atom at a level of only 39 Hz. An
error on the electron-deuteron spin-spin interaction at a level of $\sim 10$
kHz thus seems extremely unlikely.

Regarding the rotational transition, the deviations between theory and
experiment could be explained by errors in the spin-orbit ($E_{1}$) and spin-
spin tensor coefficients ($E_{6}$, $E_{7}$) coefficients, as revealed by the
adjustment described in the previous section (see last line of Table 4).
However, any error in the calculation of $E_{1}$ for the $(L=1,v=0)$ would
affect the $(L=1,v=1)$ level at a similar level. If one assumes, for example,
that $E_{1}$ is shifted by -0.5 kHz for both states with respect to our
present theoretical value, this would improve the agreement for the rotational
transition, but increase the deviation from 0.4 kHz to about 0.9 kHz for the
$v=0\to 1$ transition. Moreover, it would affect the value of the
corresponding coefficient (denoted by $c_{e}$) in H${}_{2}^{+}$, which we
calculated in our previous work Haidar22 , and cause significant tension with
rf spectroscopy data. Overall, it appears that available measurements are not
fully consistent with each other.

Apart from errors or missed contributions in theoretical hyperfine
coefficients, another possible cause of deviations between theory and
experiment is the incompleteness of the effective spin Hamiltonian introduced
in Ref. Bakalov06 . Only the terms that appear at the leading order (Breit-
Pauli Hamiltonian) have been included so far, but other spin couplings appear
at the order $m\alpha^{6}$ and higher. The largest one is the spin-spin
contact interaction between proton and deuteron
($\mathbf{I}_{p}\\!\cdot\\!\mathbf{I}_{d}$). This term is already present in
the Breit-Pauli Hamiltonian, but was not included in the effective spin
Hamiltonian because the associated coupling coefficient is proportional to the
proton-deuteron delta function expectation value, which is negligibly small
due to the strong Coulomb repulsion. However, a larger coupling appears at the
$m\alpha^{6}$ order due to the second-order contribution mediated by the
electron. This contribution was studied in Korobov22 , where it was shown that
the coupling constant is on the order of 100 Hz, but shifts the hyperfine
components of the rotational transition by less than 2 Hz due to cancellation
between the lower and upper levels.

In conclusion, we did not identify any contribution that could potentially
have the required order of magnitude to explain the observed discrepancies,
and there is a strong need for additional measurements of the HD+ hyperfine
structure to give new insight on this problem.

## Acknowledgements

The authors thank J. C. J. Koelemeij for helpful discussions.

## References

  * (1) S. Alighanbari, G. S. Giri, F. L. Constantin, V. I. Korobov, and S. Schiller, Precise test of quantum electrodynamics and determination of fundamental constants with HD+ ions, Nature 581, 152 (2020).
  * (2) S. Patra, M. Germann, J.-Ph. Karr, M. Haidar, L. Hilico, V. I. Korobov, F. M. J. Cozijn, K. S. E. Eikema, W. Ubachs, and J. C. J. Koelemeij, Proton-electron mass ratio from laser spectroscopy of HD+ at the part-per-trillion level, Science 369, 1238 (2020).
  * (3) I. V. Kortunov, S. Alighanbari, M. G. Hansen, G. S. Giri, V. I. Korobov, and S. Schiller, Proton-electron mass ratio by high-resolution optical spectroscopy of ion ensembles in the resolved-carrier regime, Nature Phys. 17, 569 (2021).
  * (4) V.I. Korobov, L. Hilico, and J.-Ph. Karr, Fundamental Transitions and Ionization Energies of the Hydrogen Molecular Ions with Few ppt Uncertainty, Phys. Rev. Lett. 118, 233001 (2017).
  * (5) V. I. Korobov and J.-Ph. Karr, Rovibrational spin-averaged transitions in the hydrogen molecular ions, Phys. Rev. A 104, 032806 (2021).
  * (6) M. Germann, S. Patra, J.-Ph. Karr, L. Hilico, V. I. Korobov, E. J. Salumbides, K. S. E. Eikema, W. Ubachs, and J.C.J. Koelemeij, Three-body QED test and fifth-force constraint from vibrations and rotations of HD+, Phys. Rev. Research 3, L022028 (2021).
  * (7) D. Bakalov, V. I. Korobov, and S. Schiller, High-Precision Calculation of the Hyperfine Structure of the HD+ Ion, Phys. Rev. Lett. 97, 243001 (2006).
  * (8) J.-Ph. Karr, M. Haidar, L. Hilico, Z.-X. Zhong, and V. I. Korobov, Higher-order corrections to spin-spin scalar interactions in HD+ and H${}_{2}^{+}$”, Phys. Rev. A 102, 052827 (2020).
  * (9) J. C. J. Koelemeij, Effect of correlated hyperfine theory errors in the determination of rotational and vibrational transition frequencies in HD+, J. Mol. Spectrosc., DOI: 10.1080/00268976.2022.2058637 (2022).
  * (10) R. F. Code and N. F. Ramsey, Molecular-Beam Magnetic Resonance Studies of HD and D2, Phys. Rev. A 4, 1945 (1971).
  * (11) M. Puchalski, J. Komasa, and K. Pachucki, Hyperfine Structure of the First Rotational Level in H2, D2 and HD Molecules and the Deuteron Quadrupole Moment, Phys. Rev. Lett. 125, 253001 (2020).
  * (12) V. I. Korobov, J.-Ph. Karr, M. Haidar, and Z.-X. Zhong, Hyperfine structure in the H${}_{2}^{+}$ and HD+ molecular ions at order $m\alpha^{6}$, Phys. Rev. A 102, 022804 (2020).
  * (13) M. Haidar, V. I. Korobov, L. Hilico, and J.-Ph. Karr, Higher-order corrections to spin-orbit and spin-spin tensor interactions in hydrogen molecular ions: theory and application to H${}_{2}^{+}$, Phys. Rev. A 106, 022816 (2022).
  * (14) Z.-X. Zhong, X. Tong, Z.-C. Yan, and T.-Y. Shi, High-precision spectroscopy of hydrogen molecular ions, Chin. Phys. B 24, 053102 (2015).
  * (15) A. J. Thakkar and V. H. Smith, Jr., Compact and accurate integral-transform wave functions. I. The 1${}^{1}S$ state of the helium-like ions from H- through Mg10+, Phys. Rev. A 15, 1 (1977).
  * (16) V. I. Korobov, Coulomb three-body bound-state problem: Variational calculations of nonrelativistic energies, Phys. Rev. A 61, 064503 (2000).
  * (17) G. T. Bodwin and D. R. Yennie, Some recoil corrections to the hydrogen hyperfine splitting, Phys. Rev. D 37, 498 (1988).
  * (18) V. I. Korobov, Precision spectroscopy of the hydrogen molecular ions: present status of theory and experiment, Phys. Part. Nucl. 53, 787 (2022).
  * (19) P. J. Mohr and B. N. Taylor, CODATA recommended values of the fundamental physical constants: 1998, Rev. Mod. Phys. 72, 351 (2000).

##  Appendix A Other coefficients of the effective spin Hamiltonian

We give in Table 5 the values of all the coefficients of the HD+ effective
spin Hamiltonian (see Eq. (3) of Bakalov06 ) in the rovibrational states
considered in Sec. III (see Table 4).

$(v,L)$ |  $E_{2}$ |  $E_{3}$ |  $E_{4}$ |  $E_{5}$ |  $E_{8}$ |  $E_{9}$
---|---|---|---|---|---|---
$(0,0)$ |  |  |  925 394.159(860) |  142 287.556(84) |  |
$(0,1)$ |  -31.345(2) |  -4.809 |  924 567.718(859) |  142 160.670(84) |  -3.057 |  5.660
$(1,1)$ |  -30.463(2) |  -4.664 |  903 366.501(839) |  138 910.266(82) |  -2.945 |  5.653
$(0,3)$ |  -30.832(2) |  -4.733 |  920 479.981(855) |  141 533.075(83) |  -0.335 |  0.612
$(9,3)$ |  -21.304(1) |  -3.225 |  775 706.122(721) |  119 431.933(73) |  -0.219 |  0.501
Table 5: Hyperfine coefficients for a few rovibrational states of HD+ (in
kHz). The values of $E_{4}$, $E_{5}$ (resp. $E_{2}$, $E_{3}$, $E_{8}$) are
taken from Karr20  (resp. Bakalov06 ). The value of $E_{9}$ has been updated
with the latest determination of the deuteron’s quadrupole moment Puchalski20
. Uncertainties are discussed in Sec. III.1 of the main text.

##  Appendix B Dependence of hyperfine shifts on the hyperfine coefficients

The values of the derivatives $\gamma_{i,k}$ (as defined in Eq. (1)), for the
hyperfine states probed in spectroscopy of the $(L=3,v=0)\to(L=3,v=9)$
transition Patra20 , are given in Table 6. These quantities are required to
estimate the uncertainty of the theoretical hyperfine interval.

$(L,v,F,S,J)$ |  $\gamma_{1,v}$ |  $\gamma_{2,v}$ |  $\gamma_{3,v}$ |  $\gamma_{4,v}$ |  $\gamma_{5,v}$ |  $\gamma_{6,v}$ |  $\gamma_{7,v}$ |  $\gamma_{8,v}$ |  $\gamma_{9,v}$
---|---|---|---|---|---|---|---|---|---
(3,0,0,1,4) |  -0.414 |  0.424 |  2.986 |  -0.726 |  -0.215 |  0.274 |  3.079 |  -3.328 |  -7.400
(3,9,0,1,4) |  -0.359 |  0.372 |  2.984 |  -0.730 |  -0.201 |  0.193 |  2.474 |  -2.722 |  -7.385
(3,0,1,2,5) |  1.50 |  1.50 |  3.00 |  0.250 |  0.500 |  -7.50 |  -15.00 |  -15.00 |  -7.50
(3,9,1,2,5) |  1.50 |  1.50 |  3.00 |  0.250 |  0.500 |  -7.50 |  -15.00 |  -15.00 |  -7.50
Table 6: Sensitivity coefficients for the hyperfine states involved in the
$v=0-9$ transition in HD+. Here $\gamma_{n,v}=\partial
E_{hfs}(v,L,F,S,J)/\partial E_{n}$ . Note that for the $F=1$ level, the
sensitivities are the same in $v=0$ and $v=9$, which is due to $(F=1,S=2,J=5)$
being a “pure” state.

##  Appendix C Correlations between the hyperfine coefficients

In this Appendix, we present and justify our hypotheses regarding correlations
between theoretical errors in the hyperfine coefficients, which are summarized
in Table 7.

Let us first discuss the correlations between errors on a given coefficient
$E_{k}$ in different rovibrational states. In the case of the largest
hyperfine coefficients $E_{4}$ and $E_{5}$, related to the Fermi contact
interaction between the electron and both nuclei, uncertainties and
correlations can be controlled very well, thanks to the fact that these
interactions have been studied in great depth in the context of the hydrogen
and deuterium ground-state hyperfine splitting. This has allowed us to
estimate the unevaluated QED terms from their value in the hydrogen
(deuterium) atom ground state using the LCAO approximation. The uncertainties
of $E_{4}$ and $E_{5}$ are conservatively taken as equal to the full estimated
term Karr20 . Then, for the evaluation of uncertainties of hyperfine shifts in
rovibrational transitions, we assume perfect correlations between different
rovibrational states, i.e., the uncertainty of a transition frequency is
obtained by subtracting the uncertainties due to the upper and lower states.
We have tested this assumption by applying it to lower-order terms that have
been actually calculated by us, and found that the uncertainty estimated in
this way matches well the magnitude of the actual correction.

For other hyperfine coefficients, we have used cruder uncertainty estimates
based on the expected order of magnitude of the largest unevaluated QED terms.
For example, for coefficients calculated at the Breit-Pauli, we estimate the
relative uncertainty to be $\alpha^{2}$. In this case, it is harder to
correctly evaluate the degree of correlation between theory errors, but all
correlations are expected to be positive. In this work, we have assumed that
the theory errors of the lower and upper levels are uncorrelated, which
provides an upper limit of uncertainties on transition frequencies.

Let us now discuss correlations between errors in different coefficients
$E_{k}$, $E_{l}$ for the same rovibrational state. Such correlations occur if
the uncalculated QED terms that limit the theoretical precision are of the
same nature. This is the case for the $E_{2}$ and $E_{3}$ coefficients that
correspond to the proton ($E_{2}$) and deuteron ($E_{3}$) spin-orbit
interactions. Both coefficients have been calculated in the framework of the
Breit-Pauli Hamiltonien Bakalov06 , and their theoretical errors are
associated with $(Z\alpha^{2})$-order relativistic corrections. These
corrections will be described by the same effective potentials, the only
difference being the substitution between proton and deuteron. Similarly,
errors in the $E_{6}$ and $E_{7}$ coefficients that correspond to the
electron-proton ($E_{6}$) and electron-deuteron ($E_{7}$) spin-spin tensor
interactions are caused by the same QED term, that is the nonlogarithmic
$m\alpha^{7}$-order radiative corrections (see Sec. II). We thus assume
perfect correlations between theory errors of $E_{2}$ and $E_{3}$, as well as
between $E_{6}$ and $E_{7}$.

The case of the electron-proton ($E_{4}$) and electron-deuteron ($E_{5}$)
Fermi contact interactions is slightly more involved. Indeed, errors in these
coefficients originate from both nonrecoil QED corrections, which will have
the same expressions apart from substitution of proton and deuteron, and from
recoil corrections, which depend on the nucleus Karr20 . We assume that errors
associated with nonrecoil (resp. recoil) terms are fully correlated (resp.
uncorrelated), which yields a correlation coefficient.

$r(E_{4},E_{5})=\frac{u_{non-rec}(E_{4})u_{non-
rec}(E_{5})}{u_{tot}(E_{4})u_{tot}(E_{5})}=0.4016.$ (2)

Here, $u_{non-rec}(E_{k})$ is the uncertainty of $E_{k}$ due to nonrecoil QED
corrections, and $u_{tot}(E_{k})$ its total uncertainty. These uncertainties
are estimated to $u_{non-rec}(E_{k})=0.47\times 10^{-6}E_{k}^{(F)}$,
$u_{tot}(E_{4})=0.93\times 10^{-6}E_{4}^{(F)}$, and $u_{tot}(E_{5})=0.59\times
10^{-6}E_{5}^{(F)}$, where $E_{k}^{(F)}$ is the value of $E_{k}$ at the
leading-order ($m\alpha^{4}$), see Karr20  for details.

First coefficient | Second coefficient | Correlation coefficient
---|---|---
$E_{k}(v,L)$ | $E_{k}(v^{\prime},L^{\prime})$ |  $1$ if $k=4$ or $5$; $0$ otherwise.
$E_{k}(v,L)$ | $E_{l}(v,L)$ |  |  $1$ if $(k,l)=(2,3)$, $(3,2)$, $(6,7)$, or $(7,6)$;
---
$0.4016$ if $(k,l)=(4,5)$ or $(5,4)$; $0$ otherwise.
$E_{k}(v,L)$ | $E_{l}(v^{\prime},L^{\prime})$ |  |  $0.4016$ if $(k,l)=(4,5)$ or $(5,4)$; $0$ otherwise.
---
Table 7: Summary of our assumptions on correlations between errors in the
theoretical hyperfine coefficients.

##  Appendix D Detailed comparison between theory and experiment for the
hyperfine splitting of the rotational transition

A full comparison between the experimental results of Ref. Alighanbari20  and
our theoretical predictions is shown in Table 8.

$(v=0,L=0)\to(v^{\prime}=0,L^{\prime}=1)$
---
$i$ | $FSJ\rightarrow F^{\prime}S^{\prime}J^{\prime}$ | $j$ | $FSJ\rightarrow F^{\prime}S^{\prime}J^{\prime}$ |  $f_{ij}^{exp}$ Alighanbari20  |  $f_{ij}^{theor}$ Alighanbari20  |  $f_{ij}^{theor}$ (this work) | $\Delta_{ij}$ | $\Delta_{ij}/\sigma_{c}$
12 | $122\rightarrow 121$ | 14 | $100\rightarrow 101$ | 24 134.211(75) | 24 134.5(1.1) | 24 134.465(23) | -0.254 | -3.2
12 | $122\rightarrow 121$ | 16 | $011\rightarrow 012$ | 31 073.752(43) | 31 073.7(1.2) | 31 074.102(55) | -0.350 | -5.0
12 | $122\rightarrow 121$ | 19 | $122\rightarrow 123$ | 43 283.419(54) | 43 283.4(2.2) | 43 284.10(12) | -0.677 | -5.0
12 | $122\rightarrow 121$ | 20 | $122\rightarrow 122$ | 44 944.338(72) | 44 944.6(1.9) | 44 945.289(64) | -0.951 | -9.9
12 | $122\rightarrow 121$ | 21 | $111\rightarrow 112$ | 44 996.486(61) | 44 996.4(1.9) | 44 997.14(11) | -0.652 | -5.4
14 | $100\rightarrow 101$ | 16 | $011\rightarrow 012$ | 6 939.541(66) | 6 939.2(1.2) | 6 939.636(42) | -0.095 | -1.2
14 | $100\rightarrow 101$ | 19 | $122\rightarrow 123$ | 19 149.208(74) | 19 148.8(2.0) | 19 149.63(11) | -0.423 | -3.2
14 | $100\rightarrow 101$ | 20 | $122\rightarrow 122$ | 20 810.127(88) | 20 810.1(1.8) | 20 810.823(62) | -0.696 | -6.5
14 | $100\rightarrow 101$ | 21 | $111\rightarrow 112$ | 20 862.275(79) | 20 861.9(2.0) | 20 862.673(91) | -0.398 | -3.3
16 | $011\rightarrow 012$ | 19 | $122\rightarrow 123$ | 12 209.667(41) | 12 209.7(1.5) | 12 209.994(72) | -0.327 | -4.0
16 | $011\rightarrow 012$ | 20 | $122\rightarrow 122$ | 13 870.586(62) | 13 870.9(1.4) | 13 871.187(43) | -0.601 | -8.0
16 | $011\rightarrow 012$ | 21 | $111\rightarrow 112$ | 13 922.734(49) | 13 922.7(1.0) | 13 923.037(51) | -0.303 | -4.3
19 | $122\rightarrow 123$ | 20 | $122\rightarrow 122$ | 1 660.919(70) | 1 661.2(2.6) | 1 661.193(99) | -0.274 | -2.2
19 | $122\rightarrow 123$ | 21 | $111\rightarrow 112$ | 1 713.067(59) | 1 713.0(0.7) | 1 713.042(25) | 0.025 | 0.4
20 | $122\rightarrow 122$ | 21 | $111\rightarrow 112$ | 52.148(76) | 51.8(1.9) | 51.850(75) | 0.298 | 2.8
Table 8: Comparison between experimental and theoretical hyperfine intervals
(unit: kHz). All definitions are identical to those of Table 4.

##  Appendix E Least-squares adjustment of the spin-orbit and spin-spin tensor
coefficients in the $(v=0,L=1)$ state from experimental data

The frequencies of the hyperfine components of the rotational transition
(measured in Alighanbari20 ) can be written as a sum of the “spin-averaged”
transition frequency, $f_{SA}$ and a hyperfine shift. In Alighanbari20  (see
also Koelemeij22 ), a value of $f_{SA}$ is extracted from the experimental
data. Here, our goal is instead to extract experimental values of three
hyperfine coefficients, $E_{1}$, $E_{6}$, and $E_{7}$. We do this in a way
that is completely independent of the value of $f_{SA}$, which we describe in
the following.

From the six measured transition frequencies, we deduce five independent
hyperfine intervals by choosing one of the lines, $f_{i}$, as reference and
computing the differences $f_{ij}=f_{j}-f_{i}$ for $j\neq i$. We then perform
a least-square adjustment taking into account experimental and theoretical
uncertainties. The latter are due to uncertainties in hyperfine coefficients
(other than $E_{1}$, $E_{6}$, $E_{7}$), and are accounted for by additive
corrections $\delta(E_{k}(v,L))$ that are treated as additional adjusted
parameters, as done in the CODATA adjustments Mohr00 . Correlations between
theoretical errors are taken as described in the Appendix C. Experimental
uncertainties are assumed to be uncorrelated. For illustration, the input data
used in the adjustment are shown in Table 9 for the case where the line
labeled 12 is chosen as reference line.

$f_{14}-f_{12}=24\,134.211(75)$ kHz
---
$f_{16}-f_{12}=31\,073.752(43)$ kHz
$f_{19}-f_{12}=43\,283.419(54)$ kHz
$f_{20}-f_{12}=44\,944.338(72)$ kHz
$f_{21}-f_{12}=44\,996.486(61)$ kHz
$\delta(E_{4}(v=0,L=0))=0.000(860)$ kHz
$\delta(E_{5}(v=0,L=0))=0.000(84)$ kHz
$\delta(E_{2}(v=0,L=1))=0.0000(17)$ kHz
$\delta(E_{8}(v=0,L=1))=0.00000(16)$ kHz
$\delta(E_{9}(v=0,L=1))=0.00000(55)$ kHz
Table 9: Input data for the least-squares adjustment of the $E_{1}$, $E_{6}$,
and $E_{7}$ coefficients, where the line 12 is chosen as reference line.

We repeated this procedure for all possible choices of a reference line, and
found that the adjusted values slightly depend on the chosen line. In our
final results (reported in Table 4), we increased the error bars to make the
values compatible with all possible choices. We also checked that the
dependence of our results on the hypotheses made on correlations between
theoretical errors (see Appendix C) is negligibly small.

It is worth stressing that the adjustments give satisfactory results (with
residuals below $1\sigma$), indicating that experimental data are well
explained by the effective spin Hamiltonian introduced in Bakalov06  and in
good agreement with the theoretical values of the six other coefficients. The
only exception is the hyperfine interval involving line 16 (such as
$f_{16}-f_{12}$ in Table 9), which has $2-3\sigma$ residuals in all cases.
This might be an indication of some experimental problem for this specific
line.

##  Appendix F Numerical results

We give here the numerical values of all contributions of orders $m\alpha^{6}$
and $m\alpha^{7}\ln(\alpha)$ to the spin-orbit (Table 10) and spin-spin tensor
(Tables 11 and 12) coefficients for a range of rovibrational states
($v=0-9,L=1-4$) of HD+.

$(L,v)$ | $\displaystyle E_{1}^{(BP)}$ | $\displaystyle\mathcal{U}_{Y_{1}}$ | $\displaystyle\mathcal{U}_{W}$ | $\displaystyle\mathcal{U}_{CM}$ | $\displaystyle\Delta E_{so\hbox{-}H_{B}}$ | $\displaystyle\Delta E_{so\hbox{-}so}^{(1)}$ | $\displaystyle\Delta E_{so\hbox{-}ret}$ | $\displaystyle\Delta E_{1}^{(6)}$ | $\displaystyle\mathcal{U}_{Y_{2}}$ | $\displaystyle\mathcal{U}_{q^{2}}$ | $\displaystyle\Delta E_{so\hbox{-}H_{(5)}}^{\ln}$ | $\displaystyle\Delta E_{1}^{(7\ln)}$ |  $\displaystyle E_{1}$ (this work)
---|---|---|---|---|---|---|---|---|---|---|---|---|---
(1,0) | 31984.645 | 1.170 | -2.736 | 0.021 | 2.088 | 0.312 | 0.257 | 1.112 | -0.026 | 0.045 | -0.367 | -0.348 | 31985.41(12)
(1,1) | 30280.027 | 1.110 | -2.631 | 0.027 | 1.986 | 0.299 | 0.243 | 1.035 | -0.025 | 0.044 | -0.345 | -0.326 | 30280.74(11)
(1,2) | 28644.571 | 1.051 | -2.525 | 0.032 | 1.886 | 0.286 | 0.230 | 0.960 | -0.024 | 0.042 | -0.325 | -0.306 | 28645.23(10)
(1,3) | 27070.453 | 0.993 | -2.419 | 0.036 | 1.789 | 0.269 | 0.217 | 0.886 | -0.022 | 0.040 | -0.305 | -0.288 | 27071.05(10)
(1,4) | 25550.256 | 0.941 | -2.313 | 0.039 | 1.689 | 0.257 | 0.205 | 0.819 | -0.021 | 0.038 | -0.286 | -0.269 | 25550.81(9)
(1,5) | 24076.867 | 0.893 | -2.205 | 0.042 | 1.592 | 0.233 | 0.193 | 0.747 | -0.020 | 0.036 | -0.268 | -0.252 | 24077.36(8)
(1,6) | 22643.449 | 0.834 | -2.097 | 0.044 | 1.498 | 0.219 | 0.181 | 0.679 | -0.019 | 0.035 | -0.251 | -0.235 | 22643.89(8)
(1,7) | 21243.261 | 0.781 | -1.988 | 0.045 | 1.403 | 0.207 | 0.170 | 0.617 | -0.018 | 0.033 | -0.235 | -0.220 | 21243.66(7)
(1,8) | 19869.606 | 0.727 | -1.877 | 0.045 | 1.311 | 0.193 | 0.159 | 0.558 | -0.017 | 0.031 | -0.219 | -0.204 | 19869.96(7)
(1,9) | 18515.725 | 0.669 | -1.764 | 0.045 | 1.217 | 0.179 | 0.147 | 0.493 | -0.016 | 0.029 | -0.203 | -0.190 | 18516.03(6)
(2,0) | 31840.910 | 1.164 | -2.719 | 0.020 | 2.069 | 0.311 | 0.256 | 1.101 | -0.026 | 0.045 | -0.364 | -0.345 | 31841.67(12)
(2,1) | 30142.781 | 1.105 | -2.614 | 0.026 | 1.969 | 0.296 | 0.243 | 1.024 | -0.025 | 0.043 | -0.343 | -0.324 | 30143.48(11)
(2,2) | 28513.394 | 1.047 | -2.509 | 0.031 | 1.868 | 0.281 | 0.229 | 0.948 | -0.023 | 0.042 | -0.322 | -0.304 | 28514.04(10)
(2,3) | 26944.952 | 0.992 | -2.404 | 0.035 | 1.771 | 0.266 | 0.216 | 0.877 | -0.022 | 0.040 | -0.303 | -0.285 | 26945.54(10)
(2,4) | 25430.065 | 0.937 | -2.298 | 0.039 | 1.673 | 0.252 | 0.204 | 0.807 | -0.021 | 0.038 | -0.284 | -0.267 | 25430.61(9)
(2,5) | 23961.639 | 0.889 | -2.191 | 0.041 | 1.558 | 0.239 | 0.192 | 0.729 | -0.020 | 0.036 | -0.266 | -0.250 | 23962.12(8)
(2,6) | 22532.852 | 0.836 | -2.083 | 0.043 | 1.463 | 0.224 | 0.180 | 0.664 | -0.019 | 0.034 | -0.249 | -0.233 | 22533.28(8)
(2,7) | 21136.966 | 0.785 | -1.974 | 0.044 | 1.369 | 0.211 | 0.169 | 0.604 | -0.018 | 0.032 | -0.232 | -0.217 | 21137.35(7)
(2,8) | 19767.295 | 0.733 | -1.863 | 0.044 | 1.270 | 0.197 | 0.158 | 0.539 | -0.017 | 0.031 | -0.216 | -0.202 | 19767.63(7)
(2,9) | 18417.096 | 0.669 | -1.751 | 0.044 | 1.184 | 0.185 | 0.146 | 0.477 | -0.015 | 0.029 | -0.201 | -0.187 | 18417.39(6)
(3,0) | 31627.352 | 1.156 | -2.694 | 0.019 | 2.042 | 0.308 | 0.255 | 1.086 | -0.026 | 0.045 | -0.360 | -0.341 | 31628.10(11)
(3,1) | 29938.871 | 1.097 | -2.590 | 0.025 | 1.942 | 0.293 | 0.241 | 1.008 | -0.025 | 0.043 | -0.339 | -0.320 | 29939.56(11)
(3,2) | 28318.505 | 1.040 | -2.486 | 0.030 | 1.843 | 0.279 | 0.228 | 0.934 | -0.023 | 0.041 | -0.319 | -0.301 | 28319.14(10)
(3,3) | 26758.499 | 0.985 | -2.381 | 0.034 | 1.745 | 0.264 | 0.215 | 0.863 | -0.022 | 0.039 | -0.299 | -0.282 | 26759.08(9)
(3,4) | 25251.499 | 0.930 | -2.275 | 0.037 | 1.649 | 0.250 | 0.203 | 0.795 | -0.021 | 0.038 | -0.281 | -0.264 | 25252.03(9)
(3,5) | 23790.452 | 0.880 | -2.169 | 0.040 | 1.549 | 0.237 | 0.191 | 0.728 | -0.020 | 0.036 | -0.263 | -0.247 | 23790.93(8)
(3,6) | 22368.539 | 0.830 | -2.062 | 0.042 | 1.455 | 0.223 | 0.179 | 0.666 | -0.019 | 0.034 | -0.246 | -0.231 | 22368.98(8)
(3,7) | 20979.051 | 0.780 | -1.953 | 0.043 | 1.356 | 0.209 | 0.168 | 0.602 | -0.018 | 0.032 | -0.229 | -0.215 | 20979.44(7)
(3,8) | 19615.310 | 0.731 | -1.844 | 0.043 | 1.256 | 0.195 | 0.157 | 0.540 | -0.016 | 0.030 | -0.213 | -0.199 | 19615.65(7)
(3,9) | 18270.560 | 0.680 | -1.732 | 0.043 | 1.158 | 0.182 | 0.146 | 0.478 | -0.015 | 0.028 | -0.198 | -0.185 | 18270.85(6)
(4,0) | 31346.356 | 1.146 | -2.662 | 0.018 | 2.007 | 0.305 | 0.253 | 1.067 | -0.026 | 0.044 | -0.355 | -0.336 | 31347.09(11)
(4,1) | 29670.579 | 1.087 | -2.558 | 0.023 | 1.908 | 0.290 | 0.239 | 0.990 | -0.024 | 0.042 | -0.334 | -0.316 | 29671.25(11)
(4,2) | 28062.088 | 1.031 | -2.455 | 0.028 | 1.810 | 0.276 | 0.226 | 0.916 | -0.023 | 0.041 | -0.314 | -0.296 | 28062.71(10)
(4,3) | 26513.184 | 0.975 | -2.351 | 0.032 | 1.715 | 0.261 | 0.213 | 0.846 | -0.022 | 0.039 | -0.295 | -0.278 | 26513.96(9)
(4,4) | 25016.561 | 0.923 | -2.246 | 0.036 | 1.618 | 0.248 | 0.201 | 0.779 | -0.021 | 0.037 | -0.276 | -0.260 | 25017.25(9)
(4,5) | 23565.204 | 0.873 | -2.141 | 0.038 | 1.518 | 0.234 | 0.189 | 0.712 | -0.020 | 0.035 | -0.259 | -0.243 | 23565.83(8)
(4,6) | 22152.327 | 0.822 | -2.034 | 0.040 | 1.427 | 0.220 | 0.178 | 0.653 | -0.018 | 0.034 | -0.242 | -0.227 | 22152.87(8)
(4,7) | 20771.242 | 0.772 | -1.927 | 0.041 | 1.328 | 0.207 | 0.166 | 0.588 | -0.017 | 0.032 | -0.226 | -0.211 | 20771.71(7)
(4,8) | 19415.281 | 0.722 | -1.818 | 0.041 | 1.232 | 0.193 | 0.155 | 0.527 | -0.016 | 0.030 | -0.210 | -0.196 | 19415.67(7)
(4,9) | 18077.684 | 0.674 | -1.707 | 0.041 | 1.138 | 0.180 | 0.144 | 0.471 | -0.015 | 0.028 | -0.194 | -0.181 | 18078.00(6)
Table 10: Numerical results for the spin-orbit coefficient $E_{1}$ in HD+ for
the range of rovibrational states $(L=1-4)$ and $(v=0-9)$ (in kHz). All
definitions are identical to those given in Table 1.

$(L,v)$ | $\displaystyle E_{6}^{(BP)}$ | $\displaystyle\mathcal{U}_{2d}^{(2)}$ | $\displaystyle\mathcal{U}_{5b}^{(2)}$ | $\displaystyle\Delta E_{ss}^{(2)}$ | $\displaystyle\Delta E_{6}^{(6)}$ | $\displaystyle\Delta E_{6}^{(7)}$ |  $\displaystyle E_{6}$ (this work)
---|---|---|---|---|---|---|---
(1,0) | 8611.112 | -0.806 | 0.093 | 0.955 | 0.242 | -0.055 | 8611.299(18)
(1,1) | 8136.686 | -0.770 | 0.090 | 0.905 | 0.225 | -0.052 | 8136.859(17)
(1,2) | 7682.610 | -0.735 | 0.087 | 0.857 | 0.210 | -0.049 | 7682.771(16)
(1,3) | 7246.701 | -0.700 | 0.084 | 0.807 | 0.191 | -0.047 | 7246.846(16)
(1,4) | 6826.898 | -0.666 | 0.082 | 0.762 | 0.178 | -0.044 | 6827.031(15)
(1,5) | 6421.232 | -0.632 | 0.079 | 0.726 | 0.172 | -0.041 | 6421.364(14)
(1,6) | 6027.809 | -0.599 | 0.076 | 0.678 | 0.154 | -0.039 | 6027.925(13)
(1,7) | 5644.779 | -0.565 | 0.071 | 0.633 | 0.140 | -0.037 | 5644.883(12)
(1,8) | 5270.312 | -0.531 | 0.068 | 0.590 | 0.127 | -0.034 | 5270.405(11)
(1,9) | 4902.573 | -0.498 | 0.065 | 0.541 | 0.108 | -0.032 | 4902.649(11)
(2,0) | 2043.1452 | -0.1908 | 0.0219 | 0.2259 | 0.0569 | -0.0130 | 2043.1891(43)
(2,1) | 1930.5230 | -0.1823 | 0.0213 | 0.2144 | 0.0534 | -0.0124 | 1930.5640(41)
(2,2) | 1822.7196 | -0.1740 | 0.0207 | 0.2032 | 0.0499 | -0.0117 | 1822.7578(39)
(2,3) | 1719.2168 | -0.1658 | 0.0200 | 0.1911 | 0.0454 | -0.0110 | 1719.2512(37)
(2,4) | 1619.5252 | -0.1576 | 0.0193 | 0.1806 | 0.0422 | -0.0105 | 1619.5570(35)
(2,5) | 1523.1776 | -0.1496 | 0.0185 | 0.1712 | 0.0401 | -0.0098 | 1523.2080(33)
(2,6) | 1429.7235 | -0.1416 | 0.0178 | 0.1577 | 0.0339 | -0.0092 | 1429.7482(31)
(2,7) | 1338.7229 | -0.1337 | 0.0170 | 0.1471 | 0.0304 | -0.0086 | 1338.7446(39)
(2,8) | 1249.7406 | -0.1257 | 0.0162 | 0.1353 | 0.0257 | -0.0081 | 1249.7582(27)
(2,9) | 1162.3394 | -0.1177 | 0.0153 | 0.1257 | 0.0233 | -0.0075 | 1162.3552(25)
(3,0) | 948.5222 | -0.0883 | 0.0101 | 0.1042 | 0.0260 | -0.0060 | 948.5421(20)
(3,1) | 896.1987 | -0.0843 | 0.0098 | 0.0992 | 0.0247 | -0.0057 | 896.2176(19)
(3,2) | 846.1054 | -0.0805 | 0.0095 | 0.0936 | 0.0227 | -0.0054 | 846.1227(18)
(3,3) | 798.0019 | -0.0767 | 0.0092 | 0.0880 | 0.0206 | -0.0051 | 798.0174(17)
(3,4) | 751.6607 | -0.0729 | 0.0089 | 0.0833 | 0.0193 | -0.0048 | 751.6752(16)
(3,5) | 706.8644 | -0.0692 | 0.0086 | 0.0793 | 0.0187 | -0.0045 | 706.8786(15)
(3,6) | 663.4035 | -0.0655 | 0.0082 | 0.0737 | 0.0165 | -0.0042 | 663.4158(14)
(3,7) | 621.0731 | -0.0618 | 0.0079 | 0.0690 | 0.0151 | -0.0040 | 621.0842(13)
(3,8) | 579.6702 | -0.0581 | 0.0075 | 0.0645 | 0.0138 | -0.0037 | 579.6803(12)
(3,9) | 538.9906 | -0.0544 | 0.0070 | 0.0593 | 0.0119 | -0.0035 | 538.9991(12)
(4,0) | 550.5184 | -0.0510 | 0.0058 | 0.0600 | 0.0148 | -0.0035 | 550.5298(12)
(4,1) | 520.1196 | -0.0487 | 0.0056 | 0.0570 | 0.0139 | -0.0033 | 520.1302(11)
(4,2) | 491.0099 | -0.0465 | 0.0055 | 0.0536 | 0.0126 | -0.0031 | 491.0194(10)
(4,3) | 463.0496 | -0.0443 | 0.0053 | 0.0508 | 0.0118 | -0.0029 | 463.0585(10)
(4,4) | 436.1067 | -0.0421 | 0.0052 | 0.0480 | 0.0110 | -0.0028 | 436.1150(9)
(4,5) | 410.0548 | -0.0399 | 0.0049 | 0.0446 | 0.0096 | -0.0026 | 410.0618(9)
(4,6) | 384.7717 | -0.0378 | 0.0047 | 0.0424 | 0.0094 | -0.0024 | 384.7786(8)
(4,7) | 360.1380 | -0.0356 | 0.0045 | 0.0398 | 0.0086 | -0.0023 | 360.1443(8)
(4,8) | 336.0351 | -0.0335 | 0.0043 | 0.0370 | 0.0078 | -0.0022 | 336.0407(7)
(4,9) | 312.3437 | -0.0313 | 0.0040 | 0.0335 | 0.0063 | -0.0020 | 312.3480(7)
Table 11: Numerical results for the electron-proton spin-spin tensor
coefficient $E_{6}$ in HD+ for the range of rovibrational states $(L=1-4)$ and
$(v=0-9)$ (in kHz). All definitions are identical to those given in Table 2.

$(L,v)$ | $\displaystyle E_{7}^{(BP)}$ | $\displaystyle\mathcal{U}_{2d}^{(2)}$ | $\displaystyle\mathcal{U}_{5b}^{(2)}$ | $\displaystyle\Delta E_{ss}^{(2)}$ | $\displaystyle\Delta E_{7}^{(6)}$ | $\displaystyle\Delta E_{7}^{(7)}$ |  $\displaystyle E_{7}$ (this work)
---|---|---|---|---|---|---|---
(1,0) | 1321.7673 | -0.1237 | 0.0142 | 0.1467 | 0.0371 | -0.0085 | 1321.7960(28)
(1,1) | 1248.9359 | -0.1182 | 0.0138 | 0.1390 | 0.0346 | -0.0080 | 1248.9624(27)
(1,2) | 1179.2284 | -0.1128 | 0.0135 | 0.1316 | 0.0322 | -0.0076 | 1179.2530(25)
(1,3) | 1112.3096 | -0.1075 | 0.0130 | 0.1239 | 0.0293 | -0.0072 | 1112.3318(24)
(1,4) | 1047.8631 | -0.1023 | 0.0125 | 0.1170 | 0.0273 | -0.0068 | 1047.8836(23)
(1,5) | 985.5868 | -0.0971 | 0.0121 | 0.1115 | 0.0264 | -0.0063 | 985.6069(21)
(1,6) | 925.1895 | -0.0919 | 0.0115 | 0.1040 | 0.0237 | -0.0060 | 925.2072(20)
(1,7) | 866.3873 | -0.0867 | 0.0111 | 0.0972 | 0.0215 | -0.0056 | 866.4032(19)
(1,8) | 808.8992 | -0.0816 | 0.0105 | 0.0906 | 0.0196 | -0.0053 | 808.9135(18)
(1,9) | 752.4434 | -0.0764 | 0.0099 | 0.0830 | 0.0166 | -0.0049 | 752.4552(16)
(2,0) | 696.7322 | -0.0293 | 0.0033 | 0.0347 | 0.0087 | -0.0020 | 313.6203(7)
(2,1) | 296.3244 | -0.0280 | 0.0032 | 0.0329 | 0.0082 | -0.0019 | 296.3307(6)
(2,2) | 279.7750 | -0.0267 | 0.0031 | 0.0312 | 0.0077 | -0.0018 | 279.7808(6)
(2,3) | 263.8857 | -0.0254 | 0.0030 | 0.0293 | 0.0070 | -0.0017 | 263.8910(6)
(2,4) | 248.5814 | -0.0242 | 0.0029 | 0.0277 | 0.0065 | -0.0016 | 248.5863(5)
(2,5) | 233.7905 | -0.0230 | 0.0028 | 0.0263 | 0.0062 | -0.0015 | 233.7951(5)
(2,6) | 219.4436 | -0.0217 | 0.0027 | 0.0242 | 0.0052 | -0.0014 | 219.4474(5)
(2,7) | 205.4734 | -0.0205 | 0.0026 | 0.0226 | 0.0047 | -0.0013 | 205.4767(4)
(2,8) | 191.8128 | -0.0193 | 0.0024 | 0.0208 | 0.0040 | -0.0012 | 191.8155(4)
(2,9) | 178.3949 | -0.0181 | 0.0023 | 0.0193 | 0.0036 | -0.0011 | 178.3973(4)
(3,0) | 145.5939 | -0.0136 | 0.0016 | 0.0160 | 0.0040 | -0.0009 | 145.5969(3)
(3,1) | 137.5614 | -0.0129 | 0.0016 | 0.0152 | 0.0038 | -0.0009 | 137.5643(3)
(3,2) | 129.8713 | -0.0124 | 0.0014 | 0.0144 | 0.0035 | -0.0008 | 129.8740(3)
(3,3) | 122.4867 | -0.0118 | 0.0014 | 0.0135 | 0.0032 | -0.0008 | 122.4891(3)
(3,4) | 115.3726 | -0.0112 | 0.0013 | 0.0128 | 0.0030 | -0.0007 | 115.3748(2)
(3,5) | 108.4956 | -0.0106 | 0.0013 | 0.0122 | 0.0029 | -0.0007 | 108.4978(2)
(3,6) | 101.8236 | -0.0100 | 0.0013 | 0.0113 | 0.0025 | -0.0006 | 101.8255(2)
(3,7) | 95.3251 | -0.0095 | 0.0012 | 0.0106 | 0.0023 | -0.0006 | 95.3268(2)
(3,8) | 88.9689 | -0.0089 | 0.0012 | 0.0099 | 0.0021 | -0.0006 | 88.9705(2)
(3,9) | 82.7237 | -0.0083 | 0.0010 | 0.0091 | 0.0018 | -0.0005 | 82.7250(2)
(4,0) | 84.5020 | -0.0078 | 0.0009 | 0.0092 | 0.0023 | -0.0005 | 84.5038(2)
(4,1) | 79.8354 | -0.0075 | 0.0009 | 0.0088 | 0.0021 | -0.0005 | 79.8370(2)
(4,2) | 75.3666 | -0.0071 | 0.0008 | 0.0082 | 0.0019 | -0.0005 | 75.3680(2)
(4,3) | 71.0742 | -0.0068 | 0.0008 | 0.0078 | 0.0018 | -0.0005 | 71.0756(2)
(4,4) | 66.9381 | -0.0065 | 0.0007 | 0.0074 | 0.0017 | -0.0004 | 66.9393(1)
(4,5) | 62.9387 | -0.0061 | 0.0007 | 0.0068 | 0.0015 | -0.0004 | 62.9397(1)
(4,6) | 59.0573 | -0.0058 | 0.0007 | 0.0065 | 0.0014 | -0.0004 | 59.0583(1)
(4,7) | 55.2755 | -0.0055 | 0.0007 | 0.0061 | 0.0013 | -0.0004 | 55.2765(1)
(4,8) | 51.5753 | -0.0051 | 0.0006 | 0.0057 | 0.0012 | -0.0003 | 51.5761(1)
(4,9) | 47.9381 | -0.0048 | 0.0006 | 0.0051 | 0.0010 | -0.0003 | 47.9388(1)
Table 12: Numerical results for the electron-deuteron spin-spin tensor
coefficient $E_{7}$ in HD+ for the range of rovibrational states $(L=1-4)$ and
$(v=0-9)$ (in kHz). All definitions are identical to those given in Table 2.


# On Age of Information for Discrete Time Status Updating System With
Ber/G/1/1 Queues

Jixiang Zhang and Yinfei Xu  School of Information Science and Engineering
Southeast University, Nanjing 210096, China  Email: {230179077,
<EMAIL_ADDRESS>
###### Abstract

In this paper, we consider the age of information (AoI) of a discrete time
status updating system, focusing on finding the stationary AoI distribution
assuming that the Ber/G/1/1 queue is used. Following the standard queueing
theory, we show that by invoking a two-dimensional state vector which tracks
the AoI and packet age in system simultaneously, the stationary AoI
distribution can be derived by analyzing the steady state of the constituted
two-dimensional stochastic process. We give the general formula of the AoI
distribution and calculate the explicit expression when the service time is
also geometrically distributed. The discrete and continuous AoI are compared,
we depict the mean of discrete AoI and that of continuous time AoI for system
with M/M/1/1 queue. Although the stationary AoI distribution of some
continuous time single-server system has been determined before, in this
paper, we shall prove that the standard queueing theory is still appliable to
analyze the discrete AoI, which is even stronger than the proposed methods
handling the continuous AoI.

##  I Introduction

Nowadays, real time data exchange and information transmission are getting
more and more important, especially in those IoT applications. A network node
needs fresh information to implement certain computations or make resource
scheduling. The age of information (AoI) metric was proposed in [1] to measure
the freshness of the messages obtained at the receiver, which is defined as
the time a packet experiences since it was generated at the source up to now.
Since then, lots of papers have been published considering different aspects
of AoI, both in theory and in practice [2].

In work [3], the author determined the limiting time average AoI for the
updating system with M/M/1, M/D/1 and D/M/1 queues, respectively. The mean of
the AoI for the basic queues using LCFS discipline were obtained in [4]. The
authors of paper [5] considered the system with multiple sources and derived
the average AoI of each source, using the newly proposed method which is
called the Stochastic Hybrid System (SHS) analysis. For the case the system
has two parallel transmitters, the closed-form expression of average AoI was
determined in work [6] by sophisticated random events analysis. Another
related metric called the peak age of information (PAoI) was introduced in
[7], which is defined as the peak value of the AoI just before it drops when
every time a packet is received by the reveiver. The AoI and peak AoI analysis
for three queue models was carried out in [7], including M/M/1/1, M/M/1/2 and
M/M/1/2*, where in the last queue model the packet waiting in queue can be
replaced constantly by the newly arriving packets, which is newer in time.

We observe that the analysis of the discrete time AoI was considered in paper
[8], where the authors determined the mean of the AoI and peak AoI for system
with Ber/G/1 and G/G/$\infty$ queues by the results obtained in continuous AoI
analysis. After then, the AoI of discrete time status updating system was also
analyzed in work [9]. The authors gave a new description for one sample path
of the AoI stochastic process. Provided these results, a general expression of
the generation function of the AoI and peak was derived, which is equivalent
to obtain their stationary distributions. Recently, using the queueing theory
methods, the discrete time AoI and peak AoI distributions are numerically
determined in series work [10, 11, 12]. In [13], the stationary distribution
of continuous AoI is determined for single-server system with various queue
models, and it seems that the ideas and methods used in [9] are from paper
[13].

In this paper, we consider the discrete AoI of status updating system with
Ber/G/1/1 queue, concentrate on determining the stationary distribution of the
AoI at the destination. In [9], the authors have obtained the AoI distribution
assuming that the queue used is Geo/Geo/1. Our work is different from theirs
in two aspects. Firstly, in their model, the queue buffer is unbounded while
we consider the bufferless system here. The AoI analysis for system having
infinite size is _different_ from that of finite size system. Secondly, it is
not clear whether the service time distribution can be relaxed to follow an
arbitrary distribution in [9], since the AoI distribution of the generalized
system is not given explicitly. While in this work, we show that by defining a
two-dimensional state vector which records the instantaneous AoI and the age
of the packet in system, all the random transitions can be described following
standard queueing theory. Therefore, the steady state of the established two-
dimensional process can give the stationary AoI, because the AoI is contained
as part of the state vector. It should be pointed out that our basic idea is
straightforward and different from the methods proposed and used in [9] and
[13].

The concept of multi-dimensional state vector is also introduced in[5], where
it makes up the SHS analysis for continuous time status updating system. In
fact, what we do here is to generalize SHS analysis to discrete time cases.
Due to the difficulty of solving the system of rate (or intensity)-balance
equations, only the mean or some simple moments of AoI can be calculated via
the SHS analysis in [5]. However, we will show that when the discrete AoI is
considered, the system of probability-balance equations can be solved
completely for some kinds of status updating systems, such that the stationary
distribution of any state component in defined muiti-dimensional state vector
can be obtained, not only the distribution of AoI.

The stationary AoI distribution of many single-server system has been
determined by its Laplace-Stieltjes transformation (LST) in [13].
Nevertheless, calculating the reverse transformation is not easy. The
stationary distribution usually has no closed-form expression. Due to this, it
is not easy to reduce the continuous results of AoI analysis to discrete
versions by direct discretization. On the contrary, as long as the system of
stationary equations are solved, we can get the explicit expression of
discrete AoI distribution directly. Moreover, we can approach the continuous
AoI distribution beginning with a approximate discrete time system, since any
continuous probability distribution can be approximated by a $k$-point
discrete one. Along this way, the AoI probability density function is
established, avoiding performing the complex reverse LST transformation.

The rest of this paper is organized as follows. We describe the basic status
updating system and give the necessary definition in Section II. The main
results are given in Section III. We obtain the general formula of the
stationary AoI distribution assuming that Ber/G/1/1 queue is used in system.
As a specific example, AoI distribution expression of system with Ber/Geo/1/1
queue is calculated. The distribution curves for different system parameters
are depicted in Section IV. In the concluding remark of Section V, we discuss
how to generalize our analysis method to more AoI systms.

##  II System Model and Problem Formulation

The basic status updating system consists of a source node and a destination
node. The source observes a physical process $X(t)$ and samples the states of
$X(t)$ at random times. An updating packet is composed of the sampled state
$X(t_{i})$ of the process and the sample time $t_{i}$. The source sends the
packets to destination via a transmitter, which is modeled as a queue. In
Figure 1, the model of a status updating system is depicted.

Figure 1: The model of a status updating system.

In the discrete time model, both the time and the AoI are discretized. During
the time when no packet is obtained at the receiver, the AoI increases
constantly. When the receiver gets an update, the value of the AoI reduces to
the system time of the arriving packet at the end of that time slot. In
general, the system time of a packet equals the sum of the waiting time in
queue and the service time the packet consumes at the server. However, for the
queue model used in this paper, the packet waiting time equals zero.

Now, we give the definition of the time average AoI under the discrete time
setting.

######  Definition 1.

In the discrete time model, the limiting time average age of information
$\overline{\Delta}$ is defined as

$\overline{\Delta}=\lim_{T\to\infty}\frac{1}{T}\sum\nolimits_{k=1}^{T}a(k)$
(1)

where $a(k)$ denotes the value of AoI at $k$th time slot.

We simply deal with the expression (1) as follows.

$\displaystyle\overline{\Delta}$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\sum\nolimits_{k=1}^{T}a(k)$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\sum\nolimits_{n=1}^{M}n\cdot|\\{k:a(k)=n\\}|=\sum\nolimits_{n=1}^{\infty}n\cdot\pi_{n}$

where

$\pi_{n}=\lim_{T\to\infty}\frac{|\\{k:a(k)=n\\}|}{T}\qquad(n\geq 1)$

is the probability that the AoI takes value $n$ when the observing time tends
to infinity. The number $M$ is defined to be $M=\text{max}_{1\leq k\leq
T}a(k)$.

Always the AoI processes are assumed to be ergodic. Ergodicity property
ensures that we can obtain the AoI performance metrics by investigating any
one sample path of the process. Notice that if all the probabilities
$\pi_{n}$, $n\geq 1$ are determined, then the distribution of the AoI is known
as well.

##  III Discrete AoI-distribution: Ber/G/1/1 queues

In this Section, assume that the queue model is Ber/G/1/1, we derive the
general expression of the stationary AoI distribution.

###  III-A Stationary distribution of the AoI: Ber/G/1/1 queues

Assume that at each time slot, a new updating packet comes independently and
with identical probability $p$. We use $A$ to denote the interarrival time
between successive arriving packets, then the distribution of $A$ is

$\Pr\\{A=j\\}=(1-p)^{j-1}p\qquad(j\geq 1)$

The packet service time is represented by $B$, suppose that its distribution
is written as

$\Pr\\{B=j\\}=q_{j}\qquad(j\geq 1)$

The updating packet comes at the beginning of one time slot, while the AoI is
recorded and modified at the end of each time slot. The terms state, age-state
and the state vector are used interchangeably.

Define the two-dimensional vector $(n_{k},m_{k})$, where $n_{k}$ denotes the
value of the AoI at the $k$th time slot and the second component $m_{k}$
represents the age of the packet in system at that time. Constituting the two-
dimensional stochastic process

$AoI_{Ber/G/1/1}=\left\\{(n_{k},m_{k}):n_{k}>m_{k}\geq 0,k\geq 1\right\\}$

In the following paragraphs, we analyze all the random transfers of age-states
$(n,m)$ and determine the transtion probabilities, such that the stationary
equations describing the steady state of the AoI process can be established.
Notice that the first component of the two-dimensional state vector, $n$,
denotes the AoI at the destination. Therefore, as long as all the stationary
probabilities of state $(n,m)$ are determined, the distribution of the AoI is
obtained either. An illustration of the discrete state and state transfers are
provided in Figure 2.

Figure 2: An illustration of discrete state and state transfers: Status
updating system with Ber/G/1/1 queue.

First of all, assume that the state vector is $(n,m)$ at the $k$th time slot,
where $n>m>0$. At the $(k+1)$th time slot, it shows that $(n,m)$ can only
transfer to either $(n+1,m+1)$ or $(m+1,0)$ depending on whether or not the
packet service is complete. Because $m$ is greater than zero, which means an
updating packet is served at the transmitter currently and the system is full.
Let $P_{(n,m),(i,j)}$ be the transition probability the age-state jumps to
$(i,j)$ at the next time slot from the current state $(n,m)$. We have that

$P_{(n,m),(n+1,m+1)}=\Pr\\{B>m|B>m-1\\}$ (2)

The conditional probability is used in equation (2) because when the service
time $B$ is longer than $m$ time slots, we must know the event $\\{B>m-1\\}$.

On the contrary, if the service of the packet finishes at the next time slot,
the system is emptied and the age-state changes to $(m+1,0)$. Therefore, we
have

$P_{(n,m),(m+1,0)}=\Pr\\{B=m|B>m-1\\}$ (3)

Next, consider the case $m=0$. At this time, the transmitter is idle and the
system is empty. Thus, the state transitions are totally governed by the
interarrival time $A$.

For the state vector $(n,0)$, $n\geq 1$, if no packet comes at next time slot,
then

$P_{(n,0),(n+1,0)}=1-p$ (4)

Assume that a new packet comes at next time slot. In this case, we have to
consider two sub-cases further. When the packet comes and does not leave after
one time slot, the age-state changes to $(n+1,1)$. Otherwise, if the arriving
packet finishes the service in single time slot and gets to receiver, we show
that the state vector will jump to $(1,0)$. Summarize above discussions, we
obtain

$P_{(n,0),(n+1,1)}=p\Pr\\{B>1\\},\quad P_{(n,0),(1,0)}=p\Pr\\{B=1\\}$

So far, all the state transfers are analyzed and the corresponding transition
probabilities are determined. Define $\pi_{(n,m)}$ as the stationary
probability of the state $(n,m)$. Then, we give the stationary equations of
the AoI process $AoI_{Ber/G/1/1}$.

######  Theorem 1.

Denote $\pi_{(n,m)}$ as the stationary probability of the age-state $(n,m)$,
where $n>m\geq 0$. Then, for the two-dimensional stochastic process
$AoI_{Ber/G/1/1}$, the stationary equations are determined as

$\begin{cases}\pi_{(n,m)}=\pi_{(n-1,m-1)}\Pr\\{B>m-1|B>m-2\\}\\\
\qquad\qquad\qquad\qquad\qquad\qquad\quad\quad\quad(n>m\geq 2)\\\
\pi_{(n,1)}=\pi_{(n-1,0)}p\Pr\\{B>1\\}\qquad\quad\quad(n\geq 2)\\\
\pi_{(n,0)}=\pi_{(n-1,0)}(1-p)+\left(\sum\nolimits_{k=n}^{\infty}\pi_{(k,n-1)}\right)\\\
\qquad\quad\times\Pr\\{B=n-1|B>n-2\\}\qquad(n\geq 2)\\\
\pi_{(1,0)}=\left(\sum\nolimits_{k=1}^{\infty}\pi_{(k,0)}\right)p\Pr\\{B=1\\}\end{cases}$
(5)

###### Proof.

We explain each line of (5) as follows. First of all, for the case $n>m\geq
2$, the system is full and is still full before one time slot. The newly
arriving packet cannot enter the transmitter if there exists such one, so that
the state transitions are totally determined by the random service time $B$.
Assume that the current state is $(n-1,m-1)$, then with probability
$\Pr\\{B>m-1|B>m-2\\}$ the age-state jumps to $(n,m)$ at next time slot. This
gives the first line of (5).

Consider the state $(n,1)$ in second line, notice that before one time slot
the system is empty because $m$ reduces to 0. Starting with $(n-1,0)$, let a
packet comes and stays at the system, which occurs with probability
$p\Pr\\{B>1\\}$, we show that the age-state transfers to $(n,1)$ at next time
slot.

Now, the state transitions of age-state $(n,0)$, $n\geq 2$ are analyzed. From
the state $(n-1,0)$, as long as no packet arrives in next time slot, the value
of the AoI at the receiver increases 1 and the second component $m$ remains
zero, we get the state vector $(n,0)$. On the other hand, assume that the
original state is $(k,n-1)$ and the service of the packet finishes at the next
time slot, it shows that the state will transfer to $(n,0)$ as well. Combining
all the transitions from $(k,n-1)$ to $(n,0)$ where $k\geq n$, we obtain the
stationary equations for state vectors $(n,0)$, $n\geq 2$ and explain the
third line of (5).

Finally, beginning with an empty system, an arriving packet experiences single
time slot service and then delivers to destination will make the age-state
jump to $(1,0)$, which yields the last equation in (5). ∎

As long as all the probabilities $\pi_{(n,m)}$ are determined by solving the
system of equations (5), the stationary distribution of the AoI can be
obtained since we have

$\Pr\\{\Delta=n\\}=\sum\nolimits_{m=0}^{n-1}\pi_{(n,m)}\qquad(n\geq 1)$ (6)

Furthermore, for $k\geq 1$, the distribution function of the AoI can also be
obtained as

$\Pr\\{\Delta\leq
k\\}=\sum\nolimits_{n=1}^{k}\Pr\\{\Delta=n\\}=\sum\nolimits_{n=1}^{k}\sum\nolimits_{m=0}^{n-1}\pi_{(n,m)}$

We solve the system of equations (5) and determine all the stationary
probabilities $\pi_{(n,m)}$ in Theorem 2.

######  Theorem 2.

When the AoI process $AoI_{Ber/G/1/1}$ reaches the steady state, the
stationary probabilities for all the age-states $(n,m)$ are given as follows.
Firstly,

$\pi_{(n,0)}=\frac{pq_{1}F(p,n)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}\quad(n\geq
1)$ (7)

and the probabilities $\pi_{(n,m)}$, $n>m\geq 1$ are solved as

$\pi_{(n,m)}=\frac{p^{2}q_{1}(1-q_{1})F(p,n-m)\left(\sum\nolimits_{l=m}^{\infty}q_{l}\right)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}$
(8)

where for $n\geq 2$, define that

$\displaystyle
F(p,n)=(1-p)^{n-1}+\frac{1-q_{1}}{q_{1}}\sum\nolimits_{j=0}^{n-2}(1-p)^{j}q_{n-1-j}$
(9)

and denote $F(p,1)=1$.

The proof of Theorem 2 is given in Appendix A.

Given all the stationary probabilities, the distribution of AoI can be
determined by equation (6).

######  Corollary 1.

Assume that the queue used in status updating system is Ber/G/1/1. When the
system runs in steady state, the stationary AoI-distribution is calculated as

$\displaystyle\Pr\\{\Delta=n\\}=\pi_{(n,0)}+\sum\nolimits_{m=1}^{n-1}\pi_{(n,m)}$
$\displaystyle={}$
$\displaystyle\frac{pq_{1}F(p,n)+p^{2}q_{1}(1-q_{1})\sum\nolimits_{m=1}^{n-1}F(p,n-m)\left(\sum\nolimits_{l=m}^{\infty}q_{l}\right)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}$

The AoI cumulative probabilities are determined as

$\Pr\\{\Delta\leq k\\}=\sum\nolimits_{n=1}^{k}\Pr\\{\Delta=n\\}\quad(k\geq 1)$

where $F(p,n)$ is defined in equation (9).

###### Proof.

Since the first component of the state vector denotes the AoI at the receiver,
the probability that the AoI equals $n$ can be obtained by adding up all the
stationary probabilities having $n$ as first component. We show that

$\displaystyle\Pr\\{\Delta=n\\}=\sum\nolimits_{m=0}^{n-1}\pi_{(n,m)}$
$\displaystyle={}$
$\displaystyle\pi_{(n,0)}+\sum\nolimits_{m=1}^{n-1}\pi_{(n,m)}$
$\displaystyle={}$
$\displaystyle\frac{pq_{1}F(p,n)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}$
$\displaystyle+\sum\nolimits_{m=1}^{n-1}\frac{p^{2}q_{1}(1-q_{1})F(p,n-m)\left(\sum\nolimits_{l=m}^{\infty}q_{l}\right)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}$
$\displaystyle={}$
$\displaystyle\frac{pq_{1}F(p,n)+p^{2}q_{1}(1-q_{1})\sum\nolimits_{m=1}^{n-1}F(p,n-m)\left(\sum\nolimits_{l=m}^{\infty}q_{l}\right)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}$

For $k\geq 1$, the cumulative probability of the AoI equals the sum of
probabilities the AoI takes value $n$ from $n=1$ to $k$. This completes the
proof of Corollary 1. ∎

###  III-B The AoI stationary distribution: Ber/Geo/1/1 queue

When the packet service time $B$ is also a geometric random variable, we can
calculate the explicit expression of AoI distribution. Suppose that $B$ is
distributed as

$\Pr\\{B=k\\}=(1-\gamma)^{k-1}\gamma\quad(k\geq 1)$ (10)

To keep queueing system stable, let $p<\gamma$. We first find all the
stationary probabilities by solving the system of equations (5), and then
derive the explicit expression of AoI-distribution. The results are stated in
Theorem 3 below.

######  Theorem 3.

When the service time $B$ is also geometrically distributed, the stationary
probabilities $\pi_{(n,m)}$ are solved as

$\begin{cases}\pi_{(n,0)}=\frac{p\gamma^{2}\left[(1-p)^{n}-(1-\gamma)^{n}\right]}{(p+\gamma-p\gamma)(\gamma-p)}&(n\geq
1)\\\
\pi_{(n,m)}=\frac{(p\gamma)^{2}\left[(1-p)^{n-m}(1-\gamma)^{m}-(1-\gamma)^{n}\right]}{(p+\gamma-p\gamma)(\gamma-p)}&(n>m\geq
1)\end{cases}$

The Theorem 3 can be proved either from the general formulas (7)-(9) or
solving system of equation (5) directly. We give the derivation details in
Appendix B.

By collecting all the stationary probabilities with identical AoI-component,
we can derive the stationary distribution of the AoI. The cumulative
probability distribution of the AoI is calculated accordingly. The calculation
detalis are omitted.

######  Corollary 2.

The distribution of AoI for the status updating system with Ber/Geo/1/1 queue
is given as

$\displaystyle\Pr\\{\Delta=n\\}$
$\displaystyle=\frac{p(1-p)\gamma^{3}\left[(1-p)^{n}-(1-\gamma)^{n}\right]}{(p+\gamma-p\gamma)(\gamma-p)^{2}}$
$\displaystyle\qquad-\frac{(p\gamma)^{2}n(1-\gamma)^{n}}{(p+\gamma-p\gamma)(\gamma-p)}\qquad(n\geq
1)$ (11)

and the AoI cumulative probability is determined as

$\displaystyle\Pr\\{\Delta\leq k\\}$
$\displaystyle=1-\frac{(1-p)\gamma^{2}\left[(1-p)^{k+1}\gamma-p(1-\gamma)^{k+1}\right]}{(p+\gamma-p\gamma)(\gamma-p)^{2}}$
$\displaystyle\quad+\frac{p^{2}(1+k\gamma)(1-\gamma)^{k+1}}{(p+\gamma-p\gamma)(\gamma-p)}\qquad(k\geq
1)$ (12)

(a) Stationary distribution of discrete AoI

(b) Discrete AoI cumulative probabilities

(c) Average continuous and discrete AoI

Figure 3: Numerical simulations of discrete AoI and discrete approximations of
mean of continuous AoI.

Notice that in order to obtain the stationary distribution of AoI, we invoke a
two-dimensional state vector $(n,m)$ where the first component denotes the
instantaneous AoI at the destination, while the latter parameter represents
the age of the packet in system. Since the system size is 1, so that the
packet waiting time is zero, the second component $m$ is exactly the packet
service time, which should be geometrically distributed in this case. We prove
this claim by computing the marginal distribution of $M$, which is defined as
the random variable of the second component of the age-state.

For $m\geq 1$, we show that

$\displaystyle\Pr\\{M=m\\}=\sum\nolimits_{n=m+1}^{\infty}\pi_{(n,m)}$
$\displaystyle={}$
$\displaystyle\sum\nolimits_{n=m+1}^{\infty}\frac{(p\gamma)^{2}\left[(1-p)^{n-m}(1-\gamma)^{m}-(1-\gamma)^{n}\right]}{(p+\gamma-p\gamma)(\gamma-p)}$
$\displaystyle={}$
$\displaystyle\frac{(p\gamma)^{2}}{(p+\gamma-p\gamma)(\gamma-p)}\left[\frac{1-p}{p}(1-\gamma)^{m}-\frac{1-\gamma}{\gamma}(1-\gamma)^{m}\right]$
$\displaystyle={}$
$\displaystyle\frac{p\gamma}{p+\gamma-p\gamma}(1-\gamma)^{m}$ (13)

We have used the probability expressions obtained in Theorem 3. Probability
(13) should be normalized by the probability that the transmitter is not idle,
since in that case we have $m=0$. It shows that

$\displaystyle\Pr\\{\text{system is not idle}\\}$ $\displaystyle={}$
$\displaystyle
1-\sum\nolimits_{n=1}^{\infty}\pi_{(n,0)}=\frac{p(1-\gamma)}{p+\gamma-p\gamma}$

Therefore, for $m\geq 1$, we obtain that

$\Pr\\{B=m\\}=\frac{\Pr\\{M=m\\}}{\Pr\\{\text{system is not
idle}\\}}=(1-\gamma)^{m-1}\gamma$

which is indeed the geometric distribution with parameter $\gamma$.

From (11), the mean of the discrete AoI can be calculated. We will show that
when the time slot is sufficiently small, the average discrete AoI is close to
average continuous AoI, which is denoted by $\overline{\Delta}_{M/M/1/1}$ and
is determined in [7] as

$\overline{\Delta}_{M/M/1/1}=\frac{1}{\mu}\left(1+\frac{1}{\rho}+\frac{\rho}{1+\rho}\right)$

######  Corollary 3.

For the system with Ber/Geo/1/1 queue, the limiting time average discrete AoI
is equal to

$\overline{\Delta}^{(d)}_{Ber/Geo/1/1}=\frac{1}{\gamma}\left((1-\gamma)+\frac{1}{\rho_{d}}+\frac{\rho_{d}}{\frac{1}{1-\gamma}+\rho_{d}}\right)$
(14)

and approaches $\overline{\Delta}_{M/M/1/1}$ in almost all the range of
$\rho_{d}$, which is discrete traffic intensity and is defined as
$\rho_{d}=p/\gamma$.

The Proof of Corollary 3 is shown in Appendix C.

##  IV Numerical Results

In Figure 3, we draw the stationary distribution of discrete AoI along with
its cumulative probabilities. In addition, the mean of discrete and continuous
AoI are also depicted.

We take two group of parameters $(p,\gamma)$ with fixed ratio 0.5. From Figure
3a, it sees that the AoI distribution curve is more concentrated when $p$ and
$\gamma$ are both larger, and has higher peak probability. The cumulative
probability distributions also show that when the parameters are larger, the
AoI cumulative probability gets to 1 more faster.

Average AoIs of system with M/M/1/1 and Ber/Geo/1/1 queue are given in Figure
3c. We draw both curves for three different cases, where the service rate
$\mu$ (or $\gamma$) is fixed and the traffic intensity varies. Numerical
results show that in almost all the ranges of $\rho$ (or $\rho_{d}$), the mean
of discrete AoI is very close to its continuous counterpart.

##  V Conclusion

In this paper, we find the stationary distribution of discrete AoI for the
system with Ber/G/1/1 queue. By defining a two-dimensional age-state which
contains AoI and describing all the random state transitions, we show that the
AoI distribution can be obtained if the steady state of constituted two-
dimensional stochastic process is solved. The idea of using a multiple-
dimensional state vector comes from SHS analysis of AoI. What we do is
generalizing SHS method to discrete time status updating systems. The
numerical results show that the mean of discrete AoI is very close to the mean
of continuous AoI when time slot is small enough.

To obtain the exact probability density function of continuous AoI, often a
reverse LST transformation has to be computed which is not easy. We point out
that although the proposed idea is straightforward, it is still appliable to
analyze AoI of more systems. The biggest advantage is that the exact
expression of AoI distribution can be obtained as long as the steady state of
constituted AoI stochastic process is solved.

## References

  * [1] S. Kaul, M. Gruteser, V. Rai, and J. Kenney, “Minimizing age of information in vehicular networks,” _8th Annual IEEE Communications Society Conference on Sensor, Mesh and Ad Hoc Communications and Networks_ , pp. 350–358, 2011.
  * [2] R. D. Yates, Y. Sun, D. R. B. III, S. K. Kaul, E. Modiano, and S. Ulukus, “Age of information: An introduction and survey,” _arXiv:2007.08564v1_ , 2020\.
  * [3] S. Kaul, R. Yates, and M. Gruteser, “Real-time status: How often should one update?” _2012 Proceedings IEEE INFOCOM_ , pp. 2731–2735, 2012.
  * [4] ——, “Status updates through queues,” _46th Annual Conference on Information Sciences and Systems (CISS)_ , 2012.
  * [5] R. D. Yates and S. Kaul, “The age of information: Real-time status updating by multiple sources,” _IEEE Transactions on Information Theory_ , vol. 65, no. 3, pp. 1807–1827, 2019.
  * [6] C. Kam, S. Kompella, G. D. Nguyen, and A. Ephremides, “Effect of message transmission path diversity on status age,” _IEEE Transactions on Information Theory_ , vol. 62, no. 3, pp. 1360–1374, 2016.
  * [7] M. Costa, M. Codreanu, and A. Ephremides, “Age of information with packet management,” _2014 IEEE International Symposium on Information Theory_ , pp. 1583–1587, 2014.
  * [8] V. Tripathi, R. Talak, and E. Modiano, “Age of information for discrete time queues,” _arXiv:1901.10463v1_ , 2019.
  * [9] A. Kosta, N. Pappas, A. Ephremides, and V. Angelakis, “Non-linear age of information in a discrete time queue: Stationary distribution and average performance analysis,” _ICC 2020 - 2020 IEEE International Conference on Communications (ICC)_ , pp. 1–6, 2020.
  * [10] N. Akar, O. Doğan, and E. U. Atay, “Finding the exact distribution of (peak) age of information for queues of ph/ph/1/1 and m/ph/1/2 type,” _IEEE Transactions on Communications_ , vol. 68, no. 9, pp. 5661–5672, 2020.
  * [11] N. Akar and O. Dogan, “Discrete-time queueing model of age of information with multiple information sources,” _arXiv:2007.11650v1_ , 2020.
  * [12] O. Dogan and N. Akar, “The multi-source preemptive m/ph/1/1 queue with packet errors: Exact distribution of the age of information and its peak,” _arXiv:2007.11656v1_ , 2020.
  * [13] Y. Inoue, H. Masuyama, T. Takine, and T. Tanaka, “A general formula for the stationary distribution of the age of information and its application to single-server queues,” _IEEE Transactions on Information Theory_ , vol. 65, no. 12, pp. 8305–8324, 2019.

##  Appendix A Proof of Theorem 2

In this Appendix, we solve the system of equations (5) and derive all the
stationary probabilities $\pi_{(n,m)}$. For convenience, the stationary
equations are copied in the following.

$\begin{cases}\pi_{(n,m)}=\pi_{(n-1,m-1)}\Pr\\{B>m-1|B>m-2\\}&(n>m\geq 2)\\\
\pi_{(n,1)}=\pi_{(n-1,0)}p\Pr\\{B>1\\}&(n\geq 2)\\\
\pi_{(n,0)}=\pi_{(n-1,0)}(1-p)+\left(\sum\nolimits_{k=n}^{\infty}\pi_{(k,n-1)}\right)\Pr\\{B=n-1|B>n-2\\}&(n\geq
2)\\\
\pi_{(1,0)}=\left(\sum\nolimits_{k=1}^{\infty}\pi_{(k,0)}\right)p\Pr\\{B=1\\}\end{cases}$
(15)

Firstly, we show that the equations (17) can be simplified by eliminating the
dependence on the second component $m$. For a general age-state $(n,m)$,
$n>m\geq 2$, iteratively applying the first line of (17) yields

$\displaystyle\pi_{(n,m)}$ $\displaystyle=\pi_{(n-1,m-1)}\Pr\\{B>m-1|B>m-2\\}$
$\displaystyle=\pi_{(n-2,m-2)}\Pr\\{B>m-2|B>m-3\\}\Pr\\{B>m-1|B>m-2\\}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\quad\vdots$
$\displaystyle=\pi_{(n-m+1,1)}\Pr\\{B>1\\}\times\cdots\times\Pr\\{B>m-1|B>m-2\\}$
$\displaystyle=\pi_{(n-m+1,1)}\Pr\\{B>m-1\\}$
$\displaystyle=\pi_{(n-m,0)}p(1-q_{1})\Pr\\{B>m-1\\}$ (16)

For the last step recursion, we use the second line of (17).

Equation (18) shows that probability $\pi_{(n,m)}$, $n>m\geq 1$ can be
represented by $\pi_{(n,0)}$, $n\geq 1$. As a result, to solve system of
equations (17) we only need to determine those stationary probabilities
$\pi_{(n,0)}$.

Next, the infinite sum in the third line of (17) is calculated. We show that

$\displaystyle\left(\sum\nolimits_{k=n}^{\infty}\pi_{(k,n-1)}\right)\Pr\\{B=n-1|B>n-2\\}$
$\displaystyle={}$
$\displaystyle\left(\sum\nolimits_{k=n}^{\infty}\pi_{(k-n+1,0)}p(1-q_{1})\Pr\\{B>n-2\\}\right)\Pr\\{B=n-1|B>n-2\\}$
$\displaystyle={}$
$\displaystyle\left(\sum\nolimits_{k=n}^{\infty}\pi_{(k-n+1,0)}\right)p(1-q_{1})\Pr\\{B=n-1\\}$
$\displaystyle={}$ $\displaystyle
p(1-q_{1})q_{n-1}\left(\sum\nolimits_{k=1}^{\infty}\pi_{(k,0)}\right)$
$\displaystyle={}$ $\displaystyle\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}q_{n-1}$
(17)

In order to obtain equation (19), notice that the last equation in (17)
implies that

$\sum\nolimits_{k=1}^{\infty}\pi_{(k,0)}=\frac{\pi_{(1,0)}}{pq_{1}}$ (18)

Substituting (19) into the third line of (17), we derive the following
recursive formula for the probabilities $\pi_{(n,0)}$.

$\pi_{(n,0)}=\pi_{(n-1,0)}(1-p)+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}q_{n-1}\qquad(n\geq
2)$ (19)

Repeatedly using (21) yields following equations

$\displaystyle\pi_{(n,0)}$
$\displaystyle=\pi_{(n-1,0)}(1-p)+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}q_{n-1}$
$\displaystyle=\left(\pi_{(n-2,0)}(1-p)+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}q_{n-2}\right)(1-p)+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}q_{n-1}$
$\displaystyle=\pi_{(n-2,0)}(1-p)^{2}+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}\left[(1-p)q_{n-2}+q_{n-1}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\vdots$
$\displaystyle=\pi_{(1,0)}(1-p)^{n-1}+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}\left(\sum\nolimits_{j=0}^{n-2}(1-p)^{j}q_{n-1-j}\right)$
(20)

Let $n=1$ in (22), we obtain the obvious equation $\pi_{(1,0)}=\pi_{(1,0)}$.
Thus, the expression (22) is valid for all $n\geq 1$.

Since all the probabilities $\pi_{(n,m)}$ add up to 1, we have

$\displaystyle 1$
$\displaystyle=\sum\nolimits_{m=0}^{\infty}\sum\nolimits_{n=m+1}^{\infty}\pi_{(n,m)}$
$\displaystyle=\sum\nolimits_{n=1}^{\infty}\pi_{(n,0)}+\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{n=m+1}^{\infty}\pi_{(n,m)}$
$\displaystyle=\sum\nolimits_{n=1}^{\infty}\pi_{(n,0)}+\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{n=m+1}^{\infty}\pi_{(n-m,0)}p(1-q_{1})\Pr\\{B>m-1\\}$
$\displaystyle=\sum\nolimits_{n=1}^{\infty}\pi_{(n,0)}+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\left(\sum\nolimits_{n=1}^{\infty}\pi_{(n,0)}\right)\Pr\\{B>m-1\\}$
$\displaystyle=\left(\sum\nolimits_{n=1}^{\infty}\pi_{(n,0)}\right)\left[1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\Pr\\{B>m-1\\}\right]$
$\displaystyle=\frac{\pi_{(1,0)}}{pq_{1}}\left[1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\Pr\\{B>m-1\\}\right]$
(21)

from which the first probability $\pi_{(1,0)}$ can be determined as

$\pi_{(1,0)}=\frac{pq_{1}}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\Pr\\{B>m-1\\}}$
(22)

Define

$F(p,n)=\begin{cases}1&\qquad(n=1)\\\
(1-p)^{n-1}+\frac{1-q_{1}}{q_{1}}\sum\nolimits_{j=0}^{n-2}(1-p)^{j}q_{n-1-j}&\qquad(n\geq
2)\end{cases}$

Thus, according to equation (22), we can write that

$\displaystyle\pi_{(n,0)}=\pi_{(1,0)}F(p,n)=\frac{pq_{1}F(p,n)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}\qquad(n\geq
1)$ (23)

The other probabilities $\pi_{(n,m)}$ can be determined by using recursive
expression (18) as

$\displaystyle\pi_{(n,m)}$
$\displaystyle=\pi_{(n-m,0)}p(1-q_{1})\Pr\\{B>m-1\\}$
$\displaystyle=\pi_{(1,0)}F(p,n-m)p(1-q_{1})\left(\sum\nolimits_{l=m}^{\infty}q_{l}\right)$
$\displaystyle=\frac{p^{2}q_{1}(1-q_{1})F(p,n-m)\left(\sum\nolimits_{l=m}^{\infty}q_{l}\right)}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{l=m}^{\infty}q_{l}}$
(24)

So far, we completely solve the system of equations (17) and obtain all the
stationary probabilities. This completes the proof of Theorem 2.

##  Appendix B Proof of Theorem 3

In this appendix, we solve the system of stationary equations when the queue
adopted in system is Ber/Geo/1/1 queue.

It is not difficult to obtain the solutions in Theorem 3 by using the general
results (6)-(9). Firstly, from the equation (22) we have

$\displaystyle\pi_{(n,0)}$
$\displaystyle=\pi_{(1,0)}(1-p)^{n-1}+\frac{(1-q_{1})\pi_{(1,0)}}{q_{1}}\left(\sum\nolimits_{j=0}^{n-2}(1-p)^{j}q_{n-1-j}\right)$
$\displaystyle=\pi_{(1,0)}(1-p)^{n-1}+\frac{(1-\gamma)\pi_{(1,0)}}{\gamma}\left(\sum\nolimits_{j=0}^{n-2}(1-p)^{j}(1-\gamma)^{n-2-j}\gamma\right)$
$\displaystyle=\pi_{(1,0)}\bigg{[}(1-p)^{n-1}+(1-\gamma)^{n-1}\sum\nolimits_{j=0}^{n-2}\left(\frac{1-p}{1-\gamma}\right)^{j}\bigg{]}$
$\displaystyle=\pi_{(1,0)}\bigg{[}(1-p)^{n-1}+\frac{(1-p)^{n-1}(1-\gamma)-(1-\gamma)^{n}}{\gamma-p}\bigg{]}$
$\displaystyle=\pi_{(1,0)}\frac{(1-p)^{n}-(1-\gamma)^{n}}{\gamma-p}$ (25)

The first probability $\pi_{(1,0)}$ can be determined from equation (24) as

$\displaystyle\pi_{(1,0)}$
$\displaystyle=\frac{pq_{1}}{1+p(1-q_{1})\sum\nolimits_{m=1}^{\infty}\Pr\\{B>m-1\\}}$
$\displaystyle=\frac{p\gamma}{1+p(1-\gamma)\sum\nolimits_{m=1}^{\infty}\sum\nolimits_{k=m}^{\infty}(1-\gamma)^{k-1}\gamma}$
$\displaystyle=\frac{p\gamma}{1+p(1-\gamma)\sum\nolimits_{m=1}^{\infty}(1-\gamma)^{m-1}}$
$\displaystyle=\frac{p\gamma}{1+p(1-\gamma)(1/\gamma)}$
$\displaystyle=\frac{p\gamma^{2}}{p+\gamma-p\gamma}$ (26)

Combining (27) and (28), we obtain

$\pi_{(n,0)}=\frac{p\gamma^{2}\left[(1-p)^{n}-(1-\gamma)^{n}\right]}{(p+\gamma-p\gamma)(\gamma-p)}\qquad(n\geq
2)$ (27)

Notice that expression (29) is also valid for the case $n=1$.

The other probabilities $\pi_{(n,m)}$ can be found by calculating equation
(18) directly. For $n>m\geq 1$, we have

$\displaystyle\pi_{(n,m)}$
$\displaystyle=\pi_{(n-m,0)}p(1-q_{1})\Pr\\{B>m-1\\}$
$\displaystyle=\pi_{(n-m,0)}p(1-\gamma)\sum\nolimits_{k=m}^{\infty}(1-\gamma)^{k-1}\gamma$
$\displaystyle=\pi_{(n-m,0)}p(1-\gamma)^{m}$
$\displaystyle=\frac{(p\gamma)^{2}\left[(1-p)^{n-m}(1-\gamma)^{m}-(1-\gamma)^{n}\right]}{(p+\gamma-p\gamma)(\gamma-p)}$
(28)

So far, all the probabilities are determined. This completes the proof of
Theorem 3.

##  Appendix C Proof of Corollary 3

The average discrete AoI $\overline{\Delta}^{(d)}_{Ber/Geo/1/1}$ is calculated
as

$\displaystyle\overline{\Delta}^{(d)}_{Ber/Geo/1/1}$
$\displaystyle=\sum\nolimits_{n=1}^{\infty}n\Pr\\{\Delta=n\\}$
$\displaystyle=\frac{p(1-p)\gamma^{3}}{(p+\gamma-p\gamma)(\gamma-p)^{2}}\sum\nolimits_{n=1}^{\infty}\left[n(1-p)^{n}-n(1-\gamma)^{n}\right]-\frac{(p\gamma)^{2}}{(p+\gamma-p\gamma)(\gamma-p)}\sum\nolimits_{n=1}^{\infty}n^{2}(1-\gamma)^{n}$
$\displaystyle=\frac{p(1-p)\gamma^{3}}{(p+\gamma-p\gamma)(\gamma-p)^{2}}\frac{(p+\gamma-p\gamma)(\gamma-p)}{p^{2}\gamma^{2}}-\frac{(p\gamma)^{2}}{(p+\gamma-p\gamma)(\gamma-p)}\frac{p^{2}(1-\gamma)(2-\gamma)}{\gamma^{3}}$
$\displaystyle=\frac{(1-p)\gamma}{p(\gamma-p)}-\frac{p^{2}(1-\gamma)(2-\gamma)}{(p+\gamma-p\gamma)(\gamma-p)\gamma}$
$\displaystyle=\frac{(1-\rho_{d}\gamma)\gamma}{\rho_{d}\gamma(\gamma-\rho_{d}\gamma)}-\frac{\rho_{d}^{2}\gamma^{2}(1-\gamma)(2-\gamma)}{(\rho_{d}\gamma+\gamma-\rho_{d}\gamma^{2})(\gamma-\rho_{d}\gamma)\gamma}$
(29)
$\displaystyle=\frac{1}{\gamma}\left(\frac{1}{\rho_{d}}+\frac{1-\gamma}{1-\rho_{d}}-\frac{\rho_{d}(1-\gamma)}{1-\rho_{d}}+\frac{\rho_{d}(1-\gamma)}{1+(1-\gamma)\rho_{d}}\right)$
$\displaystyle=\frac{1}{\gamma}\left((1-\gamma)+\frac{1}{\rho_{d}}+\frac{\rho_{d}}{\frac{1}{1-\gamma}+\rho_{d}}\right)$
(30)

where in (29) we substitute $\rho_{d}=p/\gamma$.

It is easy to see that

$\overline{\Delta}^{(d)}_{Ber/Geo/1/1}\to\overline{\Delta}_{M/M/1/1}=\frac{1}{\mu}\left(1+\frac{1}{\rho}+\frac{\rho}{1+\rho}\right)$

when $(1-\gamma)\to 1$, i.e., when the time slot is short enough.

This completes the proof of Corollary 3.