# Overcoming 1 part in $10^{9}$ of Earth angular rotation rate measurement
with the G Wettzell data

A. D. V. Di Virgilio and G. Terreni
INFN-Pisa largo B. Pontecorvo 3, 56127 Pisa, Italy
<EMAIL_ADDRESS>
A. Basti, N. Beverini, G. Carelli, D. Ciampini, F. Fuso, E. Maccioni and P.
Marsili
Dep.of Physics Univ. of Pisa, largo B. Pontecorvo 3, 56127 Pisa, Italy
J. Kodet and K. U. Schreiber
Technical University of Munich,
FESG, Geodetic Observatory Wettzell,
93444 Bad Kotzting, Germany

###### Abstract

The absolute measurement of the Earth angular rotation rate with ground-based
instruments becomes challenging if the 1 part in $10^{9}$ of precision has to
be obtained. This threshold is important for fundamental physics and for
geodesy, to investigate effects of General Relativity and Lorentz violation in
the gravity sector and to provide the fast variation of the Earth rotation
rate. High sensitivity Ring Laser Gyroscopes (RLG) are currently the only
promising technique to achieve this task in the near future, but their
precision has been so far limited by systematics related to the laser
operation. In this paper we analyze two different sets of observations, each
of them three days long. They were obtained from the G ring laser at the
Geodetic Observatory Wettzell. The applied method has been developed for the
GINGERINO ring laser in order to identify and extract the laser systematics.
For the available data sets the residuals show mostly white noise behavior and
the Allan deviation drops below 1 part in $10^{9}$ after an integration time
of about $10^{4}$ s.

##  1 Introduction

At present large scale ring laser gyroscopes (RLGs) are the most sensitive
instruments to measure absolute angular rotation rates[1, 2, 3, 4, 5, 6]. They
are based on high finesse square optical cavities where an active medium is
present, and two counterpropagating laser beams are generated. The frequencies
of the two beams depend on the effective time the photon takes to travel along
the perimeter, which is different when the gyro is rotated. However, non-
reciprocal effects in the cavity from the laser excitation process will cause
a systematic bias and this has to be avoided.

Furthermore, there can be small deviations from the expected rate of rotation,
due to fundamental laws of physics. One example is the precession of the frame
of reference, as predicted by General Relativity (GR), related to
gravitoelectric and gravitomagnetic effects.

For a RLG rigidly connected with the Earth surface, the observed Sagnac
frequency is the difference in angular frequency of the two laser beams
$\omega_{s}$:

$\omega_{s}=8\pi\frac{A}{P\lambda}\Omega_{\oplus}\cos\theta\;,$ (1)

where $A$ is the area and $P$ the perimeter of the ring cavity, $\lambda$ the
optical wavelength, $\Omega_{\oplus}$ the angular rotation rate, equivalent to
the Earth rotation rate. $\theta$ is the angle between the optical cavity area
vector and the $\Omega_{\oplus}$ rotation axis.

G has demonstrated a sensitivity of 12 prad/s in 1 s integration time and is
able to operate continuously and unattended for months.[8] The sensitivity is
a function of the size of the ring cavity and cavities of 77 and 121 m
perimeter,[9] have been explored, several years ago. More recently the four
component ring laser array ROMY has been installed in the geophysical
observatory of Bavaria, Germany, which comprises 4 RLGs with triangular
cavities, each with a perimeter of 36 m [2, 10]. Typically, the most sensitive
existing devices capable of long term continuous operation employ square
optical cavities with several meters on a side, rigidly attached to the
ground. The mechanical structure of the optical cavity plays a big role in the
performance of RLG and monolithic optical cavity structures have been the
first RLGs to obtain relevant performance, getting close to tens prad/s level
of sensitivity [1]. The G RLG, in operation at the Geodetic Observatory
Wettzell, is the best known example and the optimal choice to have a very
stable laser cavity. However, it is very difficult to build and cannot be
reorientated easily. Over more than 20 years, heterolythic (HL) cavities have
been developed, basically consisting of a rigid mounting frame giving hold to
different mechanical elements, which support the mirrors contained inside
metallic vacuum chambers.

GINGERINO has been the first HL RLG operating continuously with high
sensitivity [11]. GINGERINO takes advantage of the quiet location inside the
Gran Sasso, Italy, underground laboratory. Despite that, the standard
deviation of the reconstructed Sagnac frequency of GINGERINO is typically more
than 50 times larger than that of G [3, 4], mostly due to small mechanical
issues, since the composed structure is not yet rigid enough[12]. Work is
currently in progress to improve the mechanical HL design in order to increase
the sensor stability. Despite these mechanical shortcomings in GINGERINO, the
previous investigations [3, 4] indicate that the intrinic noise level is lower
than the expected shot noise limit for a RLG[13]. This model does not take the
presence of the active medium inside the cavity into account and there is work
in progress to improve the theoretical model by the Italian group of GINGER.

The main signal obtained from a RLG is the interferogram of the two
counterpropagating beams transmitted at one corner of the square cavity. The
Sagnac frequency $\omega_{s}$, which is caused by the rate of rotation, must
be not confused with the actually observed beat note $\omega_{m}$ of the
interferogram [14, 15, 16, 17], since the latter is biased by laser
systematics, usually caused by backscatter coupling and a null shift
contribution. An original analysis scheme to remove laser systematics has been
developed and successfully applied to analyse the data of GINGERINO [18, 19,
3, 4]. This effort shows how any change due to the laser dynamics can be
approximated with the sum of several terms, analytically evaluated using the
laser parameters of the Lamb theory. These in turn can be determined from the
signals provided by the RLG. It is important to remark that these laser
systematics are not a stochastic effect, since they are caused by all sorts of
small changes affecting the active optical cavity. These effects are much
smaller in monolithic RLGs.

##  2 Analysis procedure

The laser dynamics can be described by a set of differential equations
containing several parameters, which are known as Lamb parameters [14, 15, 18,
19]. They can be calculated from the diagnostic signals, available from the
ring laser intensities, in particular those providing the AC ($IS_{1,2}$) and
the DC ($PH_{1,2}$) levels of the detectors, each looking at a different laser
beam (hereafter denoted as monobeam signals) and their relative phase offset
$\epsilon$. The model refers to the intracavity power, while the signals are
taken outside the cavity; it is therefore necessary to know the values of
mirror transmission and detector gain.

The analysis proceeds in steps. The first approximation of the Sagnac
frequency, denoted as $\omega_{s0}$, is analytically evaluated as follows:

$\omega_{s0}=\frac{1}{2}\sqrt{\frac{2I_{S1}I_{S2}\omega_{m}^{2}\cos(2\epsilon)}{PH_{1}PH_{2}}+\omega_{m}^{2}}+\frac{\omega_{m}}{2}\;,$
(2)

where $\omega_{m}$ is the experimental interferometric angular frequency, that
is the beat note of the interference of the two counterpropagating beams
coming out from the cavity. It must be noted that the monobeam intensities
enter Eq. 2 as the ratio of their AC and DC components, leading to mirror
transmission and the differences in electronic gains to play a smaller role.
For many applications, $\omega_{s0}$ is indeed a good approximation of
$\omega_{s}$.

Nevertheless, the different terms of the equation are certainly affected by
errors, as for example the dark current of the applied photodiodes or any non-
linearity in the electronic gain. In order to improve the $\omega_{s0}$
evaluation, six terms, namely $\omega_{\xi}$, have been elaborated, assuming
small errors in the measurement data used in Eq. 2 and expanding the equation
to the first order. The first term was already discussed in the first paper
dedicated to the analysis method [18], the other ones have been added more
recently. We remark that the six $\omega_{\xi}$ terms do not depend on the
laser dynamics and are meant to improve $\omega_{s0}$ by correcting errors
associated with the measurements themselves.

The null shift depends on the Lamb parameters associated with the laser
excitation and is obtained by the term $\omega_{ns1}$, the first order
expansion of the theory. Higher order expansion terms can be calculated, but
in the following analysis only the first one will be used. The null shift is
strictly connected to the non-reciprocity of the optical path in the two
directions, related in turn to dissipative processes. In terms of the Lamb
parameters of the laser functions, $\omega_{ns1}$ is connected to the
difference $\mu_{c}-\mu_{cc}$ between the cavity losses, where $\mu_{c,cc}$
represent losses in clockwise and counter-clockwise propagation directions.
However, it has been demonstrated in [18] that, assuming a quasi-stationary
laser and the same value for the parameter $\beta$ for the two
beams111$\beta_{1,2}$ is the self-saturation parameter of the laser transition
for the two counter-propagating beams. In the model it is assumed
$\beta_{c}=\beta_{cc}=\beta$., the loss difference is made explicit, so that
only one of the two $\mu_{c,cc}$ is a free parameter. In the following, we
will consider only the clockwise cavity loss, which will be indicated simply
as $\mu$; $\omega_{ns1}$ and $\beta$ are proportional to $\mu$ and completely
defined by the theory in terms of the available signals. The null shift
correction has values in the region of mHz (several ppm), so the accuracy of
$\mu$, usually measured by the ring-down time of the optical cavity with
percent accuracy, can severely limit the $\omega_{s}$ reconstruction. The
analysis developed for GINGERINO assumes constant the plasma temperature and
pressure, while $\mu$ changes with time, accordingly any change is interpreted
as change of $\mu$, but it could be due to the other parameters, or to the
electronic circuit regulating the gain tube. In the future the analysis model
will be further expanded in order to better identify the origin of the
changes. A suitable procedure is developed to take into account variations of
$\mu$ with time and refine the identification of the null shift. Changes of
$\mu$ in time, however, are small and slow enough to not invalidate the
assumption of a quasi-stationary regime. It is possible to consider
$\mu(t)=\mu(t_{0})+\delta\mu(t)$, $t_{0}$ being the origin of the expansion in
the series. It is convenient to describe $\delta\mu(t)$ using the available
signals. To this aim, the laser gain monitor signal can be used, but in more
recent versions of the analysis the Lamb parameter $\beta$ has been
considered, since it is rather constant in time and proportional to $\mu$. By
using the relationship reported in the appendix of [19] and the parameters of
the RLG G, we have

$\beta=\frac{2.98772\mu}{2.98772-PH_{1}}\;,$ (3)

where $PH_{1}$ is the DC value of monobeam 1, in Volts. Equation 3 shows that
$\mu$ is a proportionality constant, accordingly the quantity that can be
evaluated does not contain $\mu$. Assuming $\mu(t)=\mu(t_{0})+\delta\mu(t)$,
and taking into account that $\beta$ is rather constant, since it is related
to the laser transition, $\delta\mu(t)$ can be evaluated at the first
approximation as follows:

$\displaystyle\bar{\beta}$ $\displaystyle=$
$\displaystyle\frac{\beta(t)}{\mu}$ $\displaystyle\delta\mu(t)$
$\displaystyle=$
$\displaystyle\mu(t_{0})\times\left(\frac{\bar{\beta}(t_{0})}{\bar{\beta}(t)}-1\right)\;.$
(4)

In this way, $\mu(t_{0})$, the loss at the time $t_{0}$, has to be determined
by statistical means, with $t_{0}$ arbitrarily chosen: in the present analysis
$t_{0}$ is the central point of the data set. The Lamb parameters can be
considered constant over short time intervals; therefore, the stationary
mathematical relationships can be considered valid for short time intervals.
Accordingly the required analytical relationships are elaborated at high
frequency rate, and decimated afterwards down to the desired low frequency
rate, details can be found in the related literature[18, 19].

The second step of the analysis determines $\omega_{s}$ with a linear
regression to optimize the subtraction of the $\omega_{\xi}$ terms,
$\omega_{ns1}$ and $\omega_{ns1}\times\delta\mu$. In this way, $\omega_{s}$ is
recovered, but it is still necessary to further identify and subtract other
known signals. The G RLG is well isolated, but in any case affected by global
and local motion of the Earth crust, such as the diurnal polar motion and the
solid Earth tides. In the present analysis the direct effect of the
environmental disturbances on the mechanical apparatus is considered
negligible, owing to the monolithic and stable ZERODUR222ZERODUR is a glass
with very low thermal expansion coefficient, that can be considered
practically negligible at room temperature. structure. It is assumed that all
external effects, hereafter denoted geodetic components, can be described by
sufficiently accurate models. The geodetic components are well known, but they
are added to the linear regression along with the other terms in order to
avoid biases induced by previous analyses. The terms relevant to describe the
main geodetic components are taken from model data provided by the G group.
Furthermore, it is also possible to take the angular rotation around the
vertical in the local reference frame at the G latitude and longitude into
account. They are based on the publications of the International Earth
Rotation and Reference Systems Service (IERS)333https://hpiers.obspm.fr/eop-
pc/index.php?index=C04&lang=en. We denote relevant terms as $FGEO$. It is
important to remark that $FGEO$ contains the Chandler and the Annual wobble,
whose values are published with several days delay. For this reason these
signals are not included in the near realtime data files of G. Finally, we
have verified that the calibration and alignment of G is at the level of
$0.3\%$.

###  2.1 Details of data analysis

The whole analysis is based on linear regression (LR) [20, 21, 22, 23]. Only
terms with p-values below $0.4$ are kept in the linear regression, carried out
by using the MATLAB routine $fitlm$, with the option $RobustOpt$ on. It has
been checked that the final p-values are always below $0.2$. The
reconstruction of the AC and DC signals is based on Hilbert transform: data
are band-pass filtered around the beat note (in a range $\pm 12$ Hz around
$\omega_{m}$) and Blackman windowing is implemented. This part of the analysis
has been validated by processing a known ideal sinusoidal signal, without
noise addition. It has been checked that the Hilbert transform routine leads
to recover the expected frequency within 7 nHz before decimation, and within 5
nHz after decimation down to the rate of $0.016667$ Hz (60 s measurement
time). Relevant terms in the analysis are evaluated at 2 kHz, the rate of the
available data, and decimated down to $0.016667$ Hz rate, taking into account
the information provided by the G group.

Gas pressure | 10 mbar
---|---
Area of the beam at waist |  $0.51\times 0.74^{-6}$ m2
Mirror transmission | 0.2 ppm
Perimeter of the cavity | 16 m
Kinetic discharge temperature | 360 K
Photodiode quantum efficiency | 0.5
Trans-impedance amplifier gain |  $1.1\times 10^{8}$ V/A
Beat note mean value | 348.516 Hz
average loss $<\mu>$ | $6.51280\times 10^{-5}$
Scale factor |  $6.3211125158\times 10^{6}$ Hz s/rad
Table 1: Data of the G RLG set up.

A set of explanatory variables is evaluated: the laser explanatory variables
$\omega_{\xi}$, $\omega_{ns1}$, $\omega_{ns1}\times\delta\mu$, and the
available geodetic terms. The LR procedure has been repeated using different
schemes. In the first scheme the LR uses the whole set of terms all at once.
In the second scheme the LR has been applied in two steps: the first step
evaluates $\omega_{s}$ using only the laser terms, leading to an initial
evaluation of the laser corrections. In the second step the geodetic terms are
used as explanatory variables. When available, $FGEO$ is added as a known
signal and used for check. The different schemes produce always very similar
results: in the following, the first scheme is used, since it is the more
conservative one. Care has been put to avoid local minima in the regression,
which sometimes occurred.

###  2.2 The data

The two data sets have 3 full days stored at 2 kHz rate, day 99, 100, and 101,
year 2020, and day 121, 122, and 123, year 2022. The two data sets are
different and the 2020 one presents more short duration glitches than the
other. The G data, all the related analysis steps, the employed geophysical
models, and environmental monitors are stored on a daily basis in a file at
0.016667 Hz rate. Details of data associated to each column of the published
files can be found in the Appendix: in the following, the column numbers will
be used in the analysis description. Column 34 provides the measured Earth
rotation with the the estimated backscatter correction applied and the
theoretical diurnal polar motion and solid Earth tides model subtracted. Since
the absolute orientation of the G structure is not known well enough, these
latter corrections have to be considered as preliminary. Several tiltmeters
are used in order to establish the long term changes of the sensor
orientation. These in turn have to be corrected for temperature related
effects and from mass attraction derived from a global weather model [7]. In
particular the latter comes with several days delay, hence it can not be
corrected in realtime. The initial value of the ring laser orientation is
taken from a local survey, which provides values with a substantial error,
much larger than 5 ppm. Typical values are 49.144802 degrees for the latitude
and 12.87630 degrees for the longitude. While the analysis is based on the
methods developed for GINGERINO, there are several differences in the
experimental readout system. In GINGERINO the interferogram is the beat note
taken at one corner and the monobeams are both measured at the same
neighboring corner. For G, the available interferogram is the sum of the two
beams at the beam combiner, and the monobeams are each measured at a different
corner in order to avoid perturbations from back reflections. In this case the
$\omega_{\xi}$ terms are effective in correcting differences in the mirror
transmission, and other differences due to the measurement at different
corners.

##  3 Analysis results

The beat note $\omega_{m}$ is evaluated with our method and compared with the
one recorded in the published file, contained in column 3.

There are small differences in the standard deviation, which in the present
analysis turns out smaller by $3-4$ $\mu$Hz, probably due to the band-pass
filter around the beat note applied before the frequency reconstruction. This
is relevant to reduce the impact of the laser systematic terms $\omega_{K1}$
and $\omega_{K2}$ [18, 19], first and second order expansion, since the second
harmonic contribution, which is very difficult to model and subtract, is
effectively eliminated in this way.

Remarkably, in the present data sets irregular signals occur in the separation
between one and the subsequent day, in fact to keep constant the number of
data for each day some points are missing at the end of each daily file,
standard procedure of the Miniseed data retrieval. In the analysis the missing
points are replaced with zeroes.

###  3.1 Days 99-101, 2020

We have analysed the days 99, 100 and 101 starting from raw data. Published
files, containing one day data, are extracted from a continuous logging
sequence with the help of a MINISEED data format, and, to prevent loss of
precision, a set of data in a MATLAB readable format has been used.
Analogously to the procedure used for GINGERINO, data are band-passed around
$\omega_{m}=348$ Hz, in a $\pm 12$ Hz interval. The relevant parameters are
evaluated (AC and DC of the monobeams, relative phase of the mode $\epsilon$,
and the beat note $\omega_{m}$) at 2 kHz rate. Different terms of the
analysis, $\omega_{s0}$, $\omega_{\xi}$, and $\omega_{ns1}$ are evaluated
using the available information summarized in Table 1, including photodiode
gain, parameters of the G RLG (size and gas pressure), the measured average
losses $<{\mu}>$. Six $\omega_{\xi}$ terms are considered, to take into
account errors in the evaluation of $\epsilon$, $I_{S1,2}$, and $PH_{1,2}$ of
the monobeams.

The first step of the analysis evaluates $\omega_{s0}$ using Eq. 2.
Remarkably, $\omega_{s0}$ values remain rather stable within each single day,
but small discontinuities are evident between one day and the other, due to
the fact that some points around midnight are missing.

For this reason data around midnight have been eliminated from the analysis:
no other cuts have been applied to the data.

Figure 1 shows $\omega_{s}$, mean value subtracted, evaluated by the LR
procedure for the three days (red solid line), and the discontinuities between
different days are no longer evident. $\omega_{m}$, mean values subtracted, is
also shown (blue solid line): small differences with $\omega_{s}$ are seen.

Figure 1: Comparison of $\omega_{s}$ (red) and $\omega_{m}$ (blue) evaluated
on the 2020 three days data set: mean values are subtracted in both cases. The
2020 three days data set is considered

The top panel of Fig. 2 reports the $\omega_{ns1}$ corrections evaluated by
the LR procedure. The total geodetic components are plotted in the bottom
panel, as evaluated by the LR procedure without including Chandler wobble
data.

Figure 2: Sum of the laser terms (top), total geodetic components as evaluated
by the procedure (bottom), without including Chandler wobble data. The 2020
three days data set is considered.

###  3.2 Days 121-123, 2022

A second, more recent, set of three days (day 121, 122, and 123, year 2022)
has been also analysed. In this second data set the G RLG was operating under
ideal conditions and at roughly 20% lower beam power. The signals are cleaner,
although small discontinuities between different days are still present, but
fewer points have been eliminated around midnight to cure the problem. The
$\omega_{s}$ has been evaluated analogously to the first data set. Data
indicate that the signal is more stable compared with the previous set of
data, accordingly the null shift terms have little impact, but remain to be
meaningful for the Allan deviation analysis. Figure 3 and Fig. 1 are very
similar with each other.

Figure 3: Same of Fig. 1 for the 2022 three days data set.

##  4 Residuals and Allan deviation of the two data sets

The two data sets exhibit very similar results, however, more points have been
eliminated in the 2020 case, probably in response to the higher beam power
setting. The amount of data retained for the analysis is $93.4\%$ of the total
for the 2020 and $95.8\%$ of the total for the 2022 data set. All cuts had to
be applied around midnight. Figure 4 shows the distributions of the residuals
for both sets; they are very similar with respect to each other. This
indicates that the noise from the laser dynamics has been effectively
identified and subtracted by the procedure in both cases.

Figure 4: Distribution of the residuals for the two data sets expressed as
angular velocity. The standard deviations are $7.73\pm 0.05$ $\mu$Hz and
$7.14\pm 0.05$ $\mu$Hz for the 2020 and 2022 data sets, respectively.

Figure 5 shows the Allan444The function allan of Matlab has been used.
deviation of the residuals for both data sets, sampled at 0.016667 Hz. The
green solid lines display results obtained by applying the LR procedure at
once on the whole three days data sets. The procedure has been also applied to
each single day separately: orange solid lines report the corresponding
results (day 99, year 2020, and day 121, year 2022, are considered as
examples). Remarkably, plots cross the 1 part in $10^{9}$ threshold in less
than one day of integration time, albeit with a larger standard deviation.

The same procedure has been applied to each single day separately and the two
panels in Fig. 5 refer to the Allan deviation of the residuals. Curves
obtained in the present analysis are always below those coming out from the
standard one, with the exception of one of the days, around $\tau\simeq
8\times 10^{3}$ s. When $FGEO$ contributions are taken into account, the
present analysis leads to even smaller Allan deviation. In any case, plots
cross the 1 part in $10^{9}$ threshold with less than 1 day of integration
time.

Figure 5: LEFT: 2020, Allan deviation obtained according to the available
analyses, for the 3 days (green) and day 99 (orange). RIGHT: 2022, Allan
deviation obtained according to the available analyses, for the 3 days (green)
and day 121 (orange). The red dashed line in the right figure shows the Allan
deviation subtracting only the geophysical components from the beat note
$\omega_{m}$.

At the present stage we do not investigate the nature of the residuals, since
the focus is here to investigate whether the Allan deviation goes below the
fundamental physics threshold, i.e., below about 1 part in $10^{9}$, while
keeping the analysis as simple as possible. Moreover, owing to its monolithic
structure, G is certainly less prone to those mechanical coupling effects
accounted for, in GINGERINO analysis, by the extra term based on the product
of the residuals and the tiltmeter signal.

The geophysical components have been subtracted to $\omega_{m}$ using the LR,
without taking into account the laser terms $\omega_{\xi}$ and $\omega_{ns1}$.
With the 2022 data set, definitely less affected by laser dynamics, the Allan
deviation is a factor 2 worse at $10^{4}$ seconds compared to Fig. 5 RIGHT,
where the Allan shown with the dashed red line is super-imposed. With the 2020
data set, the result is even higher, indicating that, despite G is based on an
extremely rigid optical cavity, it is necessary to identify and subtract the
laser systematic terms in order to reach, and go beyond, the 1 part in
$10^{9}$ precision level of the Earth rotation rate measurement.

##  5 Conclusions

The analysis based on the model developed for GINGERINO has been extended to
the G RLG, rewriting the mathematical relationships with the main aim to
account for the characteristic values of the G apparatus. In this way,
corrections to the laser systematics have been subtracted in a deterministic
way following the model based on the stationary solution of the RLG equations
and the loss $\bar{\mu}$ of the optical cavity determined by statistical
means. The signals of geodetic origin have been always subtracted assuming
linear relationships to minimize the residuals via linear regression. For the
full set of three days of 2020, the Allan deviation drops below one part in
$10^{9}$ in less than 1 day of integration time.

The data point close to midnight are eliminated, to overcome a problem in the
data retrieval subroutine and the separate analysis of each single day
exhibits Allan deviations going below 1 part in $10^{9}$ more rapidly.

The analysis has been repeated with a second, and more recent, set of three
days in June 2022, where the ring laser was operated at lower power. In this
case, the monobeam signals are cleaner, and the cavity appears to be more
stable, accordingly the effects of $\omega_{ns1}$ are smaller, but still
remain significant for the lower Allan deviation outcome.

G is based on a monolithic structure in ZERODUR, a very low thermal expansion
material, with mirrors optically contacted to the structure. For this reason
the cavity is extremely stable, nevertheless tiny laser systematic effects are
present and it is necessary to subtract them in order to improve the
performance above the 1 part in $10^{9}$ level. Certainly laser systematics
cancellation is more relevant for RLGs based on a HL design, as GINGERINO and
ROMY, since in this case cavity losses, the quantity effectively ruling laser
systematics, is affected by variations of the mirror distances and their
relative alignment.

## References

  * [1] K. U. Schreiber and K.U. Wells, Large ring lasers for rotation sensing, Rev. Scient. Instr. 84, 041101 (2013).
  * [2] Heiner Igel et al.ROMY: a multicomponent ring laser for geodesy and geophysics, Geophysical Journal International, Volume 225, Issue 1, April 2021, Pages 684-698.
  * [3] Angela D. V. Di Virgilio et al.Underground Sagnac gyroscope with sub-prad/s rotation rate sensitivity Phys. Rev. Research 2, 032069(R) (2020).
  * [4] A. D. V. Di Virgilio et al., Sensitivity limit investigation of a Sagnac gyroscope through linear regression analysis Eur. Phys. J. C, 81, 400 (2021).
  * [5] D. Zou et al., Gyroscopic performance and some seismic measurements made with a 10 meter perimeter ring laser gyro, Appl. Opt. 60, 1737-1743 (2021).
  * [6] K. Liu et al., Large-scale passive laser gyroscope for earth rotation sensing , Optics Letters, 44, 11, 2732-2735 (2019).
  * [7] Klügel, T., Wziontek, H. Correcting gravimeters and tiltmeters for atmospheric mass attraction using operational weather models, Journal of Geodynamics 48, 204–210 (2009).
  * [8] K. U. Schreiber, T. Klügel, J.-P. R. Wells, R. B. Hurst, and A. Gebauer, How to detect the Chandler and the annual wobble of the Earth with a large ring laser gyroscope, Phys. Rev. Lett. 107, 173904 (2011).
  * [9] Robert W. Dunn, Dmitry E. Shabalin, Robert J. Thirkettle,et al.,Design and initial operation of a 367-m2 rectangular ring laser  Vol. 41, APPLIED OPTICS 1685, (2002).
  * [10] A. Gebauer et al.Reconstruction of the Instantaneous Earth Rotation Vector with Sub -Arcsecond Resolution Using a Large Scale Ring Laser Array Phys. Rev. Lett. 125, 033605 (2020)
  * [11] J. Belfi, N. Beverini, G. Carelli, A. Di Virgilio, U. Giacomelli, E. Maccioni, A. Simonelli, F. Stefani, and G. Terreni, Analysis of 90 day operation of the GINGERINO gyroscope, Appl. Opt. 57, 5844-5851 (2018)
  * [12] Basti, A., Beverini, N., Bosi, F. et al., Effects of temperature variations in high-sensitivity Sagnac gyroscope, Eur. Phys. J. Plus 136, 537 (2021).
  * [13] W. Chow et al., The ring laser gyro, Rev. Mod. Phys.57.61 (1985).
  * [14] L. N. Menegozzi and W. E. Lamb. Theory of a ring laser, Phys. Rev. A, 8:4, (1973).
  * [15] F. Aronowitz and R. J. Collins. Lock-in and intensity-phase interaction in the ring laser,. Journal of Applied Physics, 41:1, (1970).
  * [16] A. Beghi, et al.. Compensation of the laser parameter fluctuations in large ring-laser gyros: a kalman flter approach. Appl. Opt., 51(31):7518-7528, (2012).
  * [17] D. Cuccato, A. Beghi, J. Belfi, N. Beverini, A. Ortolan, and A. Di Virgilio.Controlling the nonlinear inter cavity dynamics of large he-ne laser gyroscopes. Metrologia, 51:97-107, (2014).
  * [18] A. D. V. Di Virgilio, N. Beverini, G. Carelli, D. Ciampini, F. Fuso, and E. Maccioni, Eur. Phys. J. C 79, 573 (2019).
  * [19] A. D. V. Di Virgilio, N. Beverini, et al., Eur. Phys. J. C 80, 163 (2020).
  * [20] Kay, Steven M., Fundamentals of Statistical Signal Processing: Estimation Theory, Prentice-Hall, Inc., US,(1993).
  * [21] J. Neter, M. H. Kutner, C. J. Nachtsheim, and W. Wasserman, Applied Linear Statistical Models, Irwin, Chicago (1996).
  * [22] A. Sen and M. Srivastava. Regression Analysis: Theory, Methods, and Applications, Springer Texts in Statistics. Springer New York (1997).
  * [23] F. Di Renzo, Characterisation and mitigation of non-stationary noise in Advance Gravitational Wave Detectors, PhD Thesis in Physics, University of Pisa, etd-06192020-081853 , June 2020.

##  Appendix A The G daily file with data and results

On a daily basis, G data sampled at 0.0166667 Hz rate and relevant results are
published in a file.Contents of each column are listed in Table 2.

Column | Quantity | Units
---|---|---
1 | Epoch | [MJD]
2 | Epoch | [day]
3 | Sagnac (single tone extractor) | [Hz]
4 | Sagnac rms | [Hz]
5 | Phase between SW-port and Sagnac beam combiner port | [rad]
6 | Phase between SE-port and Sagnac beam combiner port | [rad]
7 | SW-AC | [V]
8 | SW-DC | [V]
9 | SW-AC/DC | -
10 | SE-AC | [V]
11 | SE-DC | [V]
12 | SE-AC/DC | -
13 | SW-SE-Phase | [rad]
14 | estimated Sagnac correction factor | -
15 | estimated Sagnac correction value | [mHz]
16 | Sagnac BS-corrected | [Hz]
17 | SE-DDS-Amplitude AC-level for driving LED | [V]
18 | SW-DDS-Amplitude AC-level for driving LED | [V]
19 | SW-Phase to DDS driving LED | [rad]
20 | SE-Phase to DDS driving LED | [rad]
21 | Sagnac geophys. model contribution subtracted | [Hz]
| (Oppolzer-Terms $\&$ Tilt-NS deformation effect) |
22 | Sagnac of SE-monobeam estimation (Single Tone extractor) | [Hz]
23 | Sagnac of SW-monobeam estimation (Single Tone extractor) | [Hz]
24 | Geophysical Model: Oppolzer terms NS | [rad]
25 | Geophysical Model: Oppolzer terms EW | [rad]
26 | Geophysical Model: Theoretical tilt NS (attraction part) | [rad]
27 | Geophysical Model: Theoretical tilt NS (deformation part) | [rad]
28 | Observed tilt NS | [rad]
29 | Geophysical Model: Theoretical tilt EW (attraction part) | [rad]
30 | Geophysical Model: Theoretical tilt EW (deformation part) | [rad]
31 | Observed tilt EW | [rad]
32 | Geophysical Model: Sagnac geophys. model contribution | [$\mu$Hz]
| (Oppolzer-Terms $\&$ Tilt-NS deformation effect) |
33 | Sagnac geophys. model reduced | [Hz]
| (Oppolzer-Terms $\&$ Tilt-NS deformation effect) |
34 | Sagnac geophys. model and BS applied | [Hz]
35 | Pressure vessel barometric pressure | [hPa]
| (Paroscientific at pressure vessel supply) |
36 | Pressure vessel temperature | [oC]
37 | Pressure vessel humidity | [$\%$rH]
38 | Ringlaser room barometric pressure | [hPa]
| (Vaisala outside pressure vessel) |
39 | Ringlaser room temperature | [oC]
40 | Ringlaser room humidity | [$\%$rH]
41 | Control room temperature | [oC]
42 | Control room humidity | [$\%$rH]
43 | Plasma brightness | [V]
Table 2: Content of columns in the G data and results file.


The $B_{c}$ meson and its scalar cousin with the QCD sum rules

Zhi-Gang Wang<EMAIL_ADDRESS>
Department of Physics, North China Electric Power University, Baoding 071003,
P. R. China

###### Abstract

In this work, we use optical theorem to calculate the next-to-leading order
corrections to the spectral densities directly in the QCD sum rules for the
pseudoscalar and scalar $B_{c}$ mesons. We take the experimental data as
guides to perform updated analysis, and obtain the masses and decay constants,
therefore the leptonic decay widths, which can be confronted to the
experimental data in the future.

PACS number: 12.38.Bx, 12.38.Lg

Key words: Next-to-leading order contributions, QCD sum rules

##  1 Introduction

In 1998, the CDF collaboration observed the pseudoscalar $B_{c}$ mesons
through the decay modes $B_{c}^{\pm}\to J/\psi\ell^{\pm}X$ and $B_{c}^{\pm}\to
J/\psi\ell^{\pm}\bar{\nu}_{\ell}$ in the $p\bar{p}$ collisions at the energy
$\sqrt{s}=1.8\,\rm{TeV}$ at the Fermilab Tevatron, the measured mass is
$6.40\pm 0.39\pm 0.13\,\rm{GeV}$ [1, 2].

In 2007, the CDF collaboration confirmed the $B_{c}$ mesons through the decay
modes $B_{c}^{\pm}\to J/\psi\pi^{\pm}$ with the measured mass $6275.6\pm
2.9\pm 2.5\,\rm{MeV}$ [3]. In 2008, the D0 collaboration reconstructed the
$B_{c}^{\pm}\to J/\psi\pi^{\pm}$ decays and confirmed the $B_{c}$ mesons with
the measured mass $6300\pm 14\pm 5\,\rm{MeV}$ [4]. Now the average value
listed in the Review of Particle Physics is $6274.47\pm 0.27\pm
0.17\,\rm{MeV}$ [5].

In 2014, the ATLAS collaboration reported the observation of a structure in
the $B_{c}^{\pm}\pi^{+}\pi^{-}$ invariant mass spectrum, which is consistent
with the predicted $B_{c}^{\prime}$ meson with a mass of $6842\pm 4\pm
5\,\rm{MeV}$ [6].

In 2019, the CMS collaboration observed two excited $\bar{b}c$ states in the
$B^{+}_{c}\pi^{+}\pi^{-}$ invariant mass spectrum, which are consistent with
the $B^{\prime+}_{c}$ and $B_{c}^{*\prime+}$, respectively [7]. The two states
are separated in mass by $29.1\pm 1.5\pm 0.7\,\rm{MeV}$, and the mass of the
$B^{\prime+}_{c}$ is measured to be $6871.0\pm 1.2\pm 0.8\pm 0.8\,\rm{MeV}$.
Also in 2019, the LHCb collaboration observed the excited $B^{\prime+}_{c}$
and $B_{c}^{*\prime+}$ mesons in the $B^{+}_{c}\pi^{+}\pi^{-}$ invariant mass
spectrum. The $B_{c}^{*\prime+}$ meson has a mass of $6841.2\pm 0.6\pm 0.1\pm
0.8\,\rm{MeV}$, which is reconstructed without the low-energy photon emitted
in the $B_{c}^{*+}\to B_{c}^{+}\gamma$ decay following the process
$B_{c}^{\prime*+}\to B_{c}^{*+}\pi^{+}\pi^{-}$, while the $B^{\prime+}_{c}$
meson has a mass of $6872.1\pm 1.3\pm 0.1\pm 0.8\,\rm{MeV}$ [8].

It is odd that the $B_{c}^{\prime}$ meson emerges as heavier than the mass of
the $B_{c}^{*\prime}$ meson, which is in conflict with all the theoretical
estimations, this maybe or maybe not due to impossibility of reconstruction of
the low-energy photon in the $B_{c}^{\ast+}\rightarrow B_{c}^{+}\gamma$ decay
[9]. Only the $B_{c}$ and $B_{c}^{\prime}$ mesons are listed in Review of
Particle Physics [5], which are in contrary to the copious (well-established)
spectroscopy of the charmonium and bottomonium states. Despite the enormous
developments on the heavy quark physics in recent years, the bottom-charm
spectroscopy remains poorly known, which calls for further investigations.

The beauty-charm mesons provide an optimal laboratory for exploring both the
perturbative and nonperturbative dynamics of the heavy quarks, due to absence
of the light quark’s contamination, and for exploring the strong and electro-
weak interactions, as they are composed of two different heavy flavors and
cannot annihilate into gluons or photons. For the excited $c\bar{b}$ states,
which lie below the $BD$ threshold, would decay into the $B_{c}$ meson through
the radiative decays or hadronic decays [10, 11], while the ground state
$B_{c}$ can only decay weakly through emitting a virtual $W$-boson.

There have been several theoretical works on the mass spectroscopy of the
bottom-charm mesons, such as the relativized (or relativistic) quark model
with an special potential [10, 11, 12, 13, 14], the nonrelativistic quark
model with an special potential [15, 16, 17, 18, 19, 20, 21, 22, 23], the
semi-relativistic quark model using the shifted large-$N$ expansion [24, 25],
the perturbative QCD [26], the nonrelativistic renormalization group [27], the
lattice QCD [28, 29, 30, 31], the Bethe-Salpeter equation [32, 33, 34], the
full QCD sum rules [35, 36, 37, 38, 39, 40, 41], the potential model combined
with the QCD sum rules [15, 16], etc.

With the continuous development in experimental techniques, we expect that
more $c\bar{b}$ states would be observed in the future. The decay constant,
which parameterizes the coupling between a current and a meson, plays an
important role in exploring the exclusive processes, because the decay
constants are not only a fundamental parameter describing the pure leptonic
decays, but also are an universal input parameter related to the distribution
amplitudes, form-factors, partial decay widths and branching fractions in many
processes. By precisely measuring the branching fractions, we can resort to
the decay constants to extract the CKM matrix element in the standard model
and search for new physics beyond the standard model [42].

Decay constants of the bottom-charm mesons have been investigated in several
theoretical approaches, such as the full QCD sum rules [35, 36, 37, 38, 39,
40, 41, 43, 44], the potential model combined with the QCD sum rules [15, 16],
the QCD sum rule combined with the heavy quark effective theory [45, 46, 47,
48, 49, 50, 51, 52], the covariant light-front quark model [53, 54], the
lattice non-relativistic QCD [31], the shifted large-$N$ expansion method
[25], the field correlator method [55], etc. However, the values from
different theoretical approaches vary in a large range, it is interesting to
extend our previous works on the vector and axialvector $B_{c}$ mesons [40] to
investigate the pseudoscalar and scalar $B_{c}$ mesons with the full QCD sum
rules by including next-to-leading order radiative corrections and choose the
updated input parameters, thus our investigations are performed in a
consistent way. We take the experimental data as guides to choose the suitable
Borel parameters and continuum threshold parameters, examine the masses and
decay constants of the pseudoscalar and scalar $B_{c}$ mesons with the QCD sum
rules, therefore we calculate the leptonic decay widths to be confronted to
the experimental data in the future.

The article is arranged: we calculate the next-to-leading order contributions
to the spectral densities and obtain the QCD sum rules in Sect.2; in Sect.3,
we present the numerical results and discussions; and Sect.4 is reserved for
our conclusions.

##  2 Explicit calculations of the spectral densities at the next-to-leading
order

We write down the two-point correlation functions firstly,

$\displaystyle\Pi_{P/S}(p^{2})$ $\displaystyle=$ $\displaystyle i\int
d^{4}xe^{ip\cdot x}\langle 0|T\left\\{J(x)J^{\dagger}(0)\right\\}|0\rangle\,,$
(1)

where $J(x)=J_{P}(x)$ and $J_{S}(x)$,

$\displaystyle J_{P}(x)$ $\displaystyle=$
$\displaystyle\bar{c}(x)i\gamma_{5}b(x)\,,$ $\displaystyle J_{S}(x)$
$\displaystyle=$ $\displaystyle\bar{c}(x)b(x)\,.$ (2)

The correlation functions can be written in the form,

$\displaystyle\Pi_{P/S}(p^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{(m_{b}+m_{c})^{2}}^{\infty}ds\frac{{\rm
Im\Pi}_{P/S}(s)}{s-p^{2}}\,,$ (3)

through dispersion relation, where

$\displaystyle\frac{{\rm Im}\Pi_{P/S}(s)}{\pi}$ $\displaystyle=$
$\displaystyle\rho^{0}_{P/S}(s)+\rho^{1}_{P/S}(s)+\rho^{2}_{P/S}(s)+\cdots\,,$
(4)

the $\rho^{0}_{P/S}(s)$, $\rho^{1}_{P/S}(s)$, $\rho^{2}_{P/S}(s)$, $\cdots$
are the spectral densities of the leading order, next-to-leading order, and
next-to-next-to-leading order, $\cdots$. At the leading order,

$\displaystyle\rho^{0}_{P/S}(s)$ $\displaystyle=$
$\displaystyle\frac{3}{8\pi^{2}}\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}{s}\left[s-(m_{b}\mp
m_{c})^{2}\right]\,,$ (5)

where

$\displaystyle\lambda(s,m_{b}^{2},m_{c}^{2})$ $\displaystyle=$ $\displaystyle
s^{2}+m_{b}^{4}+m_{c}^{4}-2sm_{b}^{2}-2sm_{c}^{2}-2m_{b}^{2}m_{c}^{2}\,.$ (6)

At the next-to-leading order, there are three Feynman diagrams contribute to
the correlation functions, see Fig.1. We calculate the imaginary parts using
the Cutkosky’s rule or optical theorem, which lead to the same results, then
use dispersion relation to obtain the correlation functions [39, 56]. There
are ten possible cuts, six cuts attribute to virtual gluon emissions and four
cuts attribute to real gluon emissions.

Figure 1: The next-to-leading order contributions to the correlation
functions.  Figure 2: Six possible cuts correspond to virtual gluon emissions.
Figure 3: The quark self-energy correction.  Figure 4: The vertex correction.

The six cuts shown in Fig.2 make attributions to virtual gluon emissions, and
correspond to the self-energy and vertex corrections, respectively. We
calculate the Feynman diagrams directly using the dimensional regularization
to regularize both the ultraviolet and infrared divergences, and resort to the
on-shell renormalization scheme to subtract the divergences so as to
accomplish the wave-function and mass renormalizations. Then we take account
of all contributions shown in Fig.2 by the simple replacements of the vertexes
in the currents,

$\displaystyle\bar{u}(p_{1})i\gamma_{5}u(p_{2})$ $\displaystyle\rightarrow$
$\displaystyle\bar{u}(p_{1})i\gamma_{5}u(p_{2})+\bar{u}(p_{1})i\widetilde{\Gamma}_{5}u(p_{2})$
(7) $\displaystyle=$
$\displaystyle\sqrt{Z_{1}}\sqrt{Z_{2}}\bar{u}(p_{1})i\gamma_{5}u(p_{2})+\bar{u}(p_{1})i\Gamma_{5}u(p_{2})$
$\displaystyle=$
$\displaystyle\bar{u}(p_{1})i\gamma_{5}u(p_{2})\left(1+\frac{1}{2}\delta
Z_{1}+\frac{1}{2}\delta Z_{2}\right)+\bar{u}(p_{1})i\Gamma_{5}u(p_{2})\,,$
$\displaystyle\bar{u}(p_{1})u(p_{2})$ $\displaystyle\rightarrow$
$\displaystyle\bar{u}(p_{1})u(p_{2})+\bar{u}(p_{1})\widetilde{\Gamma}_{0}u(p_{2})$
(8) $\displaystyle=$
$\displaystyle\sqrt{Z_{1}}\sqrt{Z_{2}}\bar{u}(p_{1})u(p_{2})+\bar{u}(p_{1})\Gamma_{0}u(p_{2})$
$\displaystyle=$ $\displaystyle\bar{u}(p_{1})u(p_{2})\left(1+\frac{1}{2}\delta
Z_{1}+\frac{1}{2}\delta Z_{2}\right)+\bar{u}(p_{1})\Gamma_{0}u(p_{2})\,,$

where

$\displaystyle Z_{i}$ $\displaystyle=$ $\displaystyle 1+\delta
Z_{i}=1+\frac{4}{3}\frac{\alpha_{s}}{\pi}\left(-\frac{1}{4\varepsilon_{\rm
UV}}+\frac{1}{2\varepsilon_{\rm
IR}}+\frac{3}{4}\log\frac{m_{i}^{2}}{4\pi\mu^{2}}+\frac{3}{4}\gamma-1\right)\,,$
(9)

is the $i$ quark’s wave-function renormalization constant from the self-energy
correction, see Fig.3, and

$\displaystyle\Gamma_{5/0}$ $\displaystyle=$
$\displaystyle\gamma_{5}\frac{4}{3}g_{s}^{2}\int_{0}^{1}dx\int_{0}^{1-x}dy\int\frac{d^{D}k_{E}}{(2\pi)^{D}}$
(10)
$\displaystyle\frac{\Gamma(3)}{\left[k^{2}_{E}+(xp_{1}+yp_{2})^{2}\right]^{3}}\Big{\\{}4k^{2}_{E}\Big{(}1-\frac{1}{2}\varepsilon_{\rm
UV}\Big{)}+2(1-x-y+2xy)(s-m_{b}^{2}-m_{c}^{2})$ $\displaystyle\pm
2(x+y)m_{b}m_{c}+2x(1-2x)m_{b}^{2}+2y(1-2y)m_{c}^{2}\Big{\\}}\,,$

for the vertex corrections after the Wick’s rotation, see Fig.4, where
$\gamma$ is the Euler constant, $\mu$ is the renormalization scale, and the
Euclidean four-momentum $k_{E}=(k_{1},k_{2},k_{3},k_{4})$. We set the
dimension $D=4-2\varepsilon_{\rm UV}=4+2\varepsilon_{\rm IR}$ to regularize
the ultraviolet and infrared divergences respectively, and would add the scale
factors $\mu^{2\varepsilon_{\rm UV}}$ or $\mu^{-2\varepsilon_{\rm IR}}$ if
necessary.

We accomplish all the integrals and observe that the ultraviolet divergences
$\frac{1}{\varepsilon_{\rm UV}}$ in the $\Gamma_{5/0}$, $\delta Z_{1}$ and
$\delta Z_{2}$ are canceled out with each other, which is an outcome of the
Ward identity. The total contributions have no ultraviolet divergence,

$\displaystyle\widetilde{\Gamma}_{5}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\frac{\alpha_{s}}{4\pi}\gamma_{5}f_{P}(s)\,,$
$\displaystyle\widetilde{\Gamma}_{0}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\frac{\alpha_{s}}{4\pi}f_{S}(s)\,,$ (11)

where

$\displaystyle f_{P/S}(s)$ $\displaystyle=$
$\displaystyle\overline{f}_{P/S}(s)+\frac{2}{\varepsilon_{\rm
IR}}+3\log\frac{m_{b}m_{c}}{4\pi\mu^{2}}+4\log\frac{4\pi\mu^{2}}{s}-\gamma+4-\frac{2(s-m_{b}^{2}-m_{c}^{2})}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega}{1-\omega}\right)$
$\displaystyle\left(\frac{1}{\varepsilon_{\rm
IR}}+\log\frac{s}{4\pi\mu^{2}}+\gamma\right)\,,$
$\displaystyle\overline{f}_{P/S}(s)$ $\displaystyle=$ $\displaystyle
4\overline{V}(s)+2(s-m_{b}^{2}-m_{c}^{2})\left[\overline{V}_{00}(s)-V_{10}(s)-V_{01}(s)+2V_{11}(s)\right]\pm
2m_{b}m_{c}$
$\displaystyle\left[V_{10}(s)+V_{01}(s)\right]+2m_{b}^{2}\left[V_{10}(s)-2V_{20}(s)\right]+2m_{c}^{2}\left[V_{01}(s)-2V_{02}(s)\right]\,,$
$\displaystyle\omega$ $\displaystyle=$
$\displaystyle\sqrt{\frac{s-(m_{b}+m_{c})^{2}}{s-(m_{b}-m_{c})^{2}}}\,,$ (12)

and $s=p^{2}$, the definitions and explicit expressions of the notations
$\overline{V}(s)$, $\overline{V}_{00}(s)$ and $V_{ij}(s)$ with $i,j=0,1,2$ are
given in the appendix.

The total contributions of the virtual gluon emissions to imaginary parts of
the correlation functions are,

$\displaystyle\frac{{\rm Im}\Pi^{V}_{P/S}(s)}{\pi}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\frac{\alpha_{s}}{4\pi}\frac{6}{\pi}\int\frac{d^{D-1}\vec{p}_{1}}{(2\pi)^{D-1}2E_{p_{1}}}\frac{d^{D-1}\vec{p}_{2}}{(2\pi)^{D-1}2E_{p_{2}}}(2\pi)^{D}\delta^{D}(p-p_{1}-p_{2})f(s)\left[s-(m_{b}\mp
m_{c})^{2}\right]\,,$

the superscript $V$ denotes the virtual gluon emissions. We accomplish the
integrals straightforwardly in $D=4+2\varepsilon_{\rm IR}$ dimension as there
is no ultraviolet divergence, and obtain the analytical expressions,

$\displaystyle\frac{{\rm Im}\Pi^{V}_{P/S}(s)}{\pi}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\frac{\alpha_{s}}{\pi}\rho_{P/S}^{0}(s)\left\\{\frac{1}{\varepsilon_{\rm
IR}}-2\log
4\pi+\frac{1}{2}\gamma+\frac{1}{2}\log\frac{\lambda^{2}(s,m_{b}^{2},m_{c}^{2})m_{b}^{3}m_{c}^{3}}{\mu^{8}s^{3}}+\frac{1}{2}\overline{f}_{P/S}(s)\right.$
$\displaystyle\left.-\frac{s-m_{b}^{2}-m_{c}^{2}}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega}{1-\omega}\right)\left[\frac{1}{\varepsilon_{\rm
IR}}-2\log
4\pi+2\gamma-2+\log\frac{\lambda(s,m_{b}^{2},m_{c}^{2})}{\mu^{4}}\right]\right\\}\,.$

Figure 5: Four possible cuts correspond to real gluon emissions.

The four cuts shown in Fig.5 make contributions to the real gluon emissions,
the corresponding scattering amplitudes are shown explicitly in Fig.6. From
Fig.6, we write down the scattering amplitudes $T^{a}_{5,\alpha}(p)$ and
$T^{a}_{0,\alpha}(p)$ ,

$\displaystyle T^{a}_{5,\alpha}(p)$ $\displaystyle=$
$\displaystyle\bar{u}(p_{1})\left\\{ig_{s}\frac{\lambda^{a}}{2}\gamma_{\alpha}\frac{i}{\\!\not\\!{p}_{1}+\\!\not\\!{k}-m_{b}}i\gamma_{5}+i\gamma_{5}\frac{i}{-\\!\not\\!{p}_{2}-\\!\not\\!{k}-m_{c}}ig_{s}\frac{\lambda^{a}}{2}\gamma_{\alpha}\right\\}v(p_{2})\,,$
$\displaystyle T^{a}_{0,\alpha}(p)$ $\displaystyle=$
$\displaystyle\bar{u}(p_{1})\left\\{ig_{s}\frac{\lambda^{a}}{2}\gamma_{\alpha}\frac{i}{\\!\not\\!{p}_{1}+\\!\not\\!{k}-m_{b}}+\frac{i}{-\\!\not\\!{p}_{2}-\\!\not\\!{k}-m_{c}}ig_{s}\frac{\lambda^{a}}{2}\gamma_{\alpha}\right\\}v(p_{2})\,,$
(15)

where $\lambda^{a}$ is the Gell-Mann matrix. Then we obtain the contributions
to the imaginary parts of the correlation functions with optical theorem,

$\displaystyle\frac{{\rm Im}\Pi^{R}_{P/S}(s)}{\pi}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\int\frac{d^{D-1}\vec{k}}{(2\pi)^{D-1}2E_{k}}\frac{d^{D-1}\vec{p}_{1}}{(2\pi)^{D-1}2E_{p_{1}}}\frac{d^{D-1}\vec{p}_{2}}{(2\pi)^{D-1}2E_{p_{2}}}(2\pi)^{D}\delta^{D}(p-k-
p_{1}-p_{2})$ (16) $\displaystyle{\rm
Tr}\left\\{T^{a}_{5/0,\alpha}(p)T^{a\dagger}_{5/0,\beta}(p)\right\\}g^{\alpha\beta}$
$\displaystyle=$
$\displaystyle-\frac{2g_{s}^{2}}{\pi}\int\frac{d^{D-1}\vec{k}}{(2\pi)^{D-1}2E_{k}}\frac{d^{D-1}\vec{p}_{1}}{(2\pi)^{D-1}2E_{p_{1}}}\frac{d^{D-1}\vec{p}_{2}}{(2\pi)^{D-1}2E_{p_{2}}}(2\pi)^{D}\delta^{D}(p-k-
p_{1}-p_{2})$ $\displaystyle\left\\{2\left[s-(m_{b}\mp
m_{c})^{2}\right]\left[\frac{m_{b}^{2}}{(k\cdot
p_{1})^{2}}+\frac{m_{c}^{2}}{(k\cdot
p_{2})^{2}}-\frac{s-m_{b}^{2}-m_{c}^{2}}{k\cdot p_{1}k\cdot
p_{2}}\right.\right.$ $\displaystyle\left.\left.+\frac{s-K^{2}}{k\cdot
p_{1}k\cdot p_{2}}\right]-\frac{(s-K^{2})^{2}}{k\cdot p_{1}k\cdot
p_{2}}\right\\}\,,$

where we have used the formulas $\sum
u(p_{1})\bar{u}(p_{1})=\\!\not\\!{p}_{1}+m_{b}$ and $\sum
v(p_{2})\bar{v}(p_{2})=\\!\not\\!{p}_{2}-m_{c}$ for the quark and antiquark
respectively, $K^{2}=(p_{1}+p_{2})^{2}$, and introduce the superscript $R$ to
denote the real gluon emissions. We accomplish the integrals in
$D=4+2\varepsilon_{\rm IR}$ dimension as there is only infrared divergence,
and obtain the contributions,

$\displaystyle\frac{{\rm Im}\Pi^{R}_{P/S}(s)}{\pi}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\frac{\alpha_{s}}{\pi}\rho_{P/S}^{0}(s)\left\\{-\frac{1}{\varepsilon_{\rm
IR}}+2\log
4\pi-2\gamma+2-\log\frac{\lambda^{3}(s,m_{b}^{2},m_{c}^{2})}{m_{b}^{2}m_{c}^{2}s^{2}\mu^{4}}+(s-m_{b}^{2}-m_{c}^{2})\overline{R}_{12}(s)\right.$
(17)
$\displaystyle-\overline{R}_{11}(s)-\overline{R}_{22}(s)-R_{12}^{1}(s)+\frac{R_{12}^{2}}{2}\frac{1}{s-(m_{b}\mp
m_{c})^{2}}+\frac{s-m_{b}^{2}-m_{c}^{2}}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}$
$\displaystyle\left.\log\left(\frac{1+\omega}{1-\omega}\right)\left[\frac{1}{\varepsilon_{\rm
IR}}-2\log
4\pi+2\gamma-2+\log\frac{\lambda^{3}(s,m_{b}^{2},m_{c}^{2})}{m_{b}^{2}m_{c}^{2}s^{2}\mu^{4}}\right]\right\\}\,,$

the definitions and explicit expressions of the $\overline{R}_{11}(s)$,
$\overline{R}_{22}(s)$, $\overline{R}_{12}(s)$, $R_{12}^{1}(s)$ and
$R^{2}_{12}(s)$ are given in the appendix.

Figure 6: The amplitudes for the real gluon emissions.

Now we obtain the total spectral densities at the next-to-leading order,

$\displaystyle\rho_{P/S}^{1}(s)$ $\displaystyle=$
$\displaystyle\frac{4}{3}\frac{\alpha_{s}}{\pi}\rho_{P/S}^{0}(s)\left\\{\frac{1}{2}\overline{f}_{P/S}(s)-\overline{R}_{11}(s)-\overline{R}_{22}(s)-R_{12}^{1}(s)+(s-m_{b}^{2}-m_{c}^{2})\,\overline{R}_{12}(s)\right.$
(18) $\displaystyle\left.+\frac{R_{12}^{2}}{2}\frac{1}{s-(m_{b}\mp
m_{c})^{2}}-\frac{3}{2}\gamma+2+\frac{1}{2}\log\frac{m_{b}^{7}m_{c}^{7}s}{\lambda^{4}(s,m_{b}^{2},m_{c}^{2})}\right.$
$\displaystyle\left.+\frac{s-m_{b}^{2}-m_{c}^{2}}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega}{1-\omega}\right)\log\frac{\lambda^{2}(s,m_{b}^{2},m_{c}^{2})}{m_{b}^{2}m_{c}^{2}s^{2}}\right\\}\,.$

The infrared divergences $\frac{1}{\varepsilon_{\rm IR}}$,
$\log\left(\frac{1+\omega}{1-\omega}\right)\frac{1}{\varepsilon_{\rm IR}}$
from the virtual and real gluon emissions are canceled out with each other,
which is guaranteed by the Lee-Nauenberg theorem [57]. The analytical
expressions are applicable in many phenomenological analysis besides the QCD
sum rules.

Then we calculate the contributions of the gluon condensate directly, the
calculations are easy and no much to say. Finally, we obtain the analytical
expressions of the spectral densities, take the quark-hadron duality below the
continuum thresholds $s^{0}_{P/S}$ and perform the Borel transforms in regard
to the variable $P^{2}=-p^{2}$ to obtain the QCD sum rules,

$\displaystyle\frac{f_{P/S}^{2}M_{P/S}^{4}}{(m_{b}\pm
m_{c})^{2}}\exp\left(-\frac{M_{P/S}^{2}}{T^{2}}\right)$ $\displaystyle=$
$\displaystyle\int_{(m_{b}+m_{c})^{2}}^{s^{0}_{P/S}}ds\left[\rho^{0}_{P/S}(s)+\rho^{1}_{P/S}(s)+\rho^{\rm
con}_{P/S}(s)\right]\exp\left(-\frac{s}{T^{2}}\right)\,,$

where

$\displaystyle\rho^{\rm con}_{P/S}(s)$ $\displaystyle=$
$\displaystyle\mp\frac{m_{b}m_{c}}{24T^{4}}\langle\frac{\alpha_{s}GG}{\pi}\rangle\int_{0}^{1}dx\left[\frac{m_{c}^{2}}{x^{3}}+\frac{m_{b}^{2}}{(1-x)^{3}}\right]\delta(s-\widetilde{m}_{Q}^{2})$
(20)
$\displaystyle\pm\frac{m_{b}m_{c}}{8T^{2}}\langle\frac{\alpha_{s}GG}{\pi}\rangle\int_{0}^{1}dx\left[\frac{1}{x^{2}}+\frac{1}{(1-x)^{2}}\right]\delta(s-\widetilde{m}_{Q}^{2})$
$\displaystyle-\frac{s}{24T^{4}}\langle\frac{\alpha_{s}GG}{\pi}\rangle\int_{0}^{1}dx\left[\frac{(1-x)m_{c}^{2}}{x^{2}}+\frac{xm_{b}^{2}}{(1-x)^{2}}\right]\delta(s-\widetilde{m}_{Q}^{2})\,,$

$\widetilde{m}_{Q}^{2}=\frac{m_{b}^{2}}{1-x}+\frac{m_{c}^{2}}{x}$, the $T^{2}$
is the Borel parameter, and the decay constants are defined by,

$\displaystyle\langle 0|J_{P}(0)|P(p)\rangle$ $\displaystyle=$
$\displaystyle\frac{f_{P}M_{P}^{2}}{m_{b}+m_{c}}\,,$ $\displaystyle\langle
0|J_{S}(0)|S(p)\rangle$ $\displaystyle=$
$\displaystyle\frac{f_{S}M_{S}^{2}}{m_{b}-m_{c}}\,,$ (21)

in other words,

$\displaystyle\langle 0|J^{\alpha}_{A}(0)|P(p)\rangle$ $\displaystyle=$
$\displaystyle if_{P}\,p^{\alpha}\,,$ $\displaystyle\langle
0|J_{V}^{\alpha}(0)|S(p)\rangle$ $\displaystyle=$ $\displaystyle
if_{S}\,p^{\alpha}\,,$ (22)

the subscripts $A$ and $V$ denote the axial-vector and vector currents,
respectively.

We eliminate the decay constants $f_{P/S}$ and obtain the QCD sum rules for
the masses of the pseudoscalar and scalar $B_{c}$ mesons,

$\displaystyle M_{P/S}^{2}$ $\displaystyle=$
$\displaystyle\frac{\int_{(m_{b}+m_{c})^{2}}^{s^{0}_{P/S}}ds\frac{d}{d\left(-1/T^{2}\right)}\left[\rho^{0}_{P/S}(s)+\rho^{1}_{P/S}(s)+\rho^{\rm
con}_{P/S}(s)\right]\exp\left(-\frac{s}{T^{2}}\right)}{\int_{(m_{b}+m_{c})^{2}}^{s^{0}_{P/S}}ds\left[\rho^{0}_{P/S}(s)+\rho^{1}_{P/S}(s)+\rho^{\rm
con}_{P/S}(s)\right]\exp\left(-\frac{s}{T^{2}}\right)}\,.$ (23)

##  3 Numerical results and discussions

The value of the gluon condensate $\langle\frac{\alpha_{s}GG}{\pi}\rangle$ has
been updated from time to time, and changes greatly, we use the updated value
$\langle\frac{\alpha_{s}GG}{\pi}\rangle=0.022\pm 0.004\,\rm{GeV}^{4}$ [58]. We
take the $\overline{MS}$ masses of the heavy quarks $m_{c}(m_{c})=1.275\pm
0.025\,\rm{GeV}$ and $m_{b}(m_{b})=4.18\pm 0.03\,\rm{GeV}$ from the Particle
Data Group [5]. In addition, we take account of the energy-scale dependence of
the $\overline{MS}$ masses,

$\displaystyle m_{Q}(\mu)$ $\displaystyle=$ $\displaystyle
m_{Q}(m_{Q})\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}(m_{Q})}\right]^{\frac{12}{33-2n_{f}}}\,,$
$\displaystyle\alpha_{s}(\mu)$ $\displaystyle=$
$\displaystyle\frac{1}{b_{0}t}\left[1-\frac{b_{1}}{b_{0}^{2}}\frac{\log
t}{t}+\frac{b_{1}^{2}(\log^{2}{t}-\log{t}-1)+b_{0}b_{2}}{b_{0}^{4}t^{2}}\right]\,,$
(24)

where $t=\log\frac{\mu^{2}}{\Lambda^{2}}$, $b_{0}=\frac{33-2n_{f}}{12\pi}$,
$b_{1}=\frac{153-19n_{f}}{24\pi^{2}}$,
$b_{2}=\frac{2857-\frac{5033}{9}n_{f}+\frac{325}{27}n_{f}^{2}}{128\pi^{3}}$,
$\Lambda=213\,\rm{MeV}$, $296\,\rm{MeV}$ and $339\,\rm{MeV}$ for the quark
flavor numbers $n_{f}=5$, $4$ and $3$, respectively [5]. We choose $n_{f}=4$
and $5$ for the $c$ and $b$ quarks, respectively, and then evolve the heavy
quark masses to the typical energy scale $\mu=2\,\rm{GeV}$.

The threshold $(m_{b}+m_{c})^{2}$ decreases quickly with increase of the
energy scale, the energy scale should be larger than $1.7\,\rm{GeV}$, which
corresponds to the squared mass of the $B_{c}$ meson, $39.4\,\rm{GeV}^{2}$. If
we take the typical energy scale $\mu=2\,\rm{GeV}$, which corresponds to the
threshold $(m_{b}+m_{c})^{2}\approx 36.0\,{\rm{GeV}}^{2}<M_{P}^{2}$, it is
reasonable.

The experimental masses of the $B_{c}$ and $B_{c}^{\prime}$ mesons are
$6274.47\pm 0.27\pm 0.17\,\rm{MeV}$ and $6871.2\pm 1.0\,\rm{MeV}$ respectively
from the Particle Data Group [5]. The scalar $B_{c}$ meson still escapes the
experimental detection, roughly speaking, the theoretical mass is $6712\pm
18\pm 7\,\rm{MeV}$ from the lattice QCD [30] or $6714\,\rm{MeV}$ from the
nonrelativistic quark model [23]. We can tentatively take the continuum
threshold parameters as $s^{0}_{P}=(39-47)\,\rm{GeV}^{2}$ and
$s^{0}_{S}=(45-55)\,\rm{GeV}^{2}$, and search for the ideal values by assuming
the energy gap between the ground state and first radial excited states is
about $0.6\,\rm{GeV}$ if lacking experimental data.

After trial and error, we obtain the ideal Borel windows and continuum
threshold parameters, and the corresponding pole contributions about
$(70-85)\%$, the pole dominance is well satisfied. On the other hand, the
gluon condensate plays a tiny important role, the operator product expansion
is well convergent. It is reliable to extract the masses and pole residues,
which are shown in Table 1 and Figs.7-8.

The predicted mass $M_{P}=6.274\pm 0.054\,\rm{GeV}$ is in very good agreement
with the experimental data $6274.47\pm 0.27\pm 0.17\,\rm{MeV}$ from the
Particle Data Group [5], while the predicted mass $M_{S}=6.702\pm
0.060\,\rm{GeV}$ is consistent with other theoretical calculations [10, 11,
12, 13, 14, 15, 16, 18, 19, 20, 21, 22, 23, 28, 29, 30, 32, 33, 34].

| $T^{2}(\rm{GeV}^{2})$ | $s_{0}(\rm{GeV}^{2})$ | pole | $M(\rm{GeV})$ | $f(\rm{GeV})$
---|---|---|---|---|---
$B_{c}({0}^{-})$ | $3.0-4.0$ | $44\pm 1$ | $(68-89)\%$ | $6.274\pm 0.054$ | $0.371\pm 0.037$
$B_{c}({0}^{+})$ | $5.4-6.4$ | $54\pm 1$ | $(69-83)\%$ | $6.702\pm 0.060$ | $0.236\pm 0.017$
Table 1:  The Borel parameters, continuum threshold parameters, pole
contributions, masses and decay constants of the pseudoscalar and scalar
$B_{c}$ mesons.

Figure 7:  The masses of the pseudoscalar ($P$) and scalar ($S$) $B_{c}$
mesons with variations of the Borel parameters $T^{2}$.

Figure 8:  The pole residues of the pseudoscalar ($P$) and scalar ($S$)
$B_{c}$ mesons with variations of the Borel parameters $T^{2}$.

While the decay constants, even for the pseudoscalar $B_{c}$ meson, the
theoretical values vary in a large range, for example, $f_{P}=460\pm
60\,\rm{MeV}$ [16], $300\pm 65\,\rm{MeV}$ [36], $360\pm 60\,\rm{MeV}$ [37],
$270\pm 30\,\rm{MeV}$ [41], $371\pm 17\,\rm{MeV}$ [43], $528\pm 19\,\rm{MeV}$
[44] from the full QCD sum rules, $410\pm 40\,\rm{MeV}$ [11], $433\,\rm{MeV}$
[12] from the relativistic quark model, $498\,\rm{MeV}$ [18], $440\,\rm{MeV}$
[19] from the non-relativistic quark model, $523\pm 62\,\rm{MeV}$ [54] from
the light-front quark model, $420\pm 13\,\rm{MeV}$ [31] from the lattice non-
relativistic QCD, $315^{+26}_{-50}\,\rm{MeV}$ [25] from the shifted
$N$-expansion method, $438\pm 10\,\rm{MeV}$ [55] from the field correlator
method, $322\pm 42\,\rm{MeV}$ [33] from the Bethe-Salpeter equation, etc. At
the present time, it is difficult to say which value is superior to others.

The present calculation $f_{P}=371\pm 37\,\rm{MeV}$ is in very good agreement
with the value $371\pm 17\,\rm{MeV}$ from the full QCD sum rules [43]. In our
previous work, we obtain the values $f_{V}=384\pm 32\,\rm{MeV}$ and
$f_{A}=373\pm 25\,\rm{MeV}$ for the vector and axial-vector $B_{c}$ mesons,
respectively. Our calculations indicate that $f_{P}\approx f_{V}\approx
f_{A}>f_{S}$. While in the QCD sum rule combined with the heavy quark
effective theory up to the order $\alpha_{s}^{3}$, the decay constants have
the relations $\tilde{f}_{P}=f_{P}>f_{V}>f_{S}>\tilde{f}_{S}>f_{A}$ [47],
where the decay constants $\tilde{f}_{P}$ and $\tilde{f}_{S}$ are defined by

$\displaystyle\langle 0|J_{P}(0)|P(p)\rangle$ $\displaystyle=$
$\displaystyle\tilde{f}_{P}M_{P}\,,$ $\displaystyle\langle
0|J_{S}(0)|S(p)\rangle$ $\displaystyle=$ $\displaystyle\tilde{f}_{S}M_{S}\,.$
(25)

From Eq.(2) and Eq.(3), we can obtain the relations,

$\displaystyle\tilde{f}_{P}$ $\displaystyle=$ $\displaystyle
f_{P}\frac{M_{P}}{m_{b}+m_{c}}\,,$ $\displaystyle\tilde{f}_{S}$
$\displaystyle=$ $\displaystyle f_{S}\frac{M_{S}}{m_{b}-m_{c}}\,,$ (26)

it is obvious that $\tilde{f}_{P}>f_{P}$ and $\tilde{f}_{S}<f_{S}$, which are
in contrary to the relations obtained in Ref.[47], so no definite conclusion
can be obtained. Naively, we expect that the vector mesons have larger decay
constants than the corresponding pseudoscalar mesons [59].

The leptonic decay widths $\Gamma_{\ell\bar{\nu}_{\ell}}$ of the pseudoscalar
and scalar $B_{c}$ mesons can be written as,

$\displaystyle\Gamma_{\ell\bar{\nu}_{\ell}}$ $\displaystyle=$
$\displaystyle\frac{G_{F}^{2}}{8\pi}|V_{bc}|^{2}f_{P/S}^{2}M_{P/S}M_{\ell}^{2}\left(1-\frac{M_{\ell}^{2}}{M_{P/S}^{2}}\right)^{2}\,,$
(27)

where $\ell=e,\mu,\tau$, the Fermi constant $G_{F}=1.16637\times
10^{-5}\,\rm{GeV}^{-2}$, the CKM matrix element $V_{cb}=40.8\times 10^{-3}$,
the masses of the leptons $m_{e}=0.511\times 10^{-3}\,\rm{GeV}$,
$m_{\mu}=105.658\times 10^{-3}\,\rm{GeV}$, $m_{\tau}=1776.86\times
10^{-3}\,\rm{GeV}$, the life time of the $B_{c}$ meson
$\tau_{B_{c}}=0.510\times 10^{-12}\,s$ [5]. We take the masses and decay
constants of the pseudoscalar and scalar $B_{c}$ mesons from the QCD sum rules
to obtain the partial decay widths,

$\displaystyle\Gamma_{P\to e\bar{\nu}_{e}}$ $\displaystyle=$ $\displaystyle
2.03\times 10^{-12}\,\rm{eV}\,,$
$\displaystyle\Gamma_{P\to\mu\bar{\nu}_{\mu}}$ $\displaystyle=$ $\displaystyle
8.68\times 10^{-8}\,\rm{eV}\,,$
$\displaystyle\Gamma_{P\to\tau\bar{\nu}_{\tau}}$ $\displaystyle=$
$\displaystyle 2.08\times 10^{-5}\,\rm{eV}\,,$ $\displaystyle\Gamma_{S\to
e\bar{\nu}_{e}}$ $\displaystyle=$ $\displaystyle 8.78\times
10^{-13}\,\rm{eV}\,,$ $\displaystyle\Gamma_{S\to\mu\bar{\nu}_{\mu}}$
$\displaystyle=$ $\displaystyle 3.75\,\times 10^{-8}\,\rm{eV}\,,$
$\displaystyle\Gamma_{S\to\tau\bar{\nu}_{\tau}}$ $\displaystyle=$
$\displaystyle 9.18\,\times 10^{-6}\,\rm{eV}\,,$ (28)

and the branching fractions,

$\displaystyle{\rm Br}_{P\to e\bar{\nu}_{e}}$ $\displaystyle=$ $\displaystyle
1.57\times 10^{-9}\,,$ $\displaystyle{\rm Br}_{P\to\mu\bar{\nu}_{\mu}}$
$\displaystyle=$ $\displaystyle 6.73\times 10^{-5}\,,$ $\displaystyle{\rm
Br}_{P\to\tau\bar{\nu}_{\tau}}$ $\displaystyle=$ $\displaystyle 1.61\times
10^{-2}\,.$ (29)

The largest branching fractions of the $B_{c}(0^{-})\to\ell\bar{\nu}_{\ell}$
are of the order $10^{-2}$, the tiny branching fractions maybe escape
experimental detections.

##  4 Conclusion

In this work, we extend our previous works on the vector and axialvector
$B_{c}$ mesons to investigate the pseudoscalar and scalar $B_{c}$ mesons with
the full QCD sum rules by including next-to-leading order corrections and
choose the updated input parameters. In calculating the next-to-leading order
corrections, we use optical theorem to obtain the spectral densities directly,
and resort to the scheme of the dimensional regularization to regularize both
the ultraviolet and infrared divergences, which are canceled out with each
other separately, the net spectral densities are free of divergences. We take
the experimental data as guides to choose the suitable Borel parameters and
continuum threshold parameters, and make reasonable predictions for the masses
and decay constants, therefore the pure leptonic decay widths, which can be
confronted to the experimental data in the future.

## Acknowledgements

This work is supported by National Natural Science Foundation, Grant Number
12175068.

## Appendix

Firstly, we write down the fundamental integrals in calculating the vertex
corrections,

$\displaystyle V_{ab}(s)$ $\displaystyle=$ $\displaystyle
16\pi^{2}\int_{0}^{1}dx\int_{0}^{1-x}dy\int\frac{d^{D}k_{E}}{(2\pi)^{D}}\frac{x^{a}y^{b}\Gamma(3)}{\left[k^{2}_{E}+(xp_{1}+yp_{2})^{2}\right]^{3}}\,,$
$\displaystyle V(s)$ $\displaystyle=$ $\displaystyle
16\pi^{2}\left(1-\frac{1}{2}\varepsilon_{\rm
UV}\right)\int_{0}^{1}dx\int_{0}^{1-x}dy\int\frac{d^{D}k_{E}}{(2\pi)^{D}}\frac{k^{2}_{E}\Gamma(3)}{\left[k^{2}_{E}+(xp_{1}+yp_{2})^{2}\right]^{3}}\,,$
(30)

and carry out the integrals to obtain the analytical expressions,

$\displaystyle V_{00}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\left\\{-\log\left(\frac{1+\omega}{1-\omega}\right)\left(\frac{1}{\varepsilon_{\rm
IR}}+\log\frac{s}{4\pi\mu^{2}}+\gamma\right)+\frac{\log^{2}(1-\omega_{1}^{2})}{4}-\log^{2}(1+\omega_{1})\right.$
$\displaystyle+\frac{\log^{2}(1-\omega_{2}^{2})}{4}-\log^{2}(1+\omega_{2})+2\log(\omega_{1}+\omega_{2})\log\left(\frac{1+\omega}{1-\omega}\right)-\log\omega_{1}\log\left(\frac{1+\omega_{2}}{1-\omega_{2}}\right)$
$\displaystyle\left.-\log\omega_{2}\log\left(\frac{1+\omega_{1}}{1-\omega_{1}}\right)-{\rm
Li_{2}}\left(\frac{2\omega_{1}}{1+\omega_{1}}\right)-{\rm
Li_{2}}\left(\frac{2\omega_{2}}{1+\omega_{2}}\right)+\pi^{2}\right\\}\,,$
$\displaystyle=$
$\displaystyle\overline{V}_{00}(s)-\frac{1}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega}{1-\omega}\right)\left(\frac{1}{\varepsilon_{\rm
IR}}+\log\frac{s}{4\pi\mu^{2}}+\gamma\right)\,,$ $\displaystyle V_{10}(s)$
$\displaystyle=$
$\displaystyle\frac{1}{s}\left\\{\frac{1}{2}\log\left(\frac{1-\omega_{1}^{2}}{1-\omega_{2}^{2}}\right)-\frac{1}{\omega_{2}}\log\left(\frac{1+\omega}{1-\omega}\right)+\log\frac{\omega_{2}}{\omega_{1}}\right\\}\,,$
$\displaystyle V_{01}(s)$ $\displaystyle=$ $\displaystyle
V_{10}(s)|_{\omega_{1}\leftrightarrow\omega_{2}}\,,$ $\displaystyle V_{20}(s)$
$\displaystyle=$
$\displaystyle\frac{1}{2s}\left\\{-\frac{\omega_{1}\omega_{2}}{\omega_{1}+\omega_{2}}\log\left(\frac{1+\omega}{1-\omega}\right)-\frac{\omega_{1}}{\omega_{2}(\omega_{1}+\omega_{2})}\log\left(\frac{1+\omega}{1-\omega}\right)+\frac{\omega_{1}}{\omega_{1}+\omega_{2}}\log\left(\frac{1-\omega_{1}^{2}}{1-\omega_{2}^{2}}\right)\right.$
$\displaystyle\left.+\frac{2\omega_{1}}{\omega_{1}+\omega_{2}}\log\frac{\omega_{2}}{\omega_{1}}+1\right\\}\,,$
$\displaystyle V_{02}(s)$ $\displaystyle=$ $\displaystyle
V_{20}(s)|_{\omega_{1}\leftrightarrow\omega_{2}}\,,$ $\displaystyle V_{11}(s)$
$\displaystyle=$
$\displaystyle\frac{1}{2s}\left\\{\frac{\omega_{1}\omega_{2}}{\omega_{1}+\omega_{2}}\log\left(\frac{1+\omega}{1-\omega}\right)-\frac{\omega_{1}-\omega_{2}}{2(\omega_{1}+\omega_{2})}\log\left(\frac{1-\omega_{1}^{2}}{1-\omega_{2}^{2}}\right)-\frac{1}{\omega_{1}+\omega_{2}}\log\left(\frac{1+\omega}{1-\omega}\right)\right.$
$\displaystyle\left.+\frac{\omega_{1}}{\omega_{1}+\omega_{2}}\log\frac{\omega_{1}}{\omega_{2}}+\frac{\omega_{2}}{\omega_{1}+\omega_{2}}\log\frac{\omega_{2}}{\omega_{1}}-1\right\\}\,,$
$\displaystyle V(s)$ $\displaystyle=$ $\displaystyle\frac{1}{\varepsilon_{\rm
UV}}+\log\frac{4\pi\mu^{2}}{s}-\gamma+2-\frac{2\omega_{1}\omega_{2}}{\omega_{1}+\omega_{2}}\log\left(\frac{1+\omega}{1-\omega}\right)-\frac{\omega_{2}}{\omega_{1}+\omega_{2}}\log(1-\omega_{1}^{2})$
(31)
$\displaystyle-\frac{\omega_{1}}{\omega_{1}+\omega_{2}}\log(1-\omega_{2}^{2})-2\frac{\omega_{1}\log\omega_{1}+\omega_{2}\log\omega_{2}}{\omega_{1}+\omega_{2}}+2\log(\omega_{1}+\omega_{2})\,,$
$\displaystyle=$ $\displaystyle\overline{V}(s)+\frac{1}{\varepsilon_{\rm
UV}}+\log\frac{4\pi\mu^{2}}{s}-\gamma+2\,,$

where

$\displaystyle\omega_{1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}{s+m_{b}^{2}-m_{c}^{2}}\,,$
$\displaystyle\omega_{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}{s+m_{c}^{2}-m_{b}^{2}}\,,$
$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{m_{b}+m_{c}}{m_{b}-m_{c}}\,,$ $\displaystyle{\rm
Li_{2}}(x)$ $\displaystyle=$
$\displaystyle-\int_{0}^{x}dt\frac{\log(1-t)}{t}\,.$ (32)

Secondly, we introduce the notation

$\displaystyle\int dps$ $\displaystyle=$
$\displaystyle\int\frac{d^{D-1}\vec{k}}{2E_{k}}\frac{d^{D-1}\vec{p}_{1}}{2E_{p_{1}}}\frac{d^{D-1}\vec{p}_{2}}{2E_{p_{2}}}\delta^{D}(p-k-
p_{1}-p_{2})\,,$

for simplicity, and write down the analytical expressions of the three-body
phase-space integrals,

$\displaystyle R_{11}(s)$ $\displaystyle=$
$\displaystyle\frac{sm_{b}^{2}}{\pi^{2}\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}(2\pi)^{-4\varepsilon_{\rm
IR}}\mu^{-2\varepsilon_{\rm IR}}\int dps\frac{1}{(k\cdot p_{1})^{2}}$
$\displaystyle=$ $\displaystyle\frac{1}{2\varepsilon_{\rm IR}}-\log
4\pi+\gamma-1+\log\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}^{3}}{m_{b}m_{c}s\mu^{2}}-\frac{s+m_{b}^{2}-m_{c}^{2}}{2\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega_{1}}{1-\omega_{1}}\right)$
$\displaystyle-\frac{m_{b}^{2}-m_{c}^{2}}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega_{1}}{1-\omega_{1}}\right)-\frac{s-m_{b}^{2}+m_{c}^{2}}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega}{1-\omega}\right)$
$\displaystyle=$ $\displaystyle\overline{R}_{11}(s)+\frac{1}{2\varepsilon_{\rm
IR}}-\log
4\pi+\gamma-1+\log\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}^{3}}{m_{b}m_{c}s\mu^{2}}\,,$
$\displaystyle R_{22}(s)$ $\displaystyle=$ $\displaystyle
R_{11}(s)|_{m_{b}\leftrightarrow m_{c}}\,,$ $\displaystyle R_{12}(s)$
$\displaystyle=$
$\displaystyle\frac{s}{\pi^{2}\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}(2\pi)^{-4\varepsilon_{\rm
IR}}\mu^{-2\varepsilon_{\rm IR}}\int dps\frac{1}{k\cdot p_{1}k\cdot p_{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\left\\{\log\left(\frac{1+\omega}{1-\omega}\right)\left[\frac{1}{\varepsilon_{\rm
IR}}-2\log
4\pi+2\gamma-2+2\log\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}^{3}}{m_{b}m_{c}s\mu^{2}}\right]\right.$
$\displaystyle-2\log\frac{m_{b}}{m_{c}}\log\left(\frac{M+\omega}{M-\omega}\right)-\log^{2}\left(\frac{1+\omega}{1-\omega}\right)+2\log\frac{s}{\bar{s}}\log\left(\frac{1+\omega}{1-\omega}\right)-4{\rm
Li_{2}}\left(\frac{2\omega}{1+\omega}\right)$ $\displaystyle+2{\rm
Li_{2}}\left(\frac{\omega-1}{\omega-M}\right)+2{\rm
Li_{2}}\left(\frac{\omega-1}{\omega+M}\right)-2{\rm
Li_{2}}\left(\frac{\omega+1}{\omega-M}\right)-2{\rm
Li_{2}}\left(\frac{\omega+1}{\omega+M}\right)-\frac{1}{2}{\rm
Li_{2}}\left(\frac{1+\omega_{1}}{2}\right)$
$\displaystyle\left.-\frac{1}{2}{\rm
Li_{2}}\left(\frac{1+\omega_{2}}{2}\right)-{\rm
Li_{2}}\left(\omega_{1}\right)-{\rm Li_{2}}\left(\omega_{2}\right)+\frac{\log
2\log\left[(1+\omega_{1})(1+\omega_{2})\right]}{2}-\frac{\log^{2}2}{2}+\frac{\pi^{2}}{12}\right\\}\,,$
$\displaystyle=$
$\displaystyle\overline{R}_{12}(s)+\frac{1}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\log\left(\frac{1+\omega}{1-\omega}\right)\left[\frac{1}{\varepsilon_{\rm
IR}}-2\log
4\pi+2\gamma-2+2\log\frac{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}^{3}}{m_{b}m_{c}s\mu^{2}}\right]\,,$
$\displaystyle R^{1}_{12}(s)$ $\displaystyle=$
$\displaystyle\frac{s}{\pi^{2}\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\int
dps\frac{s-K^{2}}{k\cdot p_{1}k\cdot p_{2}}$ $\displaystyle=$
$\displaystyle\frac{s}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\left\\{\log^{2}(1-\omega)-\log^{2}(1+\omega)+2\log\frac{2s}{\bar{s}}\log\left(\frac{1+\omega}{1-\omega}\right)+2{\rm
Li_{2}}\left(\frac{1-\omega}{2}\right)\right.$ $\displaystyle\left.-2{\rm
Li_{2}}\left(\frac{1+\omega}{2}\right)+2{\rm
Li_{2}}\left(\frac{1+\omega}{1+M}\right)+2{\rm
Li_{2}}\left(\frac{1+\omega}{1-M}\right)-2{\rm
Li_{2}}\left(\frac{1-\omega}{1-M}\right)-2{\rm
Li_{2}}\left(\frac{1-\omega}{1+M}\right)\right\\}\,,$ $\displaystyle
R^{2}_{12}(s)$ $\displaystyle=$
$\displaystyle\frac{s}{\pi^{2}\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\int
dps\frac{\left(s-K^{2}\right)^{2}}{k\cdot p_{1}k\cdot p_{2}}$ (33)
$\displaystyle=$
$\displaystyle\frac{s^{2}}{\sqrt{\lambda(s,m_{b}^{2},m_{c}^{2})}}\left\\{\log^{2}(1-\omega)-\log^{2}(1+\omega)+2\log\frac{4s}{\bar{s}}\log\left(\frac{1+\omega}{1-\omega}\right)+2{\rm
Li_{2}}\left(\frac{1-\omega}{2}\right)\right.$ $\displaystyle-2{\rm
Li_{2}}\left(\frac{1+\omega}{2}\right)+2{\rm
Li_{2}}\left(\frac{1+\omega}{1+M}\right)+2{\rm
Li_{2}}\left(\frac{1+\omega}{1-M}\right)-2{\rm
Li_{2}}\left(\frac{1-\omega}{1-M}\right)-2{\rm
Li_{2}}\left(\frac{1-\omega}{1+M}\right)\,$
$\displaystyle\left.+\frac{2\omega\bar{s}}{s}-\frac{\bar{s}}{s}(1+\omega^{2})\log\left(\frac{1+\omega}{1-\omega}\right)\right\\}\,,$

where $\bar{s}=s-(m_{b}-m_{c})^{2}$.

## References

  * [1] F. Abe et al, Phys. Rev. D58 (1998) 112004.
  * [2] F. Abe et al, Phys. Rev. Lett. 81 (1998) 2432.
  * [3] T. Aaltonen et al, Phys. Rev. Lett. 100 (2008) 182002.
  * [4] V. M. Abazov et al, Phys. Rev. Lett. 101 (2008) 012001.
  * [5] R. L. Workman et al, Prog. Theor. Exp. Phys. 2022 (2022) 083C01.
  * [6] G. Aad et al, Phys. Rev. Lett. 113 (2014) 212004.
  * [7] A. M. Sirunyan et al, Phys. Rev. Lett. 122 (2019) 132001.
  * [8] R. Aaij et al, Phys. Rev. Lett. 122 (2019) 232001.
  * [9] Z. G. Wang, Eur. Phys. J. C73 (2013) 2559.
  * [10] S. Godfrey and N. Isgur, Phys. Rev. D32 (1985) 189.
  * [11] S. Godfrey, Phys. Rev. D70 (2004) 054017.
  * [12] D. Ebert, R. N. Faustov and V. O. Galkin, Phys. Rev. D67 (2003) 014027.
  * [13] S. N. Gupta and J. M. Johnson, Phys. Rev. D53 (1996) 312.
  * [14] J. Zeng, J. W. Van Orden and W. Roberts, Phys. Rev. D52 (1995) 5229.
  * [15] S. S. Gershtein, V. V. Kiselev, A. K. Likhoded and A. V. Tkabladze, Phys. Rev. D51 (1995) 3613.
  * [16] S. S. Gershtein, V. V. Kiselev, A. K. Likhoded and A. V. Tkabladze, Phys. Usp. 38 (1995) 1.
  * [17] E. J. Eichten and C. Quigg, Phys. Rev. D49 (1994) 5845.
  * [18] E. J. Eichten and C. Quigg, Phys. Rev. D99 (2019) 054025.
  * [19] A. P. Monteiro, M. Bhat and K. B. V. Kumar, Int. J. Mod. Phys. A32 (2017) 1750021.
  * [20] L. P. Fulcher, Phys. Rev. D60 (1999) 074006.
  * [21] N. R. Soni, B. R. Joshi, R. P. Shah, H. R. Chauhan and J. N. Pandya, Eur. Phys. J. C78 (2018) 592.
  * [22] P. G. Ortega, J. Segovia, D. R. Entem and F. Fernandez, Eur. Phys. J. C80 (2020) 223.
  * [23] Q. Li, M. S. Liu, L. S. Lu, Q. F. Lu, L. C. Gui and X. H. Zhong, Phys. Rev. D99 (2019) 096020.
  * [24] S. M. Ikhdair and R. Sever, Int. J. Mod. Phys. A19 (2004) 1771.
  * [25] S. M. Ikhdair and R. Sever, Int. J. Mod. Phys. A21 (2006) 6699.
  * [26] N. Brambilla and A. Vairo, Phys. Rev. D62 (2000) 094019.
  * [27] A. A. Penin, A. Pineda, V. A. Smirnov and M. Steinhauser, Phys. Lett. B593 (2004) 124.
  * [28] E. B. Gregory et al, Phys. Rev. Lett. 104 (2010) 022001.
  * [29] C. T. H. Davies et al, Phys. Lett. B382 (1996) 131.
  * [30] N. Mathur, M. Padmanath and S. Mondal, Phys. Rev. Lett. 121 (2018) 202002.
  * [31] B. D. Jones and R. M. Woloshyn, Phys. Rev. D60 (1999) 014502.
  * [32] A. A. EI-Hady, M. A. K. Lodhi and J. P. Vary, Phys. Rev. D59 (1999) 094001.
  * [33] G. L. Wang, Phys. Lett. B650 (2007) 15.
  * [34] G. L. Wang, T. Wang, Q. Li and C. H. Chang, JHEP 05 (2022) 006.
  * [35] E. Bagan, H. G. Dosch, P. Gosdzinsky, S. Narison and J. M. Richard, Z. Phys. C64 (1994) 57.
  * [36] M. Chabab, Phys. Lett. B325 (1994) 205.
  * [37] P. Colangelo, G. Nardulli and N. Paver, Z. Phys. C57 (1993) 43.
  * [38] V. V. Kiselev and A. V. Tkabladze, Phys. Rev. D48 (1993) 5208.
  * [39] Z. G. Wang, Acta Phys. Polon. B44 (2013) 1971.
  * [40] Z. G. Wang, Eur. Phys. J. A49 (2013) 131.
  * [41] T. M. Aliev, T. Barakat and S. Bilmis, Nucl. Phys. B947 (2019) 114726.
  * [42] M. Blanke et al, Phys. Rev. D99 (2019) 075006.
  * [43] S. Narison, Phys. Lett. B802 (2020) 135221.
  * [44] M. J. Baker, J. Bordes, C. A. Dominguez, J. Penarrocha and K. Schilcher, JHEP 07 (2014) 032.
  * [45] A. I. Onishchenko and O. L. Veretin, Eur. Phys. J. C50 (2007) 801.
  * [46] J. Lee, W. L. Sang and S. Kim, JHEP 01 (2011) 113.
  * [47] W. Tao and Z. J. Xiao, arXiv: 2310.17500 [hep-ph].
  * [48] W. Tao, R. L. Zhu and Z. J. Xiao, Phys. Rev. D106 (2022) 114037.
  * [49] W. Tao, R. L. Zhu and Z. J. Xiao, Eur. Phys. J. C83 (2023) 294.
  * [50] L. B. Chen and C. F. Qiao, Phys. Lett. B748 (2015) 443.
  * [51] W. L. Sang, H. F. Zhang and M. Z. Zhou, Phys. Lett. B839 (2023) 137812.
  * [52] F. Feng, Y. Jia, Z. Mo, J. Pan, W. L. Sang and J. Y. Zhang, arXiv: 2208.04302 [hep-ph].
  * [53] R. C. Verma, J. Phys. G39 (2012) 025005.
  * [54] S. Tang, Y. Li, P. Maris and J. P. Vary, Phys. Rev. D98 (2018) 114038.
  * [55] A. M. Badalian, B. L. G. Bakker and Yu. A. Simonov, Phys. Rev. D75 (2007) 116001.
  * [56] L. J. Reinders, H. Rubinstein and S. Yazaki, Phys. Rept. 127 (1985) 1.
  * [57] T. D. Lee and M. Nauenberg, Phys. Rev. 133 (1964) B1549.
  * [58] S. Narison, Phys. Lett. B706 (2012) 412.
  * [59] Z. G. Wang, Eur. Phys. J. C75 (2015) 427.