# The assignments of the $B_{s}$ mesons within the screened potential model
and ${}^{3}P_{0}$ model

Wei Hao<EMAIL_ADDRESS>School of Physics and Microelectronics,
Zhengzhou University, Zhengzhou, Henan 450001, China  CAS Key Laboratory of
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190,China  University of Chinese Academy of Sciences
(UCAS), Beijing 100049, China  Yu Lu<EMAIL_ADDRESS>School of Physical
Sciences, University of Chinese Academy of Sciences (UCAS), Beijing 100049,
China  En Wang<EMAIL_ADDRESS>School of Physics and Microelectronics,
Zhengzhou University, Zhengzhou, Henan 450001, China

###### Abstract

We investigate the mass spectrum and the decay properties of the $B_{s}$
mesons within the screened nonrelativistic quark model and the ${}^{3}P_{0}$
model. Our results suggest that the $B_{sJ}(6064)$ and $B_{sJ}(6114)$ states,
as the first solution of the recently LHCb measurements, could be explained as
the $B_{s}(1^{3}D_{3})$ and $B_{s}(1^{3}D_{1})$, respectively. In addition,
the $B_{sJ}(6109)$ and $B_{sJ}(6158)$ states, as the second solution of the
LHCb measurements, could be explained as the $B^{\prime}_{s2}(1D)$ and
$B_{s1}(2P)$, respectively. Meanwhile, the $B_{s1}(5830)$ could be interpreted
as the candidate of the $B_{s1}(1P)$. We also calculated the decay properties
of the other excited $B_{s}$ mesons with the predicted masses, which should be
helpful for the experimental searching in future.

##  I Introduction

Recently, the LHCb Collaboration observed an excess structure 300 MeV above
the $B^{\pm}K^{\mp}$ threshold in the $B^{\pm}K^{\mp}$ mass spectrum in the
proton-proton collisions, which could be described by a two-peak hypothesis
LHCb:2020pet . Assuming they decay directly to the $B^{\pm}K^{\mp}$ final
states, the two peaks could be associated to the resonances $B_{sJ}(6064)$ and
$B_{sJ}(6114)$, with the masses and widths as follows,

$\displaystyle B_{sJ}(6064):M$ $\displaystyle=$ $\displaystyle 6063.5\pm
1.2\pm 0.8~{}{\rm MeV},$ (1) $\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle 26\pm 4\pm 4~{}{\rm MeV};$ (2) $\displaystyle B_{sJ}(6114):M$
$\displaystyle=$ $\displaystyle 6114\pm 3\pm 5~{}{\rm MeV},$ (3)
$\displaystyle\Gamma$ $\displaystyle=$ $\displaystyle 66\pm 18\pm 21~{}{\rm
MeV}.$ (4)

However, if a decay through $B^{*\pm}K^{\mp}$ with a missing photon from the
$B^{*\pm}\to B^{\pm}\gamma$ decay is assumed, the masses and widths will shift
to be,

$\displaystyle B_{sJ}(6109):M$ $\displaystyle=$ $\displaystyle 6108.8\pm
1.1\pm 0.7~{}{\rm MeV},$ (5) $\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle 22\pm 5\pm 4~{}{\rm MeV};$ (6) $\displaystyle B_{sJ}(6158):M$
$\displaystyle=$ $\displaystyle 6158\pm 4\pm 5~{}{\rm MeV},$ (7)
$\displaystyle\Gamma$ $\displaystyle=$ $\displaystyle 72\pm 18\pm 25~{}{\rm
MeV}.$ (8)

These observations of the LHCb have enriched the bottom-strange spectrum. We
have tabulated the experimental information of all the $B_{s}$ mesons in Table
1.

Although these two states were suggested to be the $D$-wave orbital excited
bottom-strange mesons, their masses are significantly lower than the quark
model predictions (Lu:2016bbk, ; Godfrey:2016nwn, ; Ebert:2009ua, ;
Sun:2014wea, ). There have been some theoretical works to study these states
(Chen:2022fye, ; li:2021hss, ; Patel:2022hhl, ; Gandhi:2022nnk, ) (A recent
review about these two states can be found in Ref. Chen:2022asf ). Based on
the nonrelativistic linear potential model, Ref. (li:2021hss, ) explained the
$B_{sJ}(6114)$ as the $B_{s}(2^{3}S_{1})-B_{s}(1^{3}D_{1})$ mixing state, and
obtained the mass $M=6114$ MeV and width $\Gamma=95\pm 15$ MeV with the mixing
angle $\theta=-(45\pm 16)^{\circ}$, which is supported by the results based on
the screened potential model (Patel:2022hhl, ). In addition, Ref. (li:2021hss,
) also supports the interpretation of the $B_{sJ}(6114)$ as a pure
$B_{s}(1^{3}D_{1})$ state if there is small mixing between $B_{s}(2^{3}S_{1})$
and $B_{s}(1^{3}D_{1})$.

For the $B_{sJ}(6064)$, there are also two possible interpretations in Ref.
(li:2021hss, ). According to the first interpretation, the narrow structure
around 6064 MeV is mainly caused by the $B_{sJ}(6109)$ resonance, which is
regarded as the $B_{s}(1^{1}D_{2})-B_{s}(1^{3}D_{2})$ mixing state. For the
second interpretation, the $B_{sJ}(6064)$ could be explained as a pure
$B_{s}(1^{3}D_{3})$ with the predicted mass $M=6067$ MeV and width $\Gamma=13$
MeV.

In addition, authors of Ref. (Gandhi:2022nnk, ) have studied the $B$ and
$B_{s}$ mesons using the heavy quark effective theory, which suggests that the
$B_{sJ}(6064)$ could be the candidate of the $B_{s}(2^{3}S_{1})$ with the
predicted mass $6033.0\pm 2.4$ MeV and width $170\pm 1.5$ MeV, and the
$B_{sJ}(6114)$ could be the candidate of the $B_{s}(1^{3}D_{3})$ with
predicted mass $6247.0\pm 2.4$ MeV and width $82.0\pm 1.0$ MeV. However, the
predicted width for $B_{sJ}(6064)$ and mass for $B_{sJ}(6114)$ are larger than
experimental values.

Besides the conventional $B_{s}$ mesons explanation, the two states are also
regarded as $b\bar{s}q\bar{q}$ tetraquark states in Ref. (Chen:2021uou, ), or
as the $\bar{B}K^{*}$ molecular state with quantum numbers $(I)J^{P}=(0)1^{+}$
in Ref. (Kong:2021ohg, ). Thus, one can find that the natures of these two
states are still in debate, and more efforts are needed to shed light on their
internal structures.

As we known, although the quenched quark models have obtained lots of success
in the last decades, the effects of the sea quarks and gluons interactions are
not taken into account. The unquenched models, considering all kinds of
additional effects, have been developed, and widely used to describe the
hadron spectra, such as the coupled channel model (Ferretti:2013faa, ;
Ferretti:2013vua, ; Ortega:2009hj, ; Ortega:2016mms, ; Hao:2022vwt, ;
Yang:2023tvc, ; Xie:2021dwe, ) and the screened potential model (Song:2015fha,
; Song:2015nia, ; Wang:2018rjg, ; Wang:2019mhs, ; Hao:2019fjg, ; Feng:2022esz,
). In the quenched potential model, the potentials mainly contain a coulomb
term at short distances and the linear confining interaction at large
distances. However this is not appreciate in the large mass range, since the
linear potential, which is expected to be dominant in large mass region, will
be screened or softened by the vacuum polarization effects of dynamical
fermions (Born:1989iv, ; Li:2009zu, ), i.e., the unquenched effects reflecting
the sea quarks or gluons contributions to some extent. Clearly, the unquenched
effects can lead to important influence for higher radial and orbital excited
hadrons, which means that the predicted masses of the higher excited states
will be smaller than the ones of the general liner potential models. Comparing
with the coupled channel model, the screened potential model is simpler, and
have been successfully used to describe the spectra of the charmed-strange
meson (Song:2015nia, ; Gao:2022bsb, ), charmed meson (Song:2015fha, ), excited
$\rho$ mesons (Feng:2022hwq, ; Li:2021qgz, ), bottom mesons Feng:2022esz ,
charmonium (Wang:2019mhs, ; Li:2009zu, ; Hao:2019fjg, ), and bottomonium
(Wang:2018rjg, ; Li:2009nr, ).

Table 1:  Experimental information of the $B_{s}$ mesons (ParticleDataGroup:2022pth, ). state | mass (MeV) | width (MeV) | $I(J^{P})$
---|---|---|---
$B_{s}$ | $5366.92\pm 0.10$ | $-$ | $0(0^{-})$
$B_{s}^{*}$ | $5415.4^{+1.8}_{-1.5}$ | $-$ | $0(1^{-})$
$B_{s1}(5830)$ | $5828.70\pm 0.20$ | $0.5\pm 0.3\pm 0.3$ | $0(1^{+})$
$B_{s2}^{*}(5840)$ | $5839.86\pm 0.12$ | $1.49\pm 0.27$ | $0(2^{+})$
$B_{sJ}(6064)$ | $6063.5\pm 1.2\pm 0.8$ | $26\pm 4\pm 4$ |  $0(?^{?})$LHCb:2020pet
$B_{sJ}(6114)$ | $6114\pm 3\pm 5$ | $66\pm 18\pm 21$ |  $0(?^{?})$LHCb:2020pet
$B_{sJ}(6109)$ | $6108.8\pm 1.1\pm 0.7$ | $22\pm 5\pm 4$ |  $0(?^{?})$LHCb:2020pet
$B_{sJ}(6158)$ | $6158\pm 4\pm 5$ | $72\pm 18\pm 25$ |  $0(?^{?})$LHCb:2020pet

In this paper, we use the screened nonrelativistic quark model and the
${}^{3}P_{0}$ model to study the spectrum and the strong decay properties of
the $B_{s}$ mesons, and also to explore the possible assignments of the two
resonances recently observed by the LHCb Collaboration.

This article is organized as follows. In Sec. II, we give a brief introduction
about the screened nonrelativestic quark model and the ${}^{3}P_{0}$ model. In
Sec. III, the numerical results and the discussions are presented. Finally,
the summary is given in Sec. IV.

##  II Theoretical models

###  II.1 Screened nonrelativistic quark model

The nonrelativistic quark model mainly includes the confinement term , the
spin-dependent term, and the one-loop correction for the spin-dependent terms
Gupta:1981pd ; Pantaleone:1985uf ; Lakhina:2006fy , and the Hamiltonian for a
$q\bar{q}$ meson system is defined as (Li:2010vx, ; Lu:2016bbk, ),

$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{sd}+C_{q\bar{q}},$ (9)

where $\mathcal{H}_{0}$ is the zeroth-order Hamiltonian, $\mathcal{H}_{sd}$ is
the spin-dependent Hamiltonian, and $C_{q\bar{q}}$ is a constant, which will
be fixed to experimental data. The $\mathcal{H}_{0}$ can be compressed as,

$\displaystyle\mathcal{H}_{0}$ $\displaystyle=$
$\displaystyle\frac{\bm{p}^{2}}{M_{r}}-\frac{4}{3}\frac{{\alpha}_{s}}{r}+br+\frac{32{\alpha}_{s}{\sigma}^{3}e^{-{\sigma}^{2}r^{2}}}{9\sqrt{\pi}m_{q}m_{\bar{q}}}{\bm{S}}_{q}\cdot{\bm{S}}_{\bar{q}},$
(10)

where the confinement interaction includes the standard Coulomb potential
$-4\alpha_{s}/3r$ and the linear scalar potential $br$. The last term is the
hyperfine interaction that could be treated nonperturbatively. $\bm{p}$ is
quark momentum in the system of $q\bar{q}$ meson, $r=|\vec{r}\,|$ is the
$q\bar{q}$ separation, $M_{r}=2m_{q}m_{\bar{q}}/(m_{q}+m_{\bar{q}})$, $m_{q}$
($m_{\bar{q}}$) and $\bm{S}_{q}$ (${\bm{S}}_{\bar{q}}$) are the reduced mass
of the $q\bar{q}$ system, the mass and spin of the constituent quark $q$
(antiquark $\bar{q}$), respectively.

The spin-dependent term $\mathcal{H}_{sd}$ is,

$\displaystyle\mathcal{H}_{sd}$ $\displaystyle=$
$\displaystyle\left(\frac{\bm{S}_{q}}{2m_{q}^{2}}+\frac{{\bm{S}}_{\bar{q}}}{2m_{\bar{q}}^{2}}\right)\cdot\bm{L}\,\left(\frac{1}{r}\frac{dV_{c}}{dr}+\frac{2}{r}\frac{dV_{1}}{dr}\right)$
(11)
$\displaystyle+\frac{{\bm{S}}_{+}\cdot\bm{L}}{m_{q}m_{\bar{q}}}\left(\frac{1}{r}\frac{dV_{2}}{r}\right)$
$\displaystyle+\frac{3{\bm{S}}_{q}\cdot\hat{\bm{r}}\,{\bm{S}}_{\bar{q}}\cdot\hat{\bm{r}}-{\bm{S}}_{q}\cdot{\bm{S}}_{\bar{q}}}{3m_{q}m_{\bar{q}}}V_{3}$
$\displaystyle+\left[\left(\frac{{\bm{S}}_{q}}{m_{q}^{2}}-\frac{{\bm{S}}_{\bar{q}}}{m_{\bar{q}}^{2}}\right)+\frac{{\bm{S}}_{-}}{m_{q}m_{\bar{q}}}\right]\cdot\bm{L}V_{4},$

with

$\displaystyle V_{c}$ $\displaystyle=$
$\displaystyle-\frac{4}{3}\frac{{\alpha}_{s}}{r}+br,$ $\displaystyle V_{1}$
$\displaystyle=$ $\displaystyle-
br-\frac{2}{9\pi}\frac{{\alpha}_{s}^{2}}{r}\left[9\,{\rm
ln}(\sqrt{m_{q}m_{\bar{q}}}r)+9{\gamma}_{E}-4\right],$ $\displaystyle V_{2}$
$\displaystyle=$
$\displaystyle-\frac{4}{3}\frac{{\alpha}_{s}}{r}-\frac{1}{9\pi}\frac{{\alpha}_{s}^{2}}{r}\left[-18\,{\rm
ln}(\sqrt{m_{q}m_{\bar{q}}}r)+54\,{\rm ln}(\mu r)\right.$
$\displaystyle\left.+36{\gamma}_{E}+29\right],$ $\displaystyle V_{3}$
$\displaystyle=$
$\displaystyle-\frac{4{\alpha}_{s}}{r^{3}}-\frac{1}{3\pi}\frac{{\alpha}_{s}^{2}}{r^{3}}\left[-36\,{\rm
ln}(\sqrt{m_{q}m_{\bar{q}}}r)+54\,{\rm ln}(\mu r)\right.$
$\displaystyle\left.+18{\gamma}_{E}+31\right],$ $\displaystyle V_{4}$
$\displaystyle=$ $\displaystyle\frac{1}{\pi}\frac{{\alpha}_{s}^{2}}{r^{3}}{\rm
ln}\left(\frac{m_{\bar{q}}}{m_{q}}\right),$ (12)

where $\bm{S}_{\pm}={\bm{S}}_{q}\pm{\bm{S}}_{\bar{q}}$, $\bm{L}$ is the
relative orbital angular momentum of the $q\bar{q}$ system. We take Euler
constant $\gamma_{E}=0.5772$, the scalar $\mu=1$ GeV, ${\alpha}_{s}=0.53$,
$b=0.135$ GeV2, $\sigma=1.13$ GeV, $m_{u}=m_{d}=0.45$ GeV, $m_{s}=0.55$ GeV
and $m_{b}=4.5$ GeV (Lakhina:2006fy, ).

The screening effects are introduced by the following replacement,

$\displaystyle br\to V^{\text{scr}}(r)=\frac{b(1-e^{-\beta r})}{\beta},$ (13)

where $V^{\text{scr}}(r)$ behaves like $br$ at short distances and constant
$b/\beta$ at large distance Song:2015nia ; Song:2015fha , $\beta$ is the
parameter which controls the power of the screening effects. One can find that
the screened potential approximates the liner potential for a small distance
$r$, and will be softened for large distance $r$. Since the distance between
the quarks in the excited bottom-strange mesons is larger than the one of the
ground bottom-strange meson, it is expected that the spectrum of the screened
potential model is more sensitive for the excited bottom-strange mesons.

The spin-orbit term $\mathcal{H}_{sd}$ can be decomposed into symmetric part
$\mathcal{H}_{sym}$ and antisymmetric part $\mathcal{H}_{anti}$, which can be
expressed as (Lu:2016bbk, )

$\displaystyle\mathcal{H}_{sym}$ $\displaystyle=$
$\displaystyle\frac{{\bm{S}}_{+}\cdot{\bm{L}}}{2}\left[\left(\frac{1}{2m_{q}^{2}}+\frac{1}{2m_{\bar{q}}^{2}}\right)\left(\frac{1}{r}\frac{dV_{c}}{dr}+\frac{2}{r}\frac{dV_{1}}{dr}\right)\right.$
(14)
$\displaystyle\left.+\frac{2}{m_{q}m_{\bar{q}}}\left(\frac{1}{r}\frac{dV_{2}}{r}\right)+\left(\frac{1}{m_{q}^{2}}-\frac{1}{m_{\bar{q}}^{2}}\right)V_{4}\right],$
$\displaystyle\mathcal{H}_{anti}$ $\displaystyle=$
$\displaystyle\frac{{\bm{S}}_{-}\cdot{\bm{L}}}{2}\left[\left(\frac{1}{2m_{q}^{2}}-\frac{1}{2m_{\bar{q}}^{2}}\right)\left(\frac{1}{r}\frac{dV_{c}}{dr}+\frac{2}{r}\frac{dV_{1}}{dr}\right)\right.$
(15)
$\displaystyle\left.+\left(\frac{1}{m_{q}^{2}}+\frac{1}{m_{\bar{q}}^{2}}+\frac{2}{m_{q}m_{\bar{q}}}\right)V_{4}\right].$

The antisymmetric part $\mathcal{H}_{anti}$ gives rise to the the spin-orbit
mixing of the heavy-light mesons with different total spins but with the same
total angular momentum, such as $B_{s}(n{}^{3}L_{L})$ and
$B_{s}(n{}^{1}L_{L})$. Hence, the mixing of the two physical states
$B_{sL}(nL)$ and $B_{sL}^{\prime}(nL)$ can be expressed as,

$\left(\begin{array}[]{cr}B_{sL}(nL)\\\
B^{\prime}_{sL}(nL)\end{array}\right)=\left(\begin{array}[]{cr}\cos\theta_{nL}&\sin\theta_{nL}\\\
-\sin\theta_{nL}&\cos\theta_{nL}\end{array}\right)\left(\begin{array}[]{cr}B_{s}(n^{1}L_{L})\\\
B_{s}(n^{3}L_{L})\end{array}\right),$ (16)

where the $\theta_{nL}$ is the mixing angles.

With above formalism, one can solve the Schrödinger equation with Hamiltonian
$\mathcal{H}$ of Eq. (9) to obtain the mass spectrum and the meson wave
functions, where the wave functions will be used to calculate the strong
decays of excited bottom-strange mesons in the ${}^{3}P_{0}$ model.

###  II.2 The ${}^{3}P_{0}$ model

The ${}^{3}P_{0}$ model was proposed by Micu Micu:1968mk  and further
developed by Le Yaouanc LeYaouanc:1972vsx ; LeYaouanc:1974cvx ;
LeYaouanc:1977fsz ; LeYaouanc:1977gm , and it has been widely used to
calculate the OZI allowed decay processes Roberts:1992esl ; Barnes:1996ff ;
Barnes:2002mu ; Close:2005se ; Barnes:2005pb ; Li:2009rka ; Li:2009qu ;
Li:2010vx ; Lu:2014zua ; Pan:2016bac ; Lu:2016bbk ; Hao:2019fjg ; Feng:2022esz
; Wang:2018rjg ; Wang:2019mhs ; Song:2015nia ; Song:2015fha ; Xue:2018jvi ;
Wang:2017pxm . In this model, the meson decay occurs through the regroupment
between the $q\bar{q}$ of the initial meson and the another $q\bar{q}$ pair
created from vacuum with the quantum numbers $J^{PC}=0^{++}$. The transition
operator $\mathcal{T}$ of the decay $A\rightarrow BC$ in the ${}^{3}P_{0}$
model is given by

$\displaystyle\mathcal{T}=-3\gamma\sum\limits_{m}\langle 1m;1-m|00\rangle\int
d^{3}\bm{p}_{3}d^{3}\bm{p}_{4}\delta^{3}(\bm{p}_{3}+\bm{p}_{4})$
$\displaystyle\mathcal{Y}_{1m}\left(\frac{\bm{p}_{3}-\bm{p}_{4}}{2}\right)\chi^{34}_{1,-m}\phi^{34}_{0}\left(\omega^{34}_{0}\right)b^{{\dagger}}_{3}(\bm{p}_{3})d^{{\dagger}}_{4}(\bm{p}_{4}),$
(17)

where
${\cal{Y}}^{m}_{1}(\bm{p})\equiv|\bm{p}|^{1}Y^{m}_{1}(\theta_{p},\phi_{p})$ is
solid harmonic polynomial in the momentum space of the created quark-antiquark
pair. $\chi^{34}_{1,-m}$, $\phi^{34}_{0}$, and $\omega^{34}_{0}$ are the spin,
flavor, and color wave functions, respectively. The paramtere $\gamma$ is the
quark pair creation strength parameter for $u\bar{u}$ and $d\bar{d}$ pairs,
and for $s\bar{s}$ we take $\gamma_{s\bar{s}}=\gamma\frac{m_{u}}{m_{s}}$
LeYaouanc:1977gm . The parameter $\gamma$ can be determined by fitting to the
experimental data. The partial wave amplitude ${\cal{M}}^{LS}(\bm{P})$ of the
decay $A\rightarrow BC$ is be given by,

$\displaystyle{\cal{M}}^{LS}(\bm{P})$ $\displaystyle=$
$\displaystyle\sum_{M_{J_{B}},M_{J_{C}},M_{S},M_{L}}\langle
LM_{L}SM_{S}|J_{A}M_{J_{A}}\rangle$ (18) $\displaystyle\langle
J_{B}M_{J_{B}}J_{C}M_{J_{C}}|SM_{S}\rangle$ $\displaystyle\int
d\Omega\,\mbox{}Y^{\ast}_{LM_{L}}{\cal{M}}^{M_{J_{A}}M_{J_{B}}M_{J_{C}}}(\bm{P}),$

where ${\cal{M}}^{M_{J_{A}}M_{J_{B}}M_{J_{C}}}(\bm{P})$ is the helicity
amplitude,

$\displaystyle\langle
BC|T|A\rangle=\delta^{3}(\bm{P}_{A}-\bm{P}_{B}-\bm{P}_{C})$
$\displaystyle{\cal{M}}^{M_{J_{A}}M_{J_{B}}M_{J_{C}}}(\bm{P}).$ (19)

Here, $|A\rangle$, $|B\rangle$, and $|C\rangle$ denote the mock meson states.
Then, the decay width $\Gamma(A\rightarrow BC)$ can be expressed as

$\displaystyle\Gamma(A\rightarrow BC)=\frac{\pi
P}{4M^{2}_{A}}\sum_{LS}|{\cal{M}}^{LS}(\bm{P})|^{2},$ (20)

where
$P=|\bm{P}|=\frac{\sqrt{[M^{2}_{A}-(M_{B}+M_{C})^{2}][M^{2}_{A}-(M_{B}-M_{C})^{2}]}}{2M_{A}}$,
$M_{A}$, $M_{B}$, and $M_{C}$ are the masses of the mesons $A$, $B$, and $C$,
respectively. The spatial wave functions of the mesons in the ${}^{3}P_{0}$
model are obtained by solving the Schr$\ddot{o}$dinger equation in Eq. (9).

##  III Results and discussions

In the calculation, the screened parameter $\beta=0.025$ GeV and the constant
$C_{q\bar{q}}=0.1035$ GeV were obtained by fitting the well-known states
$B_{s}(1^{1}S_{0})$, $B_{s}^{*}(1^{3}S_{1})$, and
$B_{s2}^{*}(5840)(1^{3}P_{2})$. The other parameters are taken from Ref.
Lakhina:2006fy . The ${}^{3}P_{0}$ model parameter $\gamma=0.354$ is obtained
by fitting the total decay width of the $B_{s2}^{*}(5840)$, which is regarded
as the $B_{s}(1^{3}P_{2})$. With these parameters, the predicted ratio of the
$B_{s2}^{*}(5840)$ decay modes,

$\frac{\Gamma(B_{s2}^{*}(5840)\to B^{*+}K^{-})}{\Gamma(B_{s2}^{*}(5840)\to
B^{+}K^{-})}=0.095,$ (21)

which is consistent with LHCb experimental data of $0.093\pm 0.013\pm 0.012$
LHCb:2012iuq .

The predicted masses of the $B_{s}$ mesons are listed in Table 2, where we
also show the predictions of other theoretical works for comparison. The
$B_{s}$[$B_{s}(1^{1}S_{0})$], $B_{s}^{*}$[$B_{s}(1^{3}S_{1})$], and
$B_{s2}^{*}(5840)$[$B_{s}(1^{3}P_{2})$] can be well described in the spectrum.
The masses of the $B_{s}(1^{3}D_{1})$ and $B_{s}(1^{3}D_{3})$ are predicted to
be 6061 MeV and 6117 MeV, in good agreement with the experimental results of
the $B_{sJ}(6064)$ ($6063.5\pm 1.2\pm 0.8$ MeV) and $B_{sJ}(6114)$ ($6114\pm
3\pm 5$ MeV), respectively, which indicates that the two states could be the
possible candidates of the $B_{s}(1^{3}D_{1})$ and $B_{s}(1^{3}D_{3})$ .

Of course, only the mass information is not enough to establish these
assignments. We also calculate the strong decay widths of the $B_{s}$ mesons,
as shown in Table 3. One can find that the predicted width of 23 MeV for
$B_{s}(1^{3}D_{3})$ is in good agreement with the measured width $26\pm 4\pm
4$ MeV of the $B_{sJ}(6064)$, which supports the $B_{s}(1^{3}D_{3})$
assignment of the $B_{sJ}(6064)$. In addition, the predicted width for
$B_{s}(1^{3}D_{1})$ is 127 MeV, reasonably consistent with the one of
$B_{sJ}(6114)$ if taking into account the large experimental uncertainties.
Thus, the $B_{sJ}(6114)$ could be explained as the $B_{s}(1^{3}D_{1})$ state,
and more precise measurements will be helpful to pin down this assignment.

As we discussed in the introduction, if a decay through $B^{*\pm}K^{\mp}$ with
a missing photon from the $B^{*\pm}\to B^{\pm}\gamma$ decay is assumed, the
masses and widths of the two states observed by LHCb will shift, and the two
states are named as $B_{sJ}(6109)$ and $B_{sJ}(6158)$ LHCb:2020pet . In this
case, the mass and width of the $B_{sJ}(6109)$ are close to the predicted mass
(6132 MeV) and width (43 MeV) of $B_{s2}^{\prime}(1D)$, which implies that
$B_{sJ}(6109)$ could be regarded as the $B_{s2}^{\prime}(1D)$ state. On the
other hand, the mass and width of the $B_{sJ}(6158)$ are close to the
predicted mass (6194 MeV) and width (75 MeV) of $B_{s1}(2P)$, respectively,
which supports the assignment of the $B_{sJ}(6158)$ as the $B_{s1}(2P)$ state.
It should be stressed that the two solutions could be not distinguished
according to the present LHCb measurements. Thus, the future precise
measurements of their masses, widths, and the quantum numbers of the spin-
parity would be helpful to shed light on this problem, and deepen our
understanding the spectra of the bottom-strange mesons.

Table 2:  The mass spectrum of the $B_{s}$ mesons predicted by the screened quark model in units of MeV. The mixing angles of $B_{sL}-B^{\prime}_{sL}$ calculated in this work are $\theta_{1P}=-55.8^{\circ}$, $\theta_{2P}=-54.3^{\circ}$ and $\theta_{1D}=-50.3^{\circ}$. State | PDG | Ours |  NRLu:2016bbk  |  GIGodfrey:2016nwn  |  DRVEbert:2009ua  |  KAGandhi:2022nnk  |  VRAPatel:2022hhl
---|---|---|---|---|---|---|---
$B_{s}(1^{1}S_{0})$ | $5366.92\pm 0.10$ | 5367 | 5362 | 5394 | 5372 |  | 5359
$B_{s}(1^{3}S_{1})$ | $5415.4^{+1.8}_{-1.5}$ | 5419 | 5413 | 5450 | 5414 |  | 5415
$B_{s}(2^{1}S_{0})$ |  | 5947 | 5977 | 5984 | 5976 | 6025 | 5980
$B_{s}(2^{3}S_{1})$ |  | 5972 | 6003 | 6012 | 5992 | 6033 | 5993
$B_{s}(1^{3}P_{0})$ |  | 5753 | 5756 | 5831 | 5833 | 5709 | 5798
$B_{s1}(1P)$ |  | 5797 | 5801 | 5857 | 5831 | 5768 | 5818
$B^{\prime}_{s1}(1P)$ | $5828.70\pm 0.20$ | 5825 | 5836 | 5861 | 5865 | 5875 | 5846
$B_{s}(1^{3}P_{2})$ | $5839.86\pm 0.12$ | 5840 | 5851 | 5876 | 5842 | 5890 | 5838
$B_{s}(2^{3}P_{0})$ |  | 6158 | 6203 | 6279 | 6318 | 6387 | 6292
$B_{s1}(2P)$ | $6158\pm 4\pm 5$ | 6194 | 6241 | 6279 | 6321 | 6393 | 6304
$B^{\prime}_{s1}(2P)$ |  | 6236 | 6297 | 6296 | 6345 | 6470 | 6320
$B_{s}(2^{3}P_{2})$ |  | 6249 | 6309 | 6295 | 6359 | 6476 | 6316
$B_{s}(1^{3}D_{1})$ | $6114\pm 3\pm 5$ | 6117 | 6142 | 6182 | 6209 | 6247 | 6144
$B_{s2}(1D)$ |  | 6053 | 6087 | 6169 | 6189 | 6256 | 6139
$B^{\prime}_{s2}(1D)$ | $6108.8\pm 1.1\pm 0.7$ | 6132 | 6159 | 6196 | 6218 | 6292 | 6135
$B_{s}(1^{3}D_{3})$ | $6063.5\pm 1.2\pm 0.8$ | 6061 | 6096 | 6179 | 6191 | 6297 | 6139
Table 3:  The decay widths of the $B_{s}$ mesons in units of MeV. The label ‘$-$’ means that the channel is forbidden or there is no experimental information. State |  | $BK$ | $B^{*}K$ | $BK^{*}$ | $B^{*}K^{*}$ | $B_{s}\eta$ | $B_{s}^{*}\eta$ | Total | Exp.
---|---|---|---|---|---|---|---|---|---
$B_{s1}(1P)$ | $-$ | $-$ | $-$ | $-$ | $-$ | $-$ | $-$ | $-$ | $-$
$B^{\prime}_{s1}(1P)$ | $B_{s1}(5830)$ | $-$ | 0.04 | $-$ | $-$ | $-$ | $-$ | 0.03 | $0.5\pm 0.3\pm 0.3$
$B_{s}(1^{3}P_{2})$ | $B_{s2}^{*}(5840)$ | 1.4 | 0.1 | $-$ | $-$ | $-$ | $-$ | 1.5 | $1.49\pm 0.27$
$B_{s}(2^{3}P_{0})$ | $-$ | 61.5 | $-$ | $-$ | $-$ | 5.3 | $-$ | 66.8 | $-$
$B_{s1}(2P)$ | $B_{sJ}(6158)$ | $-$ | 58.2 | 13.2 | $-$ | $-$ | 3.9 | 75.3 | $72\pm 18\pm 25$
$B^{\prime}_{s1}(2P)$ | $-$ | $-$ | 15.0 | 41.8 | 52.4 | $-$ | 4.8 | 113.9 | $-$
$B_{s}(2^{3}P_{2})$ | $-$ | 1.7 | 7.3 | 18.7 | 139.7 | 1.3 | 2.8 | 171.6 | $-$
$B_{s}(1^{3}D_{1})$ | $B_{sJ}(6114)$ | 67.5 | 37.4 | $-$ | $-$ | 15.3 | 6.8 | 127.0 | $66\pm 18\pm 21$
$B_{s2}(1D)$ | $-$ | $-$ | 115.5 | $-$ | $-$ | $-$ | 11.7 | 127.3 | $-$
$B^{\prime}_{s2}(1D)$ | $B_{sJ}(6109)$ | $-$ | 41.8 | $-$ | $-$ | $-$ | 1.2 | 43.0 | $22\pm 5\pm 4$
$B_{s}(1^{3}D_{3})$ | $B_{sJ}(6064)$ | 12.0 | 10.6 | $-$ | $-$ | 0.3 | 0.1 | 23.0 | $26\pm 4\pm 4$
---
Figure 1: Total decay widths of the $B_{s1}(5830)$ as the $1P^{\prime}$ versus
the mixing angle. The blue band denotes the uncertainties of the experimental
data.

According to the predicted mass spectrum of Table 2, the predicted mass 5840
MeV of the $B_{s1}^{\prime}(1P)$ is very close to the one of the
$B_{s1}(5830)$. With the mixing angle of $-55.8^{\circ}$, the decay widths of
the $B_{s1}^{\prime}(1P)$ are calculated, as shown in Table 3, and the
predicted total decay width is very small, in good agreement with the
experimental value $0.5\pm 0.3\pm 0.3$ MeV of the $B_{s1}(5830)$. We show the
total decay widths of the $B_{s1}^{\prime}(1P)$ versus the mixing angle in
Fig. 1, where one can find the total decay width is still consistent with the
experimental data with the mixing angle in the range of
$-60^{\circ}\sim-50^{\circ}$.

In addition, we also predict the decay widths of the other excited $B_{s}$
mesons with the predicted masses of Table 2. For the $B_{s}(2^{3}P_{0})$, its
total decal width is predicted to be 66.8 MeV, and the dominant decay mode is
$BK$ with the branching fraction 92%. The total decay width of
$B^{\prime}_{s1}(2P)$ is 113.9 MeV, while the dominant decay modes are
$BK^{*}/B^{*}K^{*}$. The predicted width and the dominant decay mode of the
$B_{s}(2^{3}P_{2})$ are 171.6 MeV and $B^{*}K^{*}$, respectively. In addition,
the decay widths of the $B_{s2}(1D)$ are 127.3 MeV, and the dominant decay
mode is $B^{*}K$. Our results should be helpful to search for them in
experiments, such as LHC.

##  IV Summary

Recently, the LHCb Collaboration has observed two resonances $B_{sJ}(6064)$
and $B_{sJ}(6114)$ assuming they decay directly to the $BK$ final states.
However, their masses and widths will shift if a decay through $B^{*}K$ with a
missing photon from the $B^{*}$ decay, and those two resonances are named as
$B_{sJ}(6109)$ and $B_{sJ}(6158)$.

Motivated by the recently LHCb measurements, in this paper we calculate the
spectrum of the $B_{s}$ mesons within the screened nonrelativistic quark
model, and also investigate the strong decay properties of these mesons with
the ${}^{3}P_{0}$ model.

By comparing with the experimental data, it is found that the $B_{sJ}(6064)$
and $B_{sJ}(6114)$ states, as the first solution of the recently LHCb
measurements, could be explained as the $B_{s}(1^{3}D_{3})$ and
$B_{s}(1^{3}D_{1})$, respectively. In addition, the $B_{sJ}(6109)$ and
$B_{sJ}(6158)$ states, as the second solution of the LHCb measurements, could
be explained as the $B^{\prime}_{s2}(1D)$ and $B_{s1}(2P)$, respectively.
Since those two solutions could be not distinguished according to the present
LHCb measurements, thus the future precise measurements of their masses,
widths, and the quantum numbers of the spin-parity would be helpful to shed
light on this problem, and deepen our understanding the spectra of the bottom-
strange mesons.

In addition, we suggest that the state $B_{s1}(5830)$ could be explained as
the candidate of the $B_{s1}^{\prime}(1P)$ state. With the predicted masses of
the other excited $B_{s}$ states, we also predict their decay widths as well
as the dominant decay modes, which should be helpful for experiments to search
for them.

It should be stressed that there are already many theoretical studies about
the family of the bottom-strange mesons. Comparing with those works, we have
adopted the nonrelativistic quark model by taking into account the screening
effects, which play an important role for the higher radial and orbital
excited mesons, and our results could give better descriptions for all the
existed bottom-strange mesons.

###### Acknowledgements.

This work is partly supported by the Natural Science Foundation of Henan under
Grand Nos. 222300420554 and 232300421140, the Project of Youth Backbone
Teachers of Colleges and Universities of Henan Province (2020GGJS017), the
Youth Talent Support Project of Henan (2021HYTP002), and the Open Project of
Guangxi Key Laboratory of Nuclear Physics and Nuclear Technology,
No.NLK2021-08.

## References

  * [1] Roel Aaij et al.  Observation of new excited $B^{0}_{s}$ states.  Eur. Phys. J. C, 81(7):601, 2021.
  * [2] Qi-Fang Lü, Ting-Ting Pan, Yan-Yan Wang, En Wang, and De-Min Li.  Excited bottom and bottom-strange mesons in the quark model.  Phys. Rev. D, 94(7):074012, 2016.
  * [3] Stephen Godfrey, K. Moats, and E. S. Swanson.  $B$ and $B_{s}$ Meson Spectroscopy.  Phys. Rev. D, 94(5):054025, 2016.
  * [4] D. Ebert, R. N. Faustov, and V. O. Galkin.  Heavy-light meson spectroscopy and Regge trajectories in the relativistic quark model.  Eur. Phys. J. C, 66:197–206, 2010.
  * [5] Yuan Sun, Qin-Tao Song, Dain-Yong Chen, Xiang Liu, and Shi-Lin Zhu.  Higher bottom and bottom-strange mesons.  Phys. Rev. D, 89(5):054026, 2014.
  * [6] Bing Chen, Si-Qiang Luo, Ke-Wei Wei, and Xiang Liu.  $b$-hadron spectroscopy study based on the similarity of double bottom baryon and bottom meson.  Phys. Rev. D, 105(7):074014, 2022.
  * [7] Qi li, Ru-Hui Ni, and Xian-Hui Zhong.  Towards establishing an abundant $B$ and $B_{s}$ spectrum up to the second orbital excitations.  Phys. Rev. D, 103:116010, 2021.
  * [8] Vikas Patel, Raghav Chaturvedi, and A. K. Rai.  Spectroscopic Properties of $\mathit{B}$ and $\mathit{B_{s}}$ meson using Screened Potential.  arXiv:2201.01120.
  * [9] Keval Gandhi and Ajay Kumar Rai.  Study of $B$, $B_{s}$ mesons using heavy quark effective theory.  Eur. Phys. J. C, 82(9):777, 2022.
  * [10] Hua-Xing Chen, Wei Chen, Xiang Liu, Yan-Rui Liu, and Shi-Lin Zhu.  An updated review of the new hadron states.  arXiv:2204.02649.
  * [11] Xiao yun Chen, Yue Tan, and Yuan Chen.  Study on $Z_{cs}$ and excited $B_{s}^{0}$ states in the chiral quark model.  Phys. Rev. D, 104(1):014017, 2021.
  * [12] Shu-Yi Kong, Jun-Tao Zhu, Dan Song, and Jun He.  Heavy-strange meson molecules and possible candidates $D_{s0}^{*}(2317)$, $D_{s1}(2460)$, and $X_{0}(2900)$.  Phys. Rev. D, 104(9):094012, 2021.
  * [13] J. Ferretti, G. Galatà, and E. Santopinto.  Interpretation of the $X(3872)$ as a charmonium state plus an extra component due to the coupling to the meson-meson continuum.  Phys. Rev. C, 88(1):015207, 2013.
  * [14] J. Ferretti and E. Santopinto.  Higher mass bottomonia.  Phys. Rev. D, 90(9):094022, 2014.
  * [15] P. G. Ortega, J. Segovia, D. R. Entem, and F. Fernandez.  Coupled channel approach to the structure of the $X(3872)$.  Phys. Rev. D, 81:054023, 2010.
  * [16] Pablo G. Ortega, Jorge Segovia, David R. Entem, and Francisco Fernandez.  Molecular components in P-wave charmed-strange mesons.  Phys. Rev. D, 94(7):074037, 2016.
  * [17] Wei Hao, Yu Lu, and Bing-Song Zou.  Coupled channel effects for the charmed-strange mesons.  Phys. Rev. D, 106(7):074014, 2022.
  * [18] Jing-Jing Yang, Wei Hao, Xiaoyu Wang, De-Min Li, Yu-Xiao Li, and En Wang.  The mass spectrum and strong decay properties of the charmed-strange mesons within Godfrey-Isgur model considering the coupled-channel effects.  3 2023.
  * [19] Jia-Ming Xie, Ming-Zhu Liu, and Li-Sheng Geng.  Ds0(2590) as a dominant cs¯ state with a small D*K component.  Phys. Rev. D, 104(9):094051, 2021.
  * [20] Qin-Tao Song, Dian-Yong Chen, Xiang Liu, and Takayuki Matsuki.  Higher radial and orbital excitations in the charmed meson family.  Phys. Rev. D, 92(7):074011, 2015.
  * [21] Qin-Tao Song, Dian-Yong Chen, Xiang Liu, and Takayuki Matsuki.  Charmed-strange mesons revisited: mass spectra and strong decays.  Phys. Rev. D, 91:054031, 2015.
  * [22] Jun-Zhang Wang, Zhi-Feng Sun, Xiang Liu, and Takayuki Matsuki.  Higher bottomonium zoo.  Eur. Phys. J. C, 78(11):915, 2018.
  * [23] Jun-Zhang Wang, Dian-Yong Chen, Xiang Liu, and Takayuki Matsuki.  Constructing $J/\psi$ family with updated data of charmoniumlike $Y$ states.  Phys. Rev. D, 99(11):114003, 2019.
  * [24] Wei Hao, Guan-Ying Wang, En Wang, Guan-Nan Li, and De-Min Li.  Canonical interpretation of the $X(4140)$ state within the ${}^{3}P_{0}$ model.  Eur. Phys. J. C, 80(7):626, 2020.
  * [25] Xue-Chao Feng, Wei Hao, and li-Juan Liu.  The assignments of the bottom mesons within the screened potential model and ${}^{3}P_{0}$ model.  Int. J. Mod. Phys. E, 31(07):2250066, 2022.
  * [26] K. D. Born, E. Laermann, N. Pirch, T. F. Walsh, and P. M. Zerwas.  Hadron Properties in Lattice QCD With Dynamical Fermions.  Phys. Rev. D, 40:1653–1663, 1989.
  * [27] Bai-Qing Li and Kuang-Ta Chao.  Higher Charmonia and X,Y,Z states with Screened Potential.  Phys. Rev. D, 79:094004, 2009.
  * [28] Zhuo Gao, Guan-Ying Wang, Qi-Fang Lü, Jingya Zhu, and Gao-Feng Zhao.  Canonical interpretation of the Ds0(2590)+ resonance.  Phys. Rev. D, 105(7):074037, 2022.
  * [29] Xue-Chao Feng, Zheng-Ya Li, De-Min Li, Qin-Tao Song, En Wang, and Wen-Cheng Yan.  Mass spectra and decay properties of the higher excited $\rho$ mesons.  Phys. Rev. D, 106(7):076012, 2022.
  * [30] Zheng-Ya Li, De-Min Li, En Wang, Wen-Cheng Yan, and Qin-Tao Song.  Assignments of the $Y(2040)$, $\rho$(1900), and $\rho$(2150) in the quark model.  Phys. Rev. D, 104(3):034013, 2021.
  * [31] Bai-Q Li and Kuang-Ta Chao.  Bottomonium Spectrum with Screened Potential.  Commun. Theor. Phys., 52:653–661, 2009.
  * [32] R. L. Workman et al.  Review of Particle Physics.  PTEP, 2022:083C01, 2022.
  * [33] S. N. Gupta and S. F. Radford.  Quark Quark and Quark - Anti-quark Potentials.  Phys. Rev. D, 24:2309–2323, 1981.
  * [34] James T. Pantaleone, S. H. Henry Tye, and Yee Jack Ng.  Spin Splittings in Heavy Quarkonia.  Phys. Rev. D, 33:777, 1986.
  * [35] Olga Lakhina and Eric S. Swanson.  A Canonical $D_{s}(2317)$?  Phys. Lett. B, 650:159–165, 2007.
  * [36] De-Min Li, Peng-Fei Ji, and Bing Ma.  The newly observed open-charm states in quark model.  Eur. Phys. J. C, 71:1582, 2011.
  * [37] L. Micu.  Decay rates of meson resonances in a quark model.  Nucl. Phys. B, 10:521–526, 1969.
  * [38] A. Le Yaouanc, L. Oliver, O. Pene, and J. C. Raynal.  Naive quark pair creation model of strong interaction vertices.  Phys. Rev. D, 8:2223–2234, 1973.
  * [39] A. Le Yaouanc, L. Oliver, O. Pene, and J. C. Raynal.  Resonant Partial Wave Amplitudes in $\pi n\to\pi\pi n$ According to the Naive Quark Pair Creation Model.  Phys. Rev. D, 11:1272, 1975.
  * [40] A. Le Yaouanc, L. Oliver, O. Pene, and J. C. Raynal.  Strong Decays of $\psi^{\prime\prime}(4.028)$ as a Radial Excitation of Charmonium.  Phys. Lett. B, 71:397–399, 1977.
  * [41] A. Le Yaouanc, L. Oliver, O. Pene, and J. C. Raynal.  Why Is $\psi^{\prime\prime\prime}(4.414)$ SO Narrow?  Phys. Lett. B, 72:57–61, 1977.
  * [42] W. Roberts and B. Silvestre-Brac.  General method of calculation of any hadronic decay in the ${}^{3}P_{0}$ model.  Few Body Syst., 11(4):171–193, 1992.
  * [43] Ted Barnes, F. E. Close, P. R. Page, and E. S. Swanson.  Higher quarkonia.  Phys. Rev. D, 55:4157–4188, 1997.
  * [44] T. Barnes, N. Black, and P. R. Page.  Strong decays of strange quarkonia.  Phys. Rev. D, 68:054014, 2003.
  * [45] F. E. Close and E. S. Swanson.  Dynamics and decay of heavy-light hadrons.  Phys. Rev. D, 72:094004, 2005.
  * [46] T. Barnes, S. Godfrey, and E. S. Swanson.  Higher charmonia.  Phys. Rev. D, 72:054026, 2005.
  * [47] De-Min Li and En Wang.  Canonical interpretation of the $\eta_{2}(1870)$.  Eur. Phys. J. C, 63:297–304, 2009.
  * [48] De-Min Li and Bing Ma.  Implication of BaBar’s new data on the $D_{s1}(2710)$ and $D_{sJ}(2860)$.  Phys. Rev. D, 81:014021, 2010.
  * [49] Qi-Fang Lü and De-Min Li.  Understanding the charmed states recently observed by the LHCb and BaBar Collaborations in the quark model.  Phys. Rev. D, 90(5):054024, 2014.
  * [50] Ting-Ting Pan, Qi-Fang Lü, En Wang, and De-Min Li.  Strong decays of the $X(2500)$ newly observed by the BESIII Collaboration.  Phys. Rev. D, 94(5):054030, 2016.
  * [51] Shi-Chen Xue, Guan-Ying Wang, Guan-Nan Li, En Wang, and De-Min Li.  The possible members of the $5^{1}S_{0}$ meson nonet.  Eur. Phys. J. C, 78(6):479, 2018.
  * [52] Guan-Ying Wang, Shi-Chen Xue, Guan-Nan Li, En Wang, and De-Min Li.  Strong decays of the higher isovector scalar mesons.  Phys. Rev. D, 97(3):034030, 2018.
  * [53] R Aaij et al.  First observation of the decay $B_{s2}^{*}(5840)^{0}\to B^{*+}K^{-}$ and studies of excited $B^{0}_{s}$ mesons.  Phys. Rev. Lett., 110(15):151803, 2013.


# Bayesian noise wave calibration for 21-cm global experiments

I. L. V. Roque,1 W. J. Handley1,2 and N. Razavi-Ghods1
1Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge, CB3
0HE, UK
2Kavli Institute for Cosmology, Madingley Road, Cambridge, CB3 0HA, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Detection of millikelvin-level signals from the ‘Cosmic Dawn’ requires an
unprecedented level of sensitivity and systematic calibration. We report the
theory behind a novel calibration algorithm developed from the formalism
introduced by the EDGES collaboration for use in 21-cm experiments.
Improvements over previous approaches are provided through the incorporation
of a Bayesian framework and machine learning techniques such as the use of
Bayesian evidence to determine the level of frequency variation of calibration
parameters that is supported by the data, the consideration of correlation
between calibration parameters when determining their values and the use of a
conjugate-prior based approach that results in a fast algorithm for
application in the field. In self-consistency tests using empirical data
models of varying complexity, our methodology is used to calibrate a 50
$\Omega$ ambient-temperature load. The RMS error between the calibration
solution and the measured temperature of the load is 8 mK, well within the
$1\sigma$ noise level. Whilst the methods described here are more applicable
to global 21-cm experiments, they can easily be adapted and applied to other
applications, including telescopes such as HERA and the SKA.

###### keywords:

instrumentation: interferometers – methods: statistical – dark ages,
reionization, first stars

††pubyear: 2020††pagerange: Bayesian noise wave calibration for 21-cm global
experiments–Bayesian noise wave calibration for 21-cm global experiments

##  1 Introduction

For nearly a century, scientists have been using radio-frequency instruments
to advance the study of astronomy and complement information from the visual
regime of the electromagnetic spectrum (Pritchard & Loeb, 2012). As we begin
to take measurements of the early universe, these instruments must continue to
evolve to support observations. Unexplored cosmic information from the Epoch
of Reionisation and Cosmic Dawn redshifted into the radio spectrum could
provide constraints on fundamental physics such as primordial black holes,
galaxy formation, and universal curvature as discussed in Furlanetto et al.
(2009). A unique probe of phenomena from the early cosmos is the hydrogen that
inundates the intergalactic medium (IGM). Heating and cooling of the IGM
associated with hydrogen’s absorption and emission of 21-cm photons produce a
dynamic brightness temperature relative to the cosmic microwave background
temperature, tracing the evolution of surrounding structure during the Cosmic
Dawn. The brightness temperature of this 21-cm photon signal can be described
by

$\displaystyle T_{21}(z)\approx$ $\displaystyle\ 0.023\mathrm{K}\ \times$ (1)
$\displaystyle
x_{\text{H\,{i}}}(z)\left[\left(\frac{0.15}{\Omega_{\mathrm{m}}}\right)\left(\frac{1+z}{10}\right)\right]^{\frac{1}{2}}\left(\frac{\Omega_{\mathrm{b}}h}{0.02}\right)\left[1-\frac{T_{\mathrm{R}}(z)}{T_{\mathrm{S}}(z)}\right],$

which is heavily dependent on environmental factors of the early universe such
as $x_{\text{H\,{i}}}$, the fraction of neutral hydrogen,
$\Omega_{\mathrm{m}}$ and $\Omega_{\mathrm{b}}$, the matter and baryon
densities with respect to the universal critical density for a flat universe
and Hubble’s constant. Here, the $0.023$ is a constant from atomic-line
physics. $T_{\mathrm{R}}$ is the background radiation temperature and
$T_{\mathrm{S}}$ is known as the ‘21-cm spin temperature’, which is related to
the kinetic temperature of neutral hydrogen gas in the IGM (Zaldarriaga et
al., 2004; Field, 1958). This cosmic hydrogen signature measurable in the
spectral sky has been redshifted to wavelengths under 200 MHz through the
expansion of the universe as discussed in Pritchard & Loeb (2012).

There has been a recent surge in the field of 21-cm cosmology following the
reported detection of an absorption feature consistent with a Cosmic Dawn
signature. This was reported by the Experiment to Detect the Global EoR
Signature (EDGES) in early 2018 from measurements of a sky-averaged radio
spectrum (Bowman et al., 2018). The signal, centred at 78 MHz with a width
corresponding to a period between 180 million and 270 million years after the
Big Bang, matches the theoretical position in frequency, but its depth of
$\sim 0.5$ K is a factor of two greater than the largest predictions from
theoretical models (Cohen et al., 2017). This discrepancy would suggest that
the temperature difference between the IGM and the cosmic microwave background
was much larger than previously thought and would require new physics to
explain, such as dark matter-baryon interactions (Barkana, 2018) or excess
radio backgrounds (Ewall-Wice et al., 2020).

Another possible explanation for this discrepancy is that the measured signal
is not cosmological but of systematic origin. This may be the case in EDGES
due to some of the methodology used, such as a potentially unphysical
foreground removal method and calibration of the receiver in a separate
environment from the data acquisition (Hills et al., 2018; Razavi-Ghods,
2018). In this paper, we present a novel calibration algorithm that improves
on the work of the EDGES team (Rogers & Bowman, 2012) through the utilisation
of a Bayesian framework to promote efficient use of the data to remove
systematics. Using conjugate priors and machine learning techniques, our
pipeline can be applied in the field with the collection of data with
additional capabilities for optimising individual noise wave parameters and
incorporating correlations between them.

This paper is organised as follows. In section 2 we review the methodology
behind calibration using noise waves as well as present a Bayesian framework
that provides greater flexibility in radiometer calibration. Section 3
describes the process of using mock data sets modelled after empirical
measurements of reflection coefficients with the incorporation of a realistic
noise model to evaluate our pipeline.

##  2 Methods

In this section, we detail the methodology behind radiometer calibration using
noise wave parameters. An overview of global signal measurement are outlined
in section 2.1. Section 2.2 summarises the basic procedure with some
mathematical improvements while section 2.3 describes our Bayesian framework
and its associated advantages.

###  2.1 Measuring the global signal

The noise necessitating calibration emerges during measurement-taking. In an
averaged or global experiment, the sky temperature
$T_{\mathrm{sky}}(\Omega,\nu,t)$ is a function of the direction $\Omega$,
frequency $\nu$ and time $t$. This can be broken down into two primary
components: the global 21-cm signal $T_{21}$ and astrophysical foregrounds
$T_{\mathrm{f}}$

$T_{\mathrm{sky}}(\Omega,\nu,t)=T_{21}(\nu)+T_{\mathrm{f}}(\Omega,\nu,t).$ (2)

The antenna measures the sky signal convolved with the normalised antenna
directivity $B$. The process of measurement introduces the random noise term
$N_{\mathrm{data}}$.

$D(\nu,t)=\int
T_{\mathrm{sky}}(\Omega,\nu,t)B(\Omega,\nu)\mathrm{d}\Omega+N_{\mathrm{data}}.$
(3)

Our 21-cm signature can thus be represented as

$T_{21}\approx D(\nu,t)-\int
T_{\mathrm{f}}(\Omega,\nu,t)B(\Omega,\nu)\mathrm{d}\Omega-N_{\mathrm{data}}.$
(4)

Here, the integral is assessed through foreground and beam modelling
techniques such as those discussed in Anstey et al. (2020) while modelling of
$N_{\mathrm{data}}$ from the statistical properties of $D(\nu,t)$ is
accomplished by a calibration algorithm as articulated in this paper and
outlined in Fig. 1. Having a fully Bayesian framework when modelling the beam,
the sky and the systematics has major advantages for global 21-cm experiments
such as REACH (de Lera Acedo et al., 2021), as it provides the greatest
flexibility in being able to model all effects and jointly fit for them.

Figure 1: Diagram showing the evolution of the 21-cm signal hampered by
astrophysical foregrounds, convolvution with the antenna beam and the
emergence of measurement noise before calibration to retrieve the antenna
temperature.

###  2.2 Calibration methodology

The standard calibration strategy follows the method introduced by Dicke to
characterise systematic features in radio frequency instruments (Dicke, 1946)
and is widely used in experiments such as EDGES (Monsalve et al., 2017) and
LOFAR (Bilous, A. V. et al., 2016) to evaluate the spectral index of the sky’s
diffuse radio background (Rogers & Bowman, 2012). This technique involves
measurements of two internal reference standards; a load and a noise source,
in addition to a series of external calibration sources attached to the
receiver input in lieu of the antenna. These include an ambient-temperature
‘cold’ load, a ‘hot’ load heated to $\sim 400$ K, an open-ended cable and a
shorted cable. A block diagram showing this arrangement is presented in Fig.
2.

When calibrating the receiver, reflection coefficients are taken of the
calibration source connected to the receiver input ($\Gamma_{\mathrm{cal}}$)
and of the receiver itself ($\Gamma_{\mathrm{rec}}$) as well as power spectral
densities (PSDs) of the input ($P_{\mathrm{cal}}$), the internal reference
load ($P_{\mathrm{L}}$) and the internal reference noise source
($P_{\mathrm{NS}}$) (Monsalve et al., 2017). These measurements are used to
calculate a preliminary ‘uncalibrated’ antenna temperature
$T_{\mathrm{cal}}^{*}$

Figure 2: Diagram of a typical calibration setup. For characterisation of the
receiver, a switch cycles between a calibrator connected to the input and the
two internal references.

$T_{\mathrm{cal}}^{*}=T_{\mathrm{NS}}\left(\frac{P_{\mathrm{cal}}-P_{\mathrm{L}}}{P_{\mathrm{NS}}-P_{\mathrm{L}}}\right)+T_{\mathrm{L}},$
(5)

where $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$ are assumptions for the noise
temperature of the internal reference load and excess noise temperature of the
internal noise source above ambient, respectively. This initial calculation is
used to calibrate out any time-dependent system gain that emerges from a
series of filters, amplifiers and cables, as well as the analogue-to-digital
converter within the experimental apparatus (Monsalve et al., 2017). Each PSD
measurement can be expressed in terms of specific response contributions as
detailed in Bowman et al. (2018)

$\displaystyle P_{\mathrm{cal}}=g_{\mathrm{sys}}\Bigg{[}$ $\displaystyle
T_{\mathrm{cal}}\left(1-|\Gamma_{\mathrm{cal}}|^{2}\right)\left|\frac{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\right|^{2}$
(6) $\displaystyle+$ $\displaystyle
T_{\mathrm{unc}}|\Gamma_{\mathrm{cal}}|^{2}\left|\frac{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\right|^{2}$
$\displaystyle+$ $\displaystyle
T_{\mathrm{cos}}\operatorname{Re}\left(\Gamma_{\mathrm{cal}}\frac{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\right)$
$\displaystyle+$ $\displaystyle
T_{\mathrm{sin}}\operatorname{Im}\left(\Gamma_{\mathrm{cal}}\frac{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\right)+T_{0}\Bigg{]}.$

Here, $g_{\mathrm{sys}}$ is the system gain referenced to the receiver input
and $T_{\mathrm{cal}}$ is our calibrated input temperature.
$T_{\mathrm{unc}}$, $T_{\mathrm{cos}}$, and $T_{\mathrm{sin}}$ are the ‘noise
wave parameters’ introduced by Meys (1978) to calibrate the instrument.
$T_{\mathrm{unc}}$ represents the portion of noise reflected by the antenna
that is uncorrelated with the output noise of the low noise amplifier (LNA).
$T_{\mathrm{cos}}$ and $T_{\mathrm{sin}}$ are the portions of reflected noise
correlated with noise from the LNA (Monsalve et al., 2017; Rogers & Bowman,
2012). In the EDGES experiment, these calibration quantities are modelled
using seven-term polynomials in frequency.

The PSDs for the internal reference load and noise source can similarly be
expressed as in equation 6. However, since the reflection coefficients of the
internal references are typically less than 0.005, they are taken to be zero
in order to simplify the equations

$P_{\mathrm{L}}=g_{\mathrm{sys}}^{*}[T_{\mathrm{L}}\left(1-|\Gamma_{\mathrm{rec}}|^{2}\right)+T_{0}^{*}],$
(7)
$P_{\mathrm{NS}}=g_{\mathrm{sys}}^{*}[\left(T_{\mathrm{L}}+T_{\mathrm{NS}}\right)\left(1-|\Gamma_{\mathrm{rec}}|^{2}\right)+T_{0}^{*}].$
(8)

As shown in Fig. 2, the internal references may be on a separate reference
plane than the receiver input, resulting in a system gain
$g_{\mathrm{sys}}^{*}$ and a noise offset $T_{0}^{*}$ different from those
defined in equation 6. This effect is taken into account by two additional
scale and offset parameters, $C_{1}$ and $C_{2}$, introduced by EDGES
(Monsalve et al., 2017).

Since $C_{1}$ and $C_{2}$ also correct for first-order assumptions in the
noise temperatures of the internal reference load and noise source, we have
chosen to absorb these terms into $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$. This
adjustment allows all calibration parameters, $T_{\mathrm{unc}}$,
$T_{\mathrm{cos}}$, $T_{\mathrm{sin}}$, and an ‘effective’ $T_{\mathrm{NS}}$
and $T_{\mathrm{L}}$, to be solved for in units of kelvin, facilitating a
joint solution of parameters. Expanding equation 5 using equations 6, 7 and 8
yields a linear identity providing a relationship between the uncalibrated
input temperature and a final calibrated temperature of any device connected
to the receiver input

$\displaystyle
T_{\mathrm{NS}}\left(\frac{P_{\mathrm{cal}}-P_{\mathrm{L}}}{P_{\mathrm{NS}}-P_{\mathrm{L}}}\right)+T_{\mathrm{L}}$
$\displaystyle=T_{\mathrm{cal}}\left[\frac{1-|\Gamma_{\mathrm{cal}}|^{2}}{|1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}|^{2}}\right]$
(9)
$\displaystyle+T_{\mathrm{unc}}\left[\frac{|\Gamma_{\mathrm{cal}}|^{2}}{|1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}|^{2}}\right]$
$\displaystyle+T_{\mathrm{cos}}\left[\frac{\operatorname{Re}\left(\frac{\Gamma_{\mathrm{cal}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\right)}{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}\right]$
$\displaystyle+T_{\mathrm{sin}}\left[\frac{\operatorname{Im}\left(\frac{\Gamma_{\mathrm{cal}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\right)}{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}\right],$

where all parameters are frequency-dependent. This is not explicitly shown for
simplicity of notation. For estimation of the noise wave parameters,
$T_{\mathrm{cal}}$, $\Gamma_{\mathrm{cal}}$ and $\Gamma_{\mathrm{rec}}$ are
measured along with the PSDs while $g_{\mathrm{sys}}$ and $T_{\mathrm{0}}$ are
calibrated out. The cold and hot loads exhibit the main temperature references
needed for $T_{\mathrm{L}}$ and $T_{\mathrm{NS}}$. The cables facilitate the
derivation of the noise wave parameters describing spectral ripples from the
noise properties of the receiver by acting as antennas looking at an isotropic
sky with temperatures equal to the cables’ physical temperatures (Rogers &
Bowman, 2012).

###  2.3 Bayesian calibration framework

One possible source of systematics in the calibration methodology used by
EDGES comes from measuring the response of the four external calibrators along
with the receiver reflection coefficient in a laboratory away from where the
instrument is actually deployed (Bowman et al., 2018). This process,
especially with regards to how calibration parameters change, can be non-
trivial. Furthermore, the fixed polynomial order used by EDGES for all noise
wave parameters may underfit or overfit individual parameters and thus ‘fit
out’ data useful for determining systematics or potentially even the 21-cm
signal itself if a joint fit is performed.

In response to these issues, we have developed a calibration pipeline that
improves on the strategies presented in section 2.2. We introduce a novel
Bayesian methodology using conjugate priors for a dynamic application of our
algorithm to be run with data collection regardless of system complexity. Also
included are model selection methods using machine learning techniques for the
optimisation of individual noise wave parameters to combat overfitting and
underfitting, the results of which converge with that of a least-squares
approach when wide priors are adopted. Our pipeline easily incorporates many
more calibrators than the standard four shown in Fig. 2 to increase
constraints on noise wave parameters while identifying possible correlations
between them. A schematic of the improved calibration method is shown in Fig.
3.

Figure 3: Outline of the Bayesian calibration algorithm. Blue blocks represent
data to be taken, red blocks represent calculations and green blocks represent
calculation outputs.

In order to simplify our calibration approach, we first define the following
terms

$X_{\mathrm{unc}}=-\frac{|\Gamma_{\mathrm{cal}}|^{2}}{1-|\Gamma_{\mathrm{cal}}|^{2}},$
(10)
$X_{\mathrm{L}}=\frac{|1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}|^{2}}{1-|\Gamma_{\mathrm{cal}}|^{2}},$
(11)
$X_{\mathrm{cos}}=-\operatorname{Re}\left(\frac{\Gamma_{\mathrm{cal}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\times\frac{X_{\mathrm{L}}}{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}\right),$
(12)
$X_{\mathrm{sin}}=-\operatorname{Im}\left(\frac{\Gamma_{\mathrm{cal}}}{1-\Gamma_{\mathrm{cal}}\Gamma_{\mathrm{rec}}}\times\frac{X_{\mathrm{L}}}{\sqrt{1-|\Gamma_{\mathrm{rec}}|^{2}}}\right),$
(13)
$X_{\mathrm{NS}}=\left(\frac{P_{\mathrm{cal}}-P_{\mathrm{L}}}{P_{\mathrm{NS}}-P_{\mathrm{L}}}\right)X_{\mathrm{L}},$
(14)

which represent initial calibration measurements on $D$ in the frequency
domain for the characterisation of $N_{\mathrm{data}}$ from equation 3 via our
noise wave parameters. It is expected that calibration-related deviations of
$D$ in the time domain are sufficiently curtailed through practical strategies
such as temperature control of the receiver environment. Incorporating these
into equation 9, with some rearrangement, then gives the equation

$X_{\mathrm{unc}}T_{\mathrm{unc}}+X_{\mathrm{cos}}T_{\mathrm{cos}}+X_{\mathrm{sin}}T_{\mathrm{sin}}+X_{\mathrm{NS}}T_{\mathrm{NS}}+X_{\mathrm{L}}T_{\mathrm{L}}=T_{\mathrm{cal}},$
(15)

at each frequency. Here, there are no squared or higher-order terms, allowing
us to take advantage of the linear form by grouping the data and noise wave
parameters into separate matrices

X $\displaystyle\equiv\begin{pmatrix}X_{\mathrm{unc}}\quad
X_{\mathrm{cos}}\quad X_{\mathrm{sin}}\quad X_{\mathrm{NS}}\quad
X_{\mathrm{L}}\end{pmatrix},$ $\displaystyle\boldsymbol{\Theta}$
$\displaystyle\equiv\begin{pmatrix}T_{\mathrm{unc}}\quad T_{\mathrm{cos}}\quad
T_{\mathrm{sin}}\quad T_{\mathrm{NS}}\quad
T_{\mathrm{L}}\end{pmatrix}^{\top}.$ (16)

In these equations, all of our data; the reflection coefficient measurements
and power spectral densities, are grouped in an X vector which forms a matrix
where one of the axes is frequency. The calibration parameters as frequency-
dependent polynomials of varying degree are collected into a
$\boldsymbol{\boldsymbol{\Theta}}$ vector which serves as our model describing
$N_{\mathrm{data}}$. Applying these definitions condenses the calibration
equation into

$\textbf{{T}}_{\mathrm{cal}}=\textbf{{X}}\boldsymbol{\boldsymbol{\Theta}}+\sigma,$
(17)

where $\textbf{{T}}_{\mathrm{cal}}$ is a vector over frequency and $\sigma$ is
a noise vector representing our error. Since EDGES assumes that each power
spectral density measurement is frequency independent, we have assumed that
$\sigma$ is a multivariate normal distribution. This assumption is implicit in
the EDGES analysis in which they use a least-squares minimisation approach for
solving model parameters.

For calibration of the receiver, we are concerned with the construction of
predictive models of the noise wave parameters, $\boldsymbol{\Theta}$, in the
context of some dataset, T. We can use $\boldsymbol{\Theta}$ to calculate the
probability of observing the data given a specific set of noise wave
parameters:

$\displaystyle
p\big{(}\textbf{{T}}\>\big{|}\>\boldsymbol{\Theta},\sigma^{2}\big{)}$
$\displaystyle=$ (18)
$\displaystyle\frac{1}{2\pi\sigma^{2}}^{N/2}\exp{\Bigg{\\{}-\frac{1}{2\sigma^{2}}\left(\textbf{{T}}-\textbf{{X}}\boldsymbol{\Theta}\right)^{\top}\left(\textbf{{T}}-\textbf{{X}}\boldsymbol{\Theta}\right)\Bigg{\\}}},$

where, $N$ is the number of measurements. This distribution on the data is the
likelihood. For the purposes of calibration, T may be
$\textbf{{T}}_{\mathrm{cal}}$ measurements or alternatively,
$\textbf{{T}}_{\mathrm{sky}}$ for prediction of a sky signal. Our model must
also specify a prior distribution, quantifying our initial assumptions on the
values and spread of our noise wave parameters which we specify as a
multivariate normal inverse gamma distribution:

$\displaystyle p\left(\boldsymbol{\Theta},\sigma^{2}\right)\propto$
$\displaystyle\left(\frac{1}{\sigma^{2}}\right)^{a+1+\left(d/2\right)}\times$
(19)
$\displaystyle\exp\left[-\frac{1}{\sigma^{2}}\\{b+\frac{1}{2}\left(\boldsymbol{\Theta}-\boldsymbol{\mu}_{\boldsymbol{\Theta}}\right)^{\top}\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}\left(\boldsymbol{\Theta}-\boldsymbol{\mu}_{\boldsymbol{\Theta}}\right)\\}\right],$

which is proportional up to an integration constant. Here, $a$ and $b$, which
are greater than zero, along with $\textbf{{V}}_{\boldsymbol{\Theta}}$ and
$\boldsymbol{\mu}_{\boldsymbol{\Theta}}$ represent our prior knowledge on the
noise wave parameters. $d$ is the length of our vector $\boldsymbol{\Theta}$.

Equation 18 is determined by a set of values for our model
$\boldsymbol{\Theta}$. We can marginalise out the dependence on
$\boldsymbol{\Theta}$ and our noise term by integrating over the prior
distribution by both $\boldsymbol{\Theta}$ and $\sigma^{2}$ at once. Following
the steps in Banerjee (2009)

$\displaystyle p\left(\textbf{{T}}_{\mathrm{cal}}\right)$ $\displaystyle=\int
p\left(\textbf{{T}}_{\mathrm{cal}}\>\big{|}\>\boldsymbol{\Theta},\sigma^{2}\right)p\left(\boldsymbol{\Theta},\sigma^{2}\right)\mathrm{d}\boldsymbol{\Theta}\mathrm{d}\sigma^{2}$
(20)
$\displaystyle=\frac{b^{a}\Gamma\left(a^{*}\right)\sqrt{|\textbf{{V}}^{*}|}}{{b^{*}}^{a^{*}}\Gamma\left(a\right)\sqrt{|\textbf{{V}}_{\boldsymbol{\Theta}}|}}(2\pi)^{-N/2},$

where

$\displaystyle a^{*}$ $\displaystyle=a+\frac{N}{2},$ (21) $\displaystyle
b^{*}$
$\displaystyle=b+\frac{1}{2}[\boldsymbol{\mu}_{\boldsymbol{\Theta}}^{\top}\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}\boldsymbol{\mu}_{\boldsymbol{\Theta}}+\textbf{{T}}_{\mathrm{cal}}^{\top}\textbf{{T}}_{\mathrm{cal}}-\boldsymbol{\mu}^{*\top}\textbf{{V}}^{*-1}\boldsymbol{\mu}^{*}],$
$\displaystyle\boldsymbol{\mu}^{*}$
$\displaystyle=\left(\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}+\textbf{{X}}^{\top}\textbf{{X}}\right)^{-1}\left(\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}\boldsymbol{\mu}_{\boldsymbol{\Theta}}+\textbf{{X}}^{\top}\textbf{{T}}_{\mathrm{cal}}\right),$
$\displaystyle\textbf{{V}}^{*}$
$\displaystyle=\left(\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}+\textbf{{X}}^{\top}\textbf{{X}}\right)^{-1},$

and $\Gamma\left(x\right)$ represents the Gamma function, not to be confused
with the notation for our reflection coefficients. Equation 20 is the
evidence, which gives the probability of observing the data
$\textbf{{T}}_{\mathrm{cal}}$ given our model.111It is in fact better to use
the equivalent more numerically stable expression
$b^{*}=b+\boldsymbol{q}^{\top}\boldsymbol{q}+\boldsymbol{q}^{\top}\textbf{{X}}\textbf{{V}}_{\boldsymbol{\Theta}}\textbf{{X}}^{\top}\boldsymbol{q}$,
where
$\boldsymbol{q}=\textbf{{T}}_{\mathrm{cal}}-\textbf{{X}}\boldsymbol{\mu}^{*}$
to avoid cancellation of large terms.

With the prior distribution specified, we use Bayes’ equation to invert the
conditioning of the likelihood and find the posterior using the likelihood,
prior and evidence:

$p\left(\boldsymbol{\Theta},\sigma^{2}\>\big{|}\>\textbf{{T}}_{\mathrm{cal}}\right)=\frac{p\left(\textbf{{T}}_{\mathrm{cal}}\>\big{|}\>\boldsymbol{\Theta},\sigma^{2}\right)p\left(\boldsymbol{\Theta},\sigma^{2}\right)}{p\left(\textbf{{T}}_{\mathrm{cal}}\right)}.$
(22)

Similarly from Banerjee (2009), this can be written as

$\displaystyle p\Bigl{(}\boldsymbol{\Theta},\sigma^{2}\>\big{|}$
$\displaystyle\textbf{{T}}_{\mathrm{cal}}\Bigl{)}\propto\left(\frac{1}{\sigma^{2}}\right)^{a^{*}+\frac{d}{2}+1}\times$
(23)
$\displaystyle\exp{\Bigg{\\{}-\frac{1}{\sigma^{2}}\Bigg{[}b^{*}+\frac{1}{2}\left(\boldsymbol{\Theta}-\boldsymbol{\mu}^{*}\right)^{\top}\textbf{{V}}^{*-1}\left(\boldsymbol{\Theta}-\boldsymbol{\mu}^{*}\right)\Bigg{]}\Bigg{\\}}}.$

The posterior distribution represents the uncertainty of our parameters after
analysis, reflecting the increase in information (Nagel, 2017). We highlight
the difference between the ‘likelihood-only’ least-squares approach versus the
Bayesian approach with the former being a special case of the latter with very
wide priors demonstrable when
$\textbf{{V}}_{\boldsymbol{\Theta}}\rightarrow\infty\Rightarrow\textbf{{V}}_{\boldsymbol{\Theta}}^{-1}\rightarrow
0$, and $\boldsymbol{\mu}^{*}$ becomes $\boldsymbol{\Theta}$. The transition
from ‘non-starred’ variables to ‘starred’ variables represents our ‘Bayesian
update’ of the prior to the posterior noise wave parameters in light of the
calibration data $\textbf{{T}}_{\mathrm{cal}}$.

As we can see, the posterior distribution is in the same probability
distribution family as equation 19, making our prior a conjugate prior on the
likelihood distribution. The use of conjugate priors gives a closed-form
solution for the posterior distribution through updates of the prior
hyperparameters via the likelihood function (Banerjee, 2009; Orloff & Bloom,
2013). The resulting numerical computation is many orders of magnitude faster
than MCMC methods relying on full numerical sampling and permits an in-place
calculation in the same environment as the data acquisition. This becomes
particularly useful for the speed of the algorithm as frequency dependence is
introduced in which the computations would not be manageable without conjugate
gradients.

To allow for a smooth frequency dependency, we promote each of our noise wave
parameters in section 2.3 to a vector of polynomial coefficients

$T_{\mathrm{i}}=\begin{pmatrix}T_{\mathrm{i}}^{[0]},&T_{\mathrm{i}}^{[1]},&T_{\mathrm{i}}^{[2]},&...,&T_{\mathrm{i}}^{[n]}\end{pmatrix},$
(24)

where $i$ is our noise wave parameter label; $i\in\\{\mathrm{unc,\ cos,\ sin,\
NS,\ L}\\}$, modelled using $n+1$ polynomial coefficients. Likewise

$\textbf{{X}}_{i}=\begin{pmatrix}\textbf{{X}}_{i},&\textbf{{X}}_{i}\left(\frac{\nu}{\nu_{0}}\right),&\textbf{{X}}_{i}{\left(\frac{\nu}{\nu_{0}}\right)}^{2},&...,&\textbf{{X}}_{i}{\left(\frac{\nu}{\nu_{0}}\right)}^{n}\end{pmatrix},$
(25)

where $\nu$ is a vector of input frequencies which are raised to powers up to
$n$. For a vector of $n$’s attributed to our calibration parameters, under
this notation multiplication in equation 17 is element-wise and equation 20 is
effectively $p\left(\textbf{{T}}_{\mathrm{cal}}|\textbf{{n}}\right)$. Assuming
a uniform prior on n, inverting Bayes’ theorem gives
$p\left(\textbf{{n}}|\textbf{{T}}_{\mathrm{cal}}\right)$ for use in model
comparison in which the relative probabilities of models can be evaluated in
light of the data and priors. Occam’s razor advises whether the extra
complexity of a model is needed to describe the data (Trotta, 2008),
permitting optimisation of the polynomial orders for individual noise wave
parameters as detailed in section 3.3. By taking a random sampling of the
resulting posterior, we characterise the noise wave parameters as multivariate
distributions depicted in contour plots which exhibit a peak value accompanied
by $1\sigma$ and $2\sigma$ variance as well as correlation between parameters
inferred from a covariance matrix.

Following characterisation of the receiver, we next apply the
$\textbf{{T}}_{\mathrm{cal}}$ from our calibration to a set of raw antenna
data $\hat{\textbf{{X}}}$ for prediction of our sky signal,
$\textbf{{T}}_{\mathrm{sky}}$, from equation 3. The predictions for the data
follow from the _posterior predictive distribution_

$p\left(\textbf{{T}}_{\mathrm{sky}}\>\big{|}\>\textbf{{T}}_{\mathrm{cal}}\right)=\int
p\left(\textbf{{T}}_{\mathrm{sky}}\>\big{|}\>\boldsymbol{\Theta},\sigma^{2}\right)p\left(\boldsymbol{\Theta},\sigma^{2}\>\big{|}\>\textbf{{T}}_{\mathrm{cal}}\right)\mathrm{d}\boldsymbol{\Theta}\mathrm{d}\sigma^{2}.$
(26)

The first probability in the integral is the likelihood for our antenna
measurement $\textbf{{T}}_{\mathrm{sky}}$ and the second is our posterior from
equation 23. Following the steps in Banerjee (2009), this can be shown to be a
multivariate Student’s t-distribution written as:

$\displaystyle p\Big{(}$
$\displaystyle\textbf{{T}}_{\mathrm{sky}}\>\big{|}\>\textbf{{T}}_{\mathrm{cal}}\Big{)}=\frac{\Gamma\left(a^{*}+\frac{d}{2}\right)}{\Gamma\left(a^{*}\right)\pi^{\frac{d}{2}}|2b^{*}\left(I+\hat{\textbf{{X}}}\textbf{{V}}^{*}\hat{\textbf{{X}}}^{\top}\right)|^{\frac{1}{2}}}$
(27)
$\displaystyle\times\left[1+\frac{\left(\textbf{{T}}_{\mathrm{sky}}-\hat{\textbf{{X}}}\boldsymbol{\mu}^{*}\right)^{\top}\left(I+\hat{\textbf{{X}}}\textbf{{V}}^{*}\hat{\textbf{{X}}}^{\top}\right)^{-1}\left(\textbf{{T}}_{\mathrm{sky}}-\hat{\textbf{{X}}}\boldsymbol{\mu}^{*}\right)}{2b^{*}}\right]^{-\left(a^{*}+\frac{d}{2}\right)},$

where $I$ is the $N\times N$ identity matrix and $a^{*}$, $b^{*}$,
$\boldsymbol{\mu}^{*}$ and $\textbf{{V}}^{*}$ are defined in equation 21. This
new distribution on $\textbf{{T}}_{\mathrm{sky}}$ corresponds to a set of
points with error bars and represents the calibrated sky temperature as the
output of the receiver.

##  3 Empirical modelling and simulations

To verify the performance of our pipeline and highlight features of the
algorithm, we evaluate the results of self-consistency checks using empirical
models of data based on measurements taken in the laboratory. To make this
data as realistic as possible, we used actual measurements of the reflection
coefficients of many types of calibrators (see table 1) to generated power
spectral densities using equations 6, 7 and 8 given a set of realistic model
noise wave parameters along with some assumptions about the noise, which are
described in section 3.4. The impedance of the calibrators which were measured
with a vector network analyser (VNA) and used in our pipeline are shown on a
Smith chart in Fig. 4

Table 1: Table of calibrators used in the creation of our empirical data models for analysis. Calibrators are added in pairs in the order below when increasing the number of calibration sources used by our algorithm. Calibrator | Temperature
---|---
Cold load ($50\ \Omega$) | 298 K
Hot load ($50\ \Omega$) | 373 K
Gore cable $+5\ \Omega$ | 298 K
Gore cable $+500\ \Omega$ | 298 K
Gore cable $+31\ \Omega$ | 298 K
Gore cable $+81\ \Omega$ | 298 K
$25\ \Omega$ resistor | 298 K
$100\ \Omega$ resistor | 298 K
Figure 4: Smith chart (Argand diagram) showing the measured complex impedance
of the calibrators used in the Bayesian pipeline across a range of
frequencies.

We start by demonstrating the importance of correlation between noise wave
parameters when determining their values to provide a better calibration
solution for the reduction of systematic features in the data such as
reflections (section 3.1). We then show the increased constraints on these
noise wave parameters attributed to the inclusion of more calibrators than the
standard number of four (section 3.2). Following this, we illustrate the
effectiveness of model selection for the optimisation of individual noise wave
parameters to prevent the loss of information resulting from overfitting or
underfitting of the data (section 3.3). Finally, these features are
incorporated into a calibration solution applied to a $50\ \Omega$ load
(section 3.4).

###  3.1 Correlation between noise wave parameters

In this section, we show the first major feature of our Bayesian pipeline; the
consideration of correlation between noise wave parameters when deriving their
values. This is best demonstrated when noise is introduced in an idealised way
as to retain a form matching the Gaussian form of our mathematical model. To
do this, empirical models of power spectral densities are calculated from
equations 6, 7 and 8 using measurements of $\Gamma_{\mathrm{rec}}$,
$\Gamma_{\mathrm{cal}}$ and $T_{\mathrm{cal}}$ for the cold and hot loads, as
well as a set of realistic noise wave parameters. Gaussian noise of one unit
variation is then added to the $T_{\mathrm{cal}}$ measurements after the
calculation to conserve its Gaussian form. This data is submitted to our
algorithm and the resulting posterior distributions for coefficients of the
polynomial noise wave parameters are compared to the initial values.

Such posterior distributions can be seen in Fig. 5 showing the results of
models using only the cold load (grey posterior), only the hot load (red
posterior) and using both loads in tandem (blue posterior). For these
calculations we chose a set of model noise wave parameters as constants across
the frequency band;

$\displaystyle T_{\mathrm{unc}}=250\ \mathrm{K}$ $\displaystyle
T_{\mathrm{cos}}=190\ \mathrm{K}$ $\displaystyle T_{\mathrm{sin}}=90\
\mathrm{K}$ $\displaystyle T_{\mathrm{NS}}=1200\ \mathrm{K}$ $\displaystyle
T_{\mathrm{L}}=298\ \mathrm{K}$

In Fig. 5, a strong correlation between the $T_{\mathrm{L}}$ and
$T_{\mathrm{NS}}$ is evident as the hot-load posterior is highly skewed as
expected from equations 11 and 14. The resulting intersection of posteriors
from the individual loads facilitate the derivation of noise wave parameters
as the dual-load posterior is found within the region of posterior overlap
crossing with the values of the model shown in the inset of Fig. 5. Retrieval
of the noise wave parameter values using correlations between them found in
the data demonstrate the relevance of this information which is not taken into
account in previous calibration techniques.

Figure 5: Plot showing the joint posteriors of $T_{\mathrm{L}}$ and
$T_{\mathrm{NS}}$ for models using the cold load, the hot load, and both loads
concurrently shown as the grey, red and blue posteriors respectively. The
black cross hairs mark the noise wave parameter values used to generate data
submitted to the pipeline. A zoom-in of the posterior intersection is provided
to illustrate the constraint of noise wave parameter values attributed to the
correlation between parameters.

###  3.2 Constraints with additional calibrators

A nice feature of our pipeline is the ability to include as many calibrators
as required to constrain the calibration parameters. For analysis, six more
calibrators are introduced in pairs following the order presented in table 1.
We include data generated from measurements of multiple resistors terminating
a high quality 25 m cable made by GORE®. Data for these calibrators is once
again generated using fixed terms and Gaussian noise of one unit variation
added to $T_{\mathrm{cal}}$ as discussed above. Fig. 6 shows the results of
models using four, six, and eight calibrators.

Figure 6: Posterior results of our pipeline using data from four, six and
eight calibrators shown in grey, red and blue, respectively. Cross hairs mark
the values of noise wave parameters used to generate the data. These values
fall within $1\sigma$ of the posterior mean values. We can see that the
constraint on noise wave parameter values increases with the number of
calibrators used in our pipeline which is encouraging.

As shown, the inclusion of more calibrators increases the constraint on the
resulting noise wave parameters. However, we note that after the inclusion of
four calibrators, the relative additional constraint decreases with each
additional calibrator and thus the use of more than eight calibrators would be
unnecessary. The values of noise wave parameters used to generate the data as
indicated by the cross hairs in Fig. 6 all fall within $1\sigma$ of our
pipeline’s resulting posterior averages for models using all eight
calibrators.

###  3.3 Optimisation of individual noise wave parameters

The final highlight of our Bayesian pipeline is a the use of machine learning
techniques to optimise individual noise wave parameters. This is advantageous
as a blanket set of order-seven polynomials applied to all noise wave
parameters, such as done in the EDGES experiment, may underfit or overfit
individual parameters and misidentify systematics or information about the
signal being measured.

The optimisation procedure compares the evidences (equation 20) of different
models to determine the vector of noise wave parameter polynomial coefficients
n that best describes the data as briefly mentioned at the end of section 2.3.
Since the model favoured by the data will have the highest evidence, we use a
steepest descent procedure to compare models in ‘n-space’ and determine the
direction of the gradient in ‘evidence-space’. After multiple iterations, this
brings us to the model with the maximal evidence. Since n consists of five
numbers corresponding to the number of polynomial coefficients for each of the
five noise wave parameters, models are generated by individually increasing
each index of n by 1. We expect the evidence to follow an ‘Occam’s cliff,’ in
which the evidence sharply increases preceding the optimal n with a slow fall
off following the maximum.

To demonstrate this, data is generated using measurements from all eight
calibrators of table 1 and noise wave parameters as second-order polynomials

$\displaystyle T_{\mathrm{unc}}=x^{2}-3x+250\ \mathrm{K}$ $\displaystyle
T_{\mathrm{cos}}=2x^{2}+190\ \mathrm{K}$ $\displaystyle
T_{\mathrm{sin}}=3x^{2}+8x+90\ \mathrm{K}$ $\displaystyle
T_{\mathrm{NS}}=4x^{2}+5x+1200\ \mathrm{K}$ $\displaystyle
T_{\mathrm{L}}=5x^{2}+10x+298\ \mathrm{K}$

where $x$ is our normalised frequency. Gaussian noise of one unit variation is
applied to the calibrator input temperatures as before. The evidences of
various models are plotted in Fig. 7 in which an Occam’s cliff can be seen
peaking at polynomial order two. As expected from the plot, the steepest
descent algorithm finds that noise wave parameters modelled as second-order
polynomials best describe the data.

Figure 7: Evidence of multiple models are plotted which display the Occam’s
cliff. Data is generated using noise wave parameters as order-2 polynomials.
We see that for the model with the highest evidence, that is, the model
favoured by the data, the number of polynomial coefficients matches that of
the model noise wave parameters.

###  3.4 Application with realistic noise

To demonstrate the robustness of our pipeline, we conducted self-consistency
checks using empirically modelled data with a more complicated noise model.
This data was generated using reflection coefficients of eight calibrators and
the receiver measured in the laboratory. These reflection coefficients were
then smoothed using a cubic smoothing spline (Weinert, 2009) in order to
maintain their approximate shape over frequency. The same second-order noise
wave parameters detailed in section 3.3 are used with the reflection
coefficients to generate our model power spectral densities. Following this,
we added of order 1% Gaussian noise independently to the smoothed
$\Gamma_{\mathrm{rec}}$ and $\Gamma_{\mathrm{cal}}$ as well as
$P_{\mathrm{cal}}$ to more accurately represent the instrument noise from
measurement equipment such as vector network analysers. No noise was added to
the calibrator input temperatures. This results in a model that does not match
the Gaussian form of our mathematical model as in the previous sections and
thus does not demonstrate the features of our pipeline as explicitly, but is
more representative of data set expected from measurements in the field. Data
for the receiver and the cold load generated using this noise model are shown
in Fig. 8.

Figure 8: Power spectral densities and reflection coefficients for the
receiver and the cold load generated under our realistic noise model.

Using data generated for all eight calibrators with our realistic noise model,
the calibration algorithm selects optimal polynomial orders matching those of
the model noise wave parameters whose values fall within within $1\sigma$ of
the posterior peak values as shown in Fig. 9. For these higher order tests, we
use fgivenx plots which condense noise wave parameter posteriors into samples
that can be compared to the model parameter values instead of comparing each
individual coefficient (Handley, 2018).

When this calibration model is used to calibrate an ambient-temperature $50\
\Omega$ load, the RMS error between the calibrated temperature and the
measured temperature is 8 mK, well within the $1\sigma$ noise level (bottom
right panel of Fig. 9). This level of accuracy is comparable to the 26 mK
noise floor estimated for the EDGES pipeline in 2016 (Monsalve et al., 2017).

By individually adjusting each component of noise arising in our realistic
noise model, we may determine what kind of noise our calibration algorithm is
most sensitive to, as well as calculate the maximum amount of noise
permissible for a specified level of systematic feature reduction. These
topics are intended to be explored in a future work.

Figure 9: Results from 1000 samples using data generated with our more
realistic noise model (shown in black). The second-order noise wave parameters
shown in red are used to generate the data inputted to our pipeline. The
polynomial order and values of the noise wave parameters that best suit the
data according to our algorithm match that of the empirical model. This
solution is applied to an ambient-temperature load, shown in the bottom right
panel as our predictive $\hat{y}$ from equation 27, and calibrates it to
within $1\sigma$ of ambient temperature.

##  4 Conclusions

Here we presented the development of a calibration methodology based on the
procedure used by EDGES but with key improvements to characterise reflections
arising at connections within the receiver. Our pipeline utilises the Dicke
switching technique and a Bayesian framework in order to individually optimise
calibration parameters while identifying correlations between them using a
dynamic algorithm to be applied in the same environment as the data
acquisition. In a comprehensive investigation, we have evaluated our
algorithm’s interpretation of empirical models of data which have been
generated from known noise wave parameters and a realistic noise model. The
solution, applied to an ambient-temperature $50\ \Omega$ load, produces a
calibrated temperature with an RMS residual temperature of 8 mK. Future work
for the pipeline regards application of real calibrator data, optimisation of
noise wave parameter coefficients through marginalisation techniques and
incorporation into an end-to-end simulation based on an entire experimental
apparatus to better understand error tolerances. The flexibility of the
algorithm attributed to our novel approach allows its application to any
experiment relying on similar forms of calibration such as REACH (de Lera
Acedo et al., 2021), were we intend to use this for in-the-field and on-the-
fly radiometer calibration.

## Acknowledgements

ILVR would like to thank S. M. Masur for her helpful comments. WJH was
supported by a Gonville & Caius Research Fellowship, STFC grant number
ST/T001054/1 and a Royal Society University Research Fellowship. We would like
to thank the The Cambridge-Africa ALBORADA Research Fund for their support. We
would also like to thank the Kavli Foundation for their support of the REACH
experiment.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * de Lera Acedo et al. (2021) de Lera Acedo E., et al., 2021, REACH: Radio Experiment for the Analysis of Cosmic Hydrogen, in press
  * Anstey et al. (2020) Anstey D., de Lera Acedo E., Handley W., 2020, arXiv e-prints, p. arXiv:2010.09644
  * Banerjee (2009) Banerjee S., 2009, Bayesian Linear Model: Gory Details, Online, http://www.biostat.umn.edu/~ph7440/pubh7440/BayesianLinearModelGoryDetails.pdf
  * Barkana (2018) Barkana R., 2018, Nature, 555, 71
  * Bilous, A. V. et al. (2016) Bilous, A. V. et al., 2016, A&A, 591, A134
  * Bowman et al. (2018) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67
  * Cohen et al. (2017) Cohen A., Fialkov A., Barkana R., Lotem M., 2017, Monthly Notices of the Royal Astronomical Society, 472, 1915
  * Dicke (1946) Dicke R. H., 1946, Review of Scientific Instruments, 17, 268
  * Ewall-Wice et al. (2020) Ewall-Wice A., Chang T.-C., Lazio T. J. W., 2020, MNRAS, 492, 6086
  * Field (1958) Field G. B., 1958, The Astrophysical Journal, 129, 536
  * Furlanetto et al. (2009) Furlanetto S. R., et al., 2009, in astro2010: The Astronomy and Astrophysics Decadal Survey. p. 83 (arXiv:0902.3011)
  * Handley (2018) Handley W., 2018, The Journal of Open Source Software, 3, 849
  * Hills et al. (2018) Hills R., Kulkarni G., Meerburg P. D., Puchwein E., 2018, Nature, 564, E32
  * Meys (1978) Meys R. P., 1978, IEEE Transactions on Microwave Theory and Techniques, MTT-26, 34
  * Monsalve et al. (2017) Monsalve R. A., Rogers A. E. E., Bowman J. D., Mozdzen T. J., 2017, The Astrophysical Journal, 835, 49
  * Nagel (2017) Nagel J. B., 2017, PhD thesis, ETH Zurich
  * Orloff & Bloom (2013) Orloff J., Bloom J., 2013, Conjugate priors: Beta and normal, Online, https://math.mit.edu/~dav/05.dir/class15-prep.pdf
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 8
  * Razavi-Ghods (2018) Razavi-Ghods N., 2018, A new instrument to measure the global EoR signal (Internal Engineering Memo), University of Cambridge, doi:10.5281/zenodo.4761931
  * Rogers & Bowman (2012) Rogers A. E. E., Bowman J. D., 2012, Radio Science, 47
  * Trotta (2008) Trotta R., 2008, Contemporary Physics, 49, 71
  * Weinert (2009) Weinert H. L., 2009, Computational Statistics and Data Analysis, 53, 932
  * Zaldarriaga et al. (2004) Zaldarriaga M., Furlanetto S. R., Hernquist L., 2004, ApJ, 608, 622