# Measurement of the 2$\nu\beta\beta$ Decay Half-Life of Se-82 with the Global
CUPID-0 Background Model

O. Azzolini  INFN Laboratori Nazionali di Legnaro, I-35020 Legnaro (Pd) -
Italy  J. W. Beeman  Lawrence Berkeley National Laboratory , Berkeley,
California 94720, USA  F. Bellini  Dipartimento di Fisica, Sapienza Università
di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  INFN, Sezione di Roma, P.le
Aldo Moro 2, 00185, Roma, Italy  M. Beretta  Present address: Department of
Physics, University of California, Berkeley, CA 94720, USA  Dipartimento di
Fisica, Università di Milano - Bicocca, I-20126 Milano - Italy  INFN Sezione
di Milano - Bicocca, I-20126 Milano - Italy  M. Biassoni  INFN Sezione di
Milano - Bicocca, I-20126 Milano - Italy  C. Brofferio  Dipartimento di
Fisica, Università di Milano - Bicocca, I-20126 Milano - Italy  INFN Sezione
di Milano - Bicocca, I-20126 Milano - Italy  C. Bucci  INFN Laboratori
Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy  S. Capelli
Dipartimento di Fisica, Università di Milano - Bicocca, I-20126 Milano - Italy
INFN Sezione di Milano - Bicocca, I-20126 Milano - Italy  V. Caracciolo
Present address: Dipartimento di Fisica, Università di Roma Tor Vergata,
I-00133, Rome, Italy  INFN Laboratori Nazionali del Gran Sasso, I-67100
Assergi (AQ) - Italy  L. Cardani  INFN, Sezione di Roma, P.le Aldo Moro 2,
00185, Roma, Italy  P. Carniti  Dipartimento di Fisica, Università di Milano -
Bicocca, I-20126 Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126
Milano - Italy  N. Casali  INFN, Sezione di Roma, P.le Aldo Moro 2, 00185,
Roma, Italy  E. Celi<EMAIL_ADDRESS>Gran Sasso Science Institute,
67100, L’Aquila - Italy  INFN Laboratori Nazionali del Gran Sasso, I-67100
Assergi (AQ) - Italy  D. Chiesa  Dipartimento di Fisica, Università di Milano
- Bicocca, I-20126 Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126
Milano - Italy  M. Clemenza  Dipartimento di Fisica, Università di Milano -
Bicocca, I-20126 Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126
Milano - Italy  I. Colantoni  INFN, Sezione di Roma, P.le Aldo Moro 2, 00185,
Roma, Italy  Consiglio Nazionale delle Ricerche, Istituto di Nanotecnologia,
c/o Dip. Fisica, Sapienza Università di Roma, 00185, Rome, Italy  O. Cremonesi
INFN Sezione di Milano - Bicocca, I-20126 Milano - Italy  A. Cruciani  INFN,
Sezione di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  A. D’Addabbo  INFN
Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy  I. Dafinei
Gran Sasso Science Institute, 67100, L’Aquila - Italy  INFN, Sezione di Roma,
P.le Aldo Moro 2, 00185, Roma, Italy  S. Di Domizio  Dipartimento di Fisica,
Università di Genova, I-16146 Genova - Italy  INFN Sezione di Genova, I-16146
Genova - Italy  V. Dompè  Dipartimento di Fisica, Sapienza Università di Roma,
P.le Aldo Moro 2, 00185, Roma, Italy  INFN, Sezione di Roma, P.le Aldo Moro 2,
00185, Roma, Italy  G. Fantini  Dipartimento di Fisica, Sapienza Università di
Roma, P.le Aldo Moro 2, 00185, Roma, Italy  INFN, Sezione di Roma, P.le Aldo
Moro 2, 00185, Roma, Italy  F. Ferroni  INFN, Sezione di Roma, P.le Aldo Moro
2, 00185, Roma, Italy  Gran Sasso Science Institute, 67100, L’Aquila - Italy
L. Gironi  Dipartimento di Fisica, Università di Milano - Bicocca, I-20126
Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126 Milano - Italy  A.
Giuliani  CNRS/CSNSM, Centre de Sciences Nucl$\acute{e}$aires et de Sciences
de la Mati$\grave{e}$re, 91405 Orsay, France  P. Gorla  INFN Laboratori
Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy  C. Gotti  INFN Sezione
di Milano - Bicocca, I-20126 Milano - Italy  G. Keppel  INFN Laboratori
Nazionali di Legnaro, I-35020 Legnaro (Pd) - Italy  J. Kotila  University of
Jyväskylä, Department of Physics, P. O. Box 35 (YFL), FI-40014, Finland
Finnish Institute for Educational Research, P.O.Box 35 FI-40014 University of
Jyväskylä - Finland  Center for Theoretical Physics, Sloane Physics
Laboratory, Yale University, New Haven, Connecticut 06520-8120 - USA  M.
Martinez  Present address: Centro de Astropartículas y Física de Altas
Energías, Universidad de Zaragoza, and ARAID, Fundación Agencia Aragonesa para
la Investigación y el Desarrollo, Gobierno de Aragón, Zaragoza 50018, Spain
Dipartimento di Fisica, Sapienza Università di Roma, P.le Aldo Moro 2, 00185,
Roma, Italy  INFN, Sezione di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  S.
S. Nagorny  Present address: Department of Physics $\&$ Engineering Physics
Astronomy, Queen’s University Kingston, Ontario, K7L 3N6 Kingston, Canada
INFN Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy  M.
Nastasi  Dipartimento di Fisica, Università di Milano - Bicocca, I-20126
Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126 Milano - Italy  S.
Nisi  INFN Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy
C. Nones  IRFU, CEA, Universit$\acute{e}$ Paris-Saclay, F-91191 Gif-sur-
Yvette, France  D. Orlandi  INFN Laboratori Nazionali del Gran Sasso, I-67100
Assergi (AQ) - Italy  L. Pagnanini<EMAIL_ADDRESS>Gran Sasso
Science Institute, 67100, L’Aquila - Italy  INFN Laboratori Nazionali del Gran
Sasso, I-67100 Assergi (AQ) - Italy  M. Pallavicini  Dipartimento di Fisica,
Università di Genova, I-16146 Genova - Italy  INFN Sezione di Genova, I-16146
Genova - Italy  L. Pattavina  Present address: Physik-Department and
Excellence Cluster Origins, Technische Universität München, 85747 Garching,
Germany  INFN Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ) -
Italy  M. Pavan  Dipartimento di Fisica, Università di Milano - Bicocca,
I-20126 Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126 Milano -
Italy  G. Pessina  INFN Sezione di Milano - Bicocca, I-20126 Milano - Italy
V. Pettinacci  INFN, Sezione di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  S.
Pirro  INFN Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy
S. Pozzi  Dipartimento di Fisica, Università di Milano - Bicocca, I-20126
Milano - Italy  INFN Sezione di Milano - Bicocca, I-20126 Milano - Italy  E.
Previtali  Dipartimento di Fisica, Università di Milano - Bicocca, I-20126
Milano - Italy  INFN Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ)
- Italy  A. Puiu  INFN Laboratori Nazionali del Gran Sasso, I-67100 Assergi
(AQ) - Italy  A. Ressa  Dipartimento di Fisica, Sapienza Università di Roma,
P.le Aldo Moro 2, 00185, Roma, Italy  INFN, Sezione di Roma, P.le Aldo Moro 2,
00185, Roma, Italy  C. Rusconi  INFN Laboratori Nazionali del Gran Sasso,
I-67100 Assergi (AQ) - Italy  Department of Physics and Astronomy, University
of South Carolina, Columbia, SC 29208 - USA  K. Schäffner  Present address:
Max-Planck-Institut für Physik, 80805 München - Germany  INFN Laboratori
Nazionali del Gran Sasso, I-67100 Assergi (AQ) - Italy  C. Tomei  INFN,
Sezione di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  M. Vignati
Dipartimento di Fisica, Sapienza Università di Roma, P.le Aldo Moro 2, 00185,
Roma, Italy  INFN, Sezione di Roma, P.le Aldo Moro 2, 00185, Roma, Italy  A.
S. Zolotarova  IRFU, CEA, Universit$\acute{e}$ Paris-Saclay, F-91191 Gif-sur-
Yvette, France

###### Abstract

We report on the results obtained with the global CUPID-0 background model,
which combines the data collected in the two measurement campaigns for a total
exposure of 16.6 kg$\times$yr. We identify with improved precision the
background sources in the region of interest for neutrinoless double
$\beta$-decay, making more solid the foundations for the background budget of
the next-generation CUPID experiment. Relying on the excellent data
reconstruction, we measure the two-neutrino double $\beta$-decay half-life of
82Se with unprecedented accuracy: $T_{1/2}^{2\nu}=[8.69\pm
0.05\textrm{(stat.)}~{}^{+0.09}_{-0.06}\textrm{(syst.)}]\times
10^{19}~{}\textrm{yr}$.

scintillating bolometers, background model, two-neutrino double $\beta$-decay,
spectral shape

###### pacs:

07.20.Mc, 23.40.-s, 21.10.Tg, 27.50.+e

Introduction. The comprehension of the neutrino nature is of fundamental
importance to understand the origin of neutrino masses and to have hints of
new physics beyond the Standard Model (SM). Thus far, the neutrinoless double
$\beta$-decay ($0\nu\beta\beta$) Furry (1939) is the only practical way to
probe if neutrinos are Dirac ($\nu\neq\bar{\nu}$) particles, as expected by
the SM, or they are described by Majorana’s theory ($\nu=\bar{\nu}$) Dell’Oro
_et al._ (2016). The double $\beta$-decay is a weak process involving the
decay of two neutrons into protons with the emission of two electrons.
According to the SM, such a process is allowed only with the emission of two
anti-neutrinos in the final state Goeppert-Mayer (1935), the so-called two-
neutrino double $\beta$-decay (2$\nu\beta\beta$). In the past decades, a big
experimental effort was committed to the search for $0\nu\beta\beta$ allowing
the existing detectors to reach sensitivities of the order of
$10^{25}-10^{26}$ yr on its half-life Agostini  _et al._ (2020); Abe  _et al._
(2023); Adams  _et al._ (2022). Nowadays, none of them has ever observed such
decay, however, the development of technologies able to reach excellent
sensitivities led to outstanding precision measurements of 2$\nu\beta\beta$
Saakyan (2013); Barabash (2015); Gando  _et al._ (2019); Arnold  _et al._
(2019). The properties of cryogenic calorimeters make them competitive
detectors both in the search for $0\nu\beta\beta$ and in the measurement of
2$\nu\beta\beta$. Being the signal source embedded into the detector, this
technology offers high detection efficiency and excellent energy resolution.
Furthermore, using different crystals opens the opportunity to study several
double-$\beta$ decay emitters. In the past years, the CUORE (Cryogenic
Underground Observatory for Rare Events) experiment demonstrated the
possibility to reach ton-scale exposures Adams  _et al._ (2022), while several
R&D projects Pirro  _et al._ (2006); Arnaboldi  _et al._ (2011); Beeman  _et
al._ (2013a); Artusa  _et al._ (2016); Armengaud  _et al._ (2017) developed a
powerful technology based on scintillating bolometers to enable particle
discrimination, thus reducing the radioactive background. These lay the
foundations for CUPID (CUORE Upgrade with Particle IDentification), a next-
generation experiment aiming to reach unprecedented sensitivities using a
large-size array of scintillating cryogenic calorimeters Armstrong  _et al._
(2019). In this letter, we present the results of the global background model
of CUPID-0, the first demonstrator of the CUPID technology, and a novel
measurement of the 82Se-2$\nu\beta\beta$ half-life, which results to be the
most precise so far.

Figure 1: Single hit spectra registered by CUPID-0 in phase-I (gray) and
phase-II (red). Below 3 MeV, we identify the dominant 2$\nu\beta\beta$
contribution (green), and the $\gamma$-rays peaks due to cosmogenic activation
of Zn82Se (65Zn), and setup contaminants
(${}^{40}\mathrm{K}$,${}^{208}\mathrm{Tl}$). At 3-4 MeV degraded
$\alpha$-particles produce a flat background, while the $\alpha$-peaks above 4
MeV are the distinctive signatures of the different detector contaminants. The
main differences between phase-I and phase-II traceable in the data are the
65Zn reduction due to its short half-life, the ${}^{40}\mathrm{K}$ and
${}^{208}\mathrm{Tl}$ reduction due to the 10 mK (1 cm of copper) shield,
which also introduces an additional ${}^{210}\mathrm{Po}$ component in the 3-4
MeV region.

CUPID-0 detector. CUPID-0 is composed of 24 ZnSe scintillating crystals
Dafinei  _et al._ (2017) enriched at (95.3 $\pm$ 0.3)% in 82Se Balabanov  _et
al._ (2023). The array is completed by two natural ZnSe crystals, for a total
ZnSe mass of 10.5 kg. The scintillation light is detected by thin Germanium
wafers operating as cryogenic calorimeters Beeman  _et al._ (2013b). A
Vikuiti${}^{\text{TM}}$ plastic reflective foil surrounds each ZnSe crystal to
enhance the light collection. ZnSe crystals and light detectors are
alternatively staked in 5 towers, and supported by a copper structure to which
they are connected by means of PTFE clamps. To convert the temperature
variation produced by interacting particles into a voltage signal, a germanium
thermistor is glued to each device (both light detectors and ZnSe) and bonded
to the cryogenic readout system through thin gold wires. The detector is
installed in an Oxford 1000 3He/4He dilution refrigerator working at a base
temperature of 12 mK. The whole experimental setup is located underground in
the Laboratori Nazionali del Gran Sasso (LNGS), in Italy. More details can be
found in CUPID-0 detector paper Azzolini  _et al._ (2018a). The entire data-
taking amounts to nearly three years, from June 2017 to February 2020, with a
few months off from January to June 2019 to upgrade the detector. We refer to
the two physics runs of the experiment as phase-I and phase-II, respectively.
In the phase-II, the reflective foils are removed and we add a further copper
shield at the 10 mK stage directly facing the detector. In phase-I, the
reflectors between the ZnSe crystals affect the correct identification of
surface/bulk $\alpha$-decays due to detector contaminants Azzolini  _et al._
(2021). Moreover, from the phase-I background model, we also noticed of a
possible background affecting the $0\nu\beta\beta$ region of interest from the
elements near the detector. Excluding the reflectors from the possible
background sources is fundamental to assess if they can be used in next-
generation experiments Armatol  _et al._ (2021). Without the reflecting foils,
the $\alpha$ discrimination power is still excellent Azzolini  _et al._
(2022), and we can disentangle contaminants coming from crystal surfaces and
components close to the detector. Moreover, the presence of the shield at 10
mK allows us to better study the cryostat contaminants, exploiting the
different counting rates of the $\gamma$-ray peaks in phase-I and phase-II
(see Fig.1).

Data taking and analysis. Given an energy release inside the ZnSe, the
relative voltage signal is triggered and saved in a 5 s time window. The
output of the corresponding light detectors is recorded by analyzing the
coincidence in time with the ZnSe signal. Light pulses are acquired in a 500
ms time window Di Domizio  _et al._ (2018). We apply a matched-filter
algorithm Gatti and Manfredi (1986) on the raw data to estimate the pulse
height and the shape parameters. Then, we calibrate each detector by analyzing
the data acquired with an external 232Th source. The 232Th decay chain
provides several $\gamma$-lines in the range (511-2615) keV that can be used
as reference peaks by fitting with a zero intercept parabolic function the
non-calibrated amplitudes. A dedicated 56Co calibration with $\gamma$-lines
above 2615 keV is also performed to investigate the uncertainty on the energy
scale with the distribution of residuals Azzolini  _et al._ (2019a). Finally,
we tag the $\alpha$ events by looking at the pulse shape parameters of the
light signal pulses Azzolini  _et al._ (2018b).

Different sources can be disentangled thanks to the following features: (i)
the identification of $\alpha$ particles from $\beta$/$\gamma$ down to 2 MeV;
(ii) the excellent resolution, enabling the study of the distortion in the
shape of $\alpha$ peaks, essential to tag the surface contaminants; (iii) the
granularity of the detector, exploited to distinguish interactions in which
one or multiple crystals are involved. For this purpose, we tag the events
triggering multiple crystals within a $\pm 20$ ms time window.

In this analysis, we refer to the $\mathcal{M}_{1\beta/\gamma}$ and
$\mathcal{M}_{1\alpha}$ spectra as the single-hits events of $\beta$/$\gamma$
and $\alpha$ respectively, while to $\mathcal{M}_{2}$ and $\Sigma_{2}$ spectra
as double-hit events ($\mathcal{M}_{2}$ comprises the energies detected by
each crystal, $\Sigma_{2}$ the total energy released in two crystals). Data
selection is based on several quality cuts to reject non-physical events and
multiple pulses in the same acquisition window. The relative efficiency of
these cuts is evaluated for phase-I and phase-II separately and combined with
the trigger and energy reconstruction efficiencies Azzolini  _et al._ (2019b).
The result is $\varepsilon_{1}=(95.7\pm 0.5)$ and $\varepsilon_{2}=(94.8\pm
0.7)$ for phase-I and phase-II, respectively.

Data Modeling. Initially, we reproduce the CUPID-0 geometry with the Geant4
toolkit. Following the same grouping criteria used in the phase-I background
model, we split the detector components in Crystals, Holder, Reflectors
(phase-I), and 10mK (phase-II). The light detectors are included in the
geometry but not considered as a background source, being their radioactive
contaminants and masses very small Azzolini  _et al._ (2019b). The groups
named Holder, Reflectors, and 10mK include also the elements close to the
detector (e.g. wiring, glue), whose contaminants will not be generated. The
different components of shielding and cryostat are also reproduced in detail.
However, for the purpose of this paper, we report their overall contribution
to the background, labeled as Cryostat (details in Supplemental Materials). We
separately generate for phase-I and phase-II the radioactive decays from
various background sources in each group of volumes. In particular, we
simulate the initial kinematics of the 2$\nu\beta\beta$ assuming the Single
State Dominance (SSD) mechanism, calculating the electron energy using exact
Dirac electron wave functions that take into account finite nuclear size and
electron screening Kotila and Iachello (2012), while the generation of
external muons at LNGS described in detail in Andreotti  _et al._ (2010). Each
simulation is then processed with a custom software to reproduce both detector
and data processing features, like time and energy resolution, energy
threshold, $\alpha$-particle identification, and time coincidence window
Azzolini  _et al._ (2018b). The identification of the radioactive sources and
their location is driven by material screening results and CUPID-0 data
analysis. An $\alpha$-decay occurring in crystal bulk produces a peak at the
Q-value of the decay in $\mathcal{M}_{1\alpha}$, while contaminants on the
crystal surface can be distinguished in $\mathcal{M}_{2}$ and $\Sigma_{2}$. In
phase-I, the presence of reflecting foils suppress the signature of surface
contaminations due to the absorption of $\alpha$-particles or recoil escapes.
However, in this analysis, we can recover that information by combining it
with the phase-II data. ${}^{232}\mathrm{Th}$ and ${}^{238}\mathrm{U}$ decay
chains in the close components can produce sharp peaks at the nominal
$\alpha$-particle energies if the contaminations are very shallow, or a
continuum spectrum if the decays occur deeply from the surface. To model the
surface contaminants both for crystals and near elements (Reflectors, Holder,
10mK), we parameterize its density profile as $e^{-x/\lambda}$, where
$\lambda$ is the depth parameter. Based on previous experiences Azzolini  _et
al._ (2019b), we used $\lambda=10$ nm to model shallow contamination and
$\lambda=10$ $\mu$m for deeper ones. Besides this, the $\gamma$-ray lines in
$\mathcal{M}_{1\beta/\gamma}$, $\mathcal{M}_{2}$, and $\Sigma_{2}$ are
typically the signature of the external sources (e.g. Cryostat and shields).
In this regard, comparing the spectra of phase-II (additional $1$ cm thick
10mK shield) with those of phase-I, allow us to improve the knowledge of the
cryostat contaminants exploiting the relative reduction of the $\gamma$-ray
peaks. For these sources, we also exploited the information from the
background model of the CUORE-0 experiment that was hosted in the same
facility Alduino  _et al._ (2017); Azzolini  _et al._ (2019b). In summary, we
included as bulk contamination the natural decay chains due to
${}^{232}\mathrm{Th}$, ${}^{238}\mathrm{U}$, ${}^{235}\mathrm{U}$, and
${}^{40}\mathrm{K}$ for all the simulated elements, the cosmogenic activation
products of copper (54Mn, 58Co, 60Co) and ZnSe (65Zn, 60Co), and the
82Se-2$\nu\beta\beta$ decay. Moreover, we included the ${}^{232}\mathrm{Th}$
and ${}^{238}\mathrm{U}$ contaminations on the surfaces of ZnSe crystals and
near elements. Being the environmental $\gamma$-rays and the neutrons reduced
to a negligible level by dedicated shielding Alduino  _et al._ (2017), we
consider only muons as an external source.

We identify the breakpoints of the secular equilibrium in the radioactive
decay chains at 228Ra for the ${}^{232}\mathrm{Th}$ chain, and at 226Ra/210Pb
for the ${}^{238}\mathrm{U}$ one, as expected from the different chemical
behavior. We select the set of cosmogenically activated contaminants by
looking at the presence of their characteristic $\gamma$-peaks and their half-
lives. Finally, the muon contribution is normalized on the number of shower
events that simultaneously trigger more than three crystals. The result of
such normalization is compatible with the one obtainable using the muon flux
at LNGS Ambrosio  _et al._ (1995).

Figure 2: Experimental background indexes (black, top pad), compared to the
reconstructed ones for phase-I (blue) and phase-II (orange). The total
background index is the sum of all components with their correlations. The
error bars represent the systematic uncertainty of the background
reconstruction.

Spectral Fit. Exploiting a multivariate Bayesian fit, we reproduce the
experimental spectra with a linear combination of the simulated spectra of the
different background sources. As a result, we obtain a reliable estimation of
the activity of each source considered.

The fit is performed with the JAGS software Plummer (2003), which exploits the
Monte Carlo Markov Chain sampling to sample the joint posterior probability
density function of the model unknowns, which linearly depend on the source
activities. We performed a combined fit on 8 experimental spectra,
corresponding to $\mathcal{M}_{1\beta/\gamma}$, $\mathcal{M}_{1\alpha}$,
$\mathcal{M}_{2}$, and $\Sigma_{2}$ for both phase-I and phase-II. The fit
energy range extends from 300 keV to 5 MeV for $\mathcal{M}_{1\beta/\gamma}$,
and from 2 to 7 MeV for $\mathcal{M}_{1\alpha}$. For the $\mathcal{M}_{2}$ and
$\Sigma_{2}$, we require that both events deposit an energy higher than 150
keV, with the exception of events in time coincidence with $\alpha$-particles,
for which the threshold is 100 keV to include the nuclear recoils. We use a
non-uniform bin size to compact $\alpha$-peaks and $\gamma$-ray lines in a
single bin, avoiding the mismatch between simulations and real data induced by
the line shape of the peaks. Moreover, we grouped more bins together if the
total number of events is lower than 30 counts, while, in any other case, we
use a 15 keV fixed step binning. We use non-negative flat priors for all the
sources whose activity cannot be constrained a priori with independent
measurements. In particular, we use the results obtained in the CUORE-0
background model to set a prior on the ${}^{232}\mathrm{Th}$ and
${}^{238}\mathrm{U}$ in the Holder and the 60Co in the cryostat components.
Since the phase-I and phase-II data sets are consecutive in time, we can
reasonably assume that the activity of all the contaminants with a half-life
longer than the overall CUPID-0 lifetime is constant. For this reason, we
constrain a subset of simulations to have the same activity between phase-I
and phase-II.We have a total of 78 simulations, with 17 couples of them
constrained (see Supplemental Materials). We consider as Reference Fit the one
performed with the simulations, the energy ranges, and the binning described
above. The fit details, like the spectral reconstruction over the 8
experimental spectra and the pull distribution, are shown in the Supplemental
Materials. We reconstruct the background index in the Region Of Interest (ROI)
by selecting the events in a 400 keV interval centered at the
$\mathrm{Q}_{\beta\beta}$ of 82Se-2$\nu\beta\beta$.We also apply the delayed
coincidence cut on the ${}^{232}\mathrm{Th}$ simulations in order to reproduce
the selection criteria of the $0\nu\beta\beta$ analysis Azzolini  _et al._
(2022, 2019b). The background index reconstructed with the fit for each group
of sources in phase-I and phase-II is reported in Fig. 2.

Systematics. To evaluate the systematic uncertainty affecting the background
index, we repeat the fit multiple times in different conditions. We include in
this category different energy thresholds (400, 500, 700 keV) and fixed step
binnings (30, 50, and 100 keV) in $\mathcal{M}_{1\beta/\gamma}$. Moreover, we
run the fit using a reduced list of sources, removing those sources whose
contribution results at limit in the Reference fit. We evaluate also
systematic effects due to energy miscalibration, re-scaling the data for the
residuals calculated in Azzolini  _et al._ (2018b). We probe the uncertainty
about the contaminant position in cryostat and shields removing alternatively
one of the elements from the source list. Finally, we investigate an
alternative description of surface crystal contaminants and near elements. In
the former case, we substitute the simulations at 10 nm depth with the 1 nm
equivalent ones, while in the latter case, we substitute the
${}^{232}\mathrm{Th}$ and ${}^{238}\mathrm{U}$ contamination at 10 $\mu$m with
bulk ones in the Reflectors (phase-I).

We evaluate the systematic uncertainty on the 82Se half-life including
additional tests to the previous list. In particular, we perform supplementary
fits by switching off the $\alpha$ identification and using non-informative
flat priors for all the contributions. Ultimately, we consider also a possible
contamination due to 90Sr/90Y in the ZnSe crystals. 90Sr/90Y is an
anthropogenic isotope that decays purely through $\beta$-emission with a
lifetime of 28.8 y and produces a featureless spectrum that correlates with
the 2$\nu\beta\beta$ one in the fit. Since the presence of 90Sr/90Y in the
ZnSe crystals is unclear, we consider this possibility as systematic. Finally,
we account for possible theoretical uncertainties on the 2$\nu\beta\beta$
decay model generating the initial kinematic of the electrons with an
alternative spectrum calculated under the Closure Approximation Kotila and
Iachello (2012).

Background Index. Ultimately, after the investigation of systematic
uncertainties, we reproduce a background index of $[4.0\pm 0.2(stat.)\pm
0.4(syst.)]\times 10^{-3}$ counts/keV/kg/yr for phase-I and $[3.9\pm
0.2(stat.)_{-0.3}^{+0.4}(syst.)]\times 10^{-3}$ counts/keV/kg/yr for phase-II,
compatible among them and with the observed background of $(3.5\pm 1.0)\times
10^{-3}$ counts/keV/kg/yr and $(5.5\pm 1.5)\times 10^{-3}$ counts/keV/kg/yr
for phase-I and phase-II respectively. According to our model, we can
attribute the slightly higher background observed in phase-II to a statistical
fluctuation. Thanks to the combined analysis, it is recognized that the
background due to crystal contamination is 30% lower than the previously
published result, while the contribution of the cryostat is 50% higher.
Finally, we also clarify that the previous model overestimated by a factor of
two the near elements contamination, which turns out to be compatible with the
requirements of next-generation experiments.

Figure 3: Measurement of the 82Se-2$\nu\beta\beta$ half-life obtained in the
systematic tests (bottom) and as final result of this analysis, compared to
the previous result of CUPID-0 (top). The error bars shown in the bottom plot
correspond to the statistical error of each systematic test, while the red
vertical lines correspond to the Reference mean value (solid line) and its
statistical uncertainty (dashed line). The 90Sr/90Y strongly correlates with
the 2$\nu\beta\beta$, thus broadening the uncertainty. The statistical
uncertainty is then combined with the systematic one and compared with the
previous result (top).

82Se-2$\nu\beta\beta$ decay half-life. We fit the 2$\nu\beta\beta$ spectrum in
the background model constraining its activity to be the same in phase-I and
phase-II. The reconstructed 82Se-2$\nu\beta\beta$ activity from the background
model considering 4.65 kg of active mass is [$8.62\pm 0.04$ (stat.)] mBq.
Differently from the previous analysis Azzolini  _et al._ (2019c), we include
the efficiencies as fit parameters with a Gaussian prior set at their measured
value. The resulting marginalized posterior includes the efficiencies and 82Se
isotopic abundance uncertainties as nuisance parameters. Thus, these
uncertainties are included in the statistical one ($\pm 0.6\%$). The results
of the systematic studies are summarized in Fig. 3, where we report the
percentage difference for the 82Se-2$\nu\beta\beta$ half-life between each
test and the Reference fit. For the threshold and the source location, we
quote the maximum difference among those of the different tests. We assume
each fit systematics to be single-sided and estimate the uncertainty as 68% of
the difference between the test fit result and the Reference one. Thus, the
overall systematic uncertainty of the 82Se-2$\nu\beta\beta$ half-life is
(${}^{+1.0\%}_{-0.7\%}$), evaluated as the sum in quadrature of the single
uncertainties. In summary, the final value of the 82Se-2$\nu\beta\beta$ half-
life value measured by CUPID-0 is:

$T_{1/2}^{2\nu}=[8.69\pm
0.05\textrm{(stat.)}~{}^{+0.09}_{-0.06}\textrm{(syst.)}]\times
10^{19}~{}\textrm{yr}.$

Combining in quadrature statistical and systematical uncertainty, we can write
the half-life of 82Se-2$\nu\beta\beta$ as
$T_{1/2}^{2\nu}=[8.69~{}^{+0.10}_{-0.08}]\times 10^{19}~{}\textrm{yr}$,
compatible with the 2019 measurement of CUPID-0 Azzolini  _et al._
(2019c).Finally, we can convert the half-life in terms of nuclear matrix
element (NME) using the Phase Space Factor $G^{2\nu}=(1.996\pm 0.028)\times
10^{-18}$ y-1 obtained under the SSD model Kotila and Iachello (2012); Kotila
(2023). Without making hypotheses on the axial coupling constant $g_{A}$, we
can provide a new experimental value of the effective NME,
$\mathcal{M}^{eff}_{2\nu}=\mathcal{M}_{2\nu}g_{A}^{2}$ =
$0.0760~{}^{+~{}0.0006}_{-~{}0.0007}$, to be compared with theoretical
calculations. This measurement represents the new benchmark for nuclear
models, which currently obtain significantly higher values Nomura (2022);
Barea  _et al._ (2015); Popara  _et al._ (2022); Coraggio  _et al._ (2019);
Ejiri  _et al._ (2019) and resort to a phenomenological quenching of $g_{A}$
to explain the discrepancy.

In conclusion, by combining the experimental data collected by CUPID-0 in
phase-I and phase-II, we get a better comprehension of the background
affecting the region of interest of 82Se-$0\nu\beta\beta$. Thanks to the
excellent data reconstruction, we measure the 82Se-2$\nu\beta\beta$ half-life
with an accuracy of (${}^{+1.2\%}_{-0.9\%}$), obtaining one of the most
precise measurements of 2$\nu\beta\beta$.

###### Acknowledgements.

This work was partially supported by the European Research Council
(FP7/2007-2013), Contract No. 247115. The work of J. Kotila was supported by
the Academy of Finland (Grant Nos. 314733 and 345869). L. P. research
activities are supported by European Union’s Horizon 2020306 research and
innovation program under the Marie Skłodowska-Curie grant agreement N.
10102968. We are particularly grateful to S. Grigioni and M. Veronelli for
their help in the design and construction of the sensor-to-absorber gluing
system, M. Iannone for the help in all the stages of the detector
construction, A. Pelosi for the construction of the assembly line, M. Guetti
for the assistance in the cryogenic operations, R. Gaigher for the calibration
system mechanics, M. Lindozzi for the development of cryostat monitoring
system, M. Perego for his invaluable help, the mechanical workshop of LNGS (E.
Tatananni, A. Rotilio, A. Corsi, and B. Romualdi) for the continuous help in
the overall setup design. We acknowledge the Dark Side Collaboration for the
use of the low-radon clean room. This work makes use of the DIANA data
analysis and APOLLO data acquisition software which has been developed by the
CUORICINO, CUORE, LUCIFER, and CUPID-0 Collaborations. This work makes use of
the Arby software for Geant4-based Monte Carlo simulations, which has been
developed in the framework of the Milano - Bicocca R&D activities and that is
maintained by O. Cremonesi and S. Pozzi.

## References

  * Furry (1939) W. H. Furry, Phys. Rev. 56, 1184 (1939).
  * Dell’Oro  _et al._ (2016) S. Dell’Oro, S. Marcocci, M. Viel, and F. Vissani, Adv. High Energy Phys. 2016, 2162659 (2016).
  * Goeppert-Mayer (1935) M. Goeppert-Mayer, Phys. Rev. 48, 512 (1935).
  * Agostini  _et al._ (2020) M. Agostini _et al._ (GERDA), Phys. Rev. Lett. 125, 252502 (2020), arXiv:2009.06079 [nucl-ex] .
  * Abe  _et al._ (2023) S. Abe _et al._ (KamLAND-Zen), Phys. Rev. Lett. 130, 051801 (2023), arXiv:2203.02139 [hep-ex] .
  * Adams  _et al._ (2022) D. Q. Adams _et al._ (CUORE), Nature 604, 53 (2022), arXiv:2104.06906 [nucl-ex] .
  * Saakyan (2013) R. Saakyan, Ann. Rev. Nucl. Part. Sci. 63, 503 (2013).
  * Barabash (2015) A. S. Barabash, Nucl. Phys. A935, 52 (2015), arXiv:1501.05133 [nucl-ex] .
  * Gando  _et al._ (2019) A. Gando _et al._ (KamLAND-Zen), Phys. Rev. Lett. 122, 192501 (2019), arXiv:1901.03871 [hep-ex] .
  * Arnold  _et al._ (2019) R. Arnold _et al._ (NEMO-3), Eur. Phys. J. C79, 440 (2019), arXiv:1903.08084 [nucl-ex] .
  * Pirro  _et al._ (2006) S. Pirro _et al._ , Phys. Atom. Nucl. 69, 2109 (2006).
  * Arnaboldi  _et al._ (2011) C. Arnaboldi, S. Capelli, O. Cremonesi, L. Gironi, M. Pavan, G. Pessina, and S. Pirro, Astropart. Phys. 34, 344 (2011), arXiv:1006.2721 [nucl-ex] .
  * Beeman  _et al._ (2013a) J. W. Beeman _et al._ , JINST 8, P05021 (2013a).
  * Artusa  _et al._ (2016) D. R. Artusa _et al._ , Eur. Phys. J. C 76, 364 (2016).
  * Armengaud  _et al._ (2017) E. Armengaud _et al._ , Eur. Phys. J. C 77, 785 (2017).
  * Armstrong  _et al._ (2019) W. R. Armstrong _et al._ (CUPID), (2019), arXiv:1907.09376 [physics.ins-det] .
  * Dafinei  _et al._ (2017) I. Dafinei _et al._ , J. Cryst. Growth 475, 158 (2017), arXiv:1702.05877 [physics.ins-det] .
  * Balabanov  _et al._ (2023) S. S. Balabanov _et al._ , JINST 18, P04035 (2023).
  * Beeman  _et al._ (2013b) J. W. Beeman _et al._ , JINST 8, P07021 (2013b), arXiv:1304.6289 [physics.ins-det] .
  * Azzolini  _et al._ (2018a) O. Azzolini _et al._ (CUPID), Eur. Phys. J. C78, 428 (2018a), arXiv:1802.06562 [physics.ins-det] .
  * Azzolini  _et al._ (2021) O. Azzolini _et al._ , Eur. Phys. J. C 81, 722 (2021), arXiv:2105.04409 [physics.ins-det] .
  * Armatol  _et al._ (2021) A. Armatol _et al._ (CUPID), Eur. Phys. J. C 81, 104 (2021), arXiv:2011.13656 [physics.ins-det] .
  * Azzolini  _et al._ (2022) O. Azzolini _et al._ (CUPID), Phys. Rev. Lett. 129, 111801 (2022), arXiv:2206.05130 [nucl-ex] .
  * Di Domizio  _et al._ (2018) S. Di Domizio _et al._ , JINST 13, P12003 (2018), arXiv:1807.11446 [physics.ins-det] .
  * Gatti and Manfredi (1986) E. Gatti and P. F. Manfredi, Riv. Nuovo Cimento 9, 1 (1986).
  * Azzolini  _et al._ (2019a) O. Azzolini _et al._ (CUPID), Phys. Rev. Lett. 123, 032501 (2019a), arXiv:1906.05001 [nucl-ex] .
  * Azzolini  _et al._ (2018b) O. Azzolini _et al._ , Eur. Phys. J. C78, 734 (2018b), arXiv:1806.02826 [physics.ins-det] .
  * Azzolini  _et al._ (2019b) O. Azzolini _et al._ (CUPID), Eur. Phys. J. C79, 583 (2019b), arXiv:1904.10397 [nucl-ex] .
  * Kotila and Iachello (2012) J. Kotila and F. Iachello, Phys. Rev. C85, 034316 (2012), arXiv:1209.5722 [nucl-th] .
  * Andreotti  _et al._ (2010) E. Andreotti _et al._ (CUORICINO), Astropart. Phys. 34, 18 (2010), arXiv:0912.3779 [nucl-ex] .
  * Alduino  _et al._ (2017) C. Alduino _et al._ (CUORE), Eur. Phys. J. C77, 13 (2017), arXiv:1609.01666 [nucl-ex] .
  * Ambrosio  _et al._ (1995) M. Ambrosio _et al._ (MACRO), Phys. Rev. D52, 3793 (1995).
  * Plummer (2003) M. Plummer, JAGS: A Program for Analysis of Bayesian Graphical Models using Gibbs Sampling. 3rd International Workshop on Distributed Statistical Computing (DSC 2003); Vienna, Austria. 124 (2003).
  * Azzolini  _et al._ (2019c) O. Azzolini _et al._ , Phys. Rev. Lett. 123, 262501 (2019c), arXiv:1909.03397 [nucl-ex] .
  * Kotila (2023) J. Kotila, Private communication (2023).
  * Nomura (2022) K. Nomura, Phys. Rev. C 105, 044301 (2022), arXiv:2203.13042 [nucl-th] .
  * Barea  _et al._ (2015) J. Barea, J. Kotila, and F. Iachello, Phys. Rev. C91, 034304 (2015), arXiv:1506.08530 [nucl-th] .
  * Popara  _et al._ (2022) N. Popara, A. Ravlić, and N. Paar, Phys. Rev. C 105, 064315 (2022), arXiv:2107.08747 [nucl-th] .
  * Coraggio  _et al._ (2019) L. Coraggio, L. De Angelis, T. Fukui, A. Gargano, N. Itaco, and F. Nowacki, Phys. Rev. C 100, 014316 (2019), arXiv:1812.04292 [nucl-th] .
  * Ejiri  _et al._ (2019) H. Ejiri, J. Suhonen, and K. Zuber, Phys. Rept. 797, 1 (2019).


# Production of dileptonic bound states in the Higgs boson decay

F. A. Martynenko<EMAIL_ADDRESS>Samara University  A. P. Martynenko
<EMAIL_ADDRESS>Samara University  A. V. Eskin
<EMAIL_ADDRESS>Samara University

###### Abstract

The production of single and paired lepton bound states in the decay of the
Higgs boson has been studied. We explore different decay mechanisms that
contribute significantly to the decay width. The decay widths are calculated
taking into account relativistic corrections in the decay amplitude and in the
wave function of the bound state of leptons.

Higgs boson, positronium, dimuonium, ditauonium

##  I Introduction

After the discovery of the Higgs boson, a period of more detailed study of its
properties and decay characteristics through various channels began. Such
theoretical research now relies on existing experimental data, the collection
of which is steadily expanding [1]. Plans have begun to emerge for various
future accelerators in which there will be a significant increase in the
production of Higgs bosons, and, therefore, the possibility of more accurately
determining various parameters of the Higgs sector will open up [2, 3]. As
often happens when formulating programs for future colliders, the task is set
to search for events that do not fit into the scheme of the Standard Model.
Among the rare decays of the Higgs boson, one can distinguish a class of
processes in which bound states of fundamental particles arise, such as mesons
and baryons [4, 5, 6, 7]. The theoretical interest in these reactions is due
to the fact that, among other things, they can additionally test models for
describing bound states of quarks. The number of such works has grown
significantly over the past 10 years, and the methods for calculating the
observed characteristics have received significant development [8]. Another
group of reactions for the formation of bound states of particles consists of
those reactions in which bound states of leptons (positronium, dimuonium,
ditauonium) can be produced. In reactions of this type there are no
uncertainties caused by the nonperturbative interaction of quarks. Various
processes of production and decay of positronium already have a rich history,
which included critical periods of divergence between theory and experiment.
Other bound states of leptons have been studied much less. Thus, it can be
noted that calculations of the characteristics of dimuonium have been carried
out quite a long time ago with great accuracy, but the bound state itself has
not yet been observed [9]. Various concepts of dimuonium and ditauonium
production are discussed in many papers [10, 11, 12]. One possibility for the
production of lepton bound states is connected with the decay of the Higgs
boson. In this work, we study the processes of both single and pair production
of dimuonium and ditauonium within the framework of relativistic method of
describing the production of bound states of particles. This method was used
previously to calculate the production of charmonium and $B_{c}$ mesons [13,
14, 15]. It allows one to obtain values of the decay width taking into account
relativistic effects. It is useful to note that the observation of
orthodimuonium or orthoditauonium is connected with subsequent decay into
three photons.

##  II Production of single leptonic bound states

Let us first consider the single production of lepton bound states in the
decay of $H\to\gamma+\bar{l}l$, $H\to Z+\bar{l}l$. The Higgs boson-lepton
vertex is determined by the factor [16]:

$m(\sqrt{2}G_{F})^{1/2}=\frac{e}{\sin 2\theta_{W}}\frac{m}{M_{Z}},$ (1)

where $m$ is the lepton mass, $M_{Z}$ is the Z-boson mass, $\theta_{W}$ is the
Weinberg angle, $G_{F}$ is the Fermi coupling constant.

Figure 1: Different mechanisms of single dilepton production in the Higgs
boson decay.

The lepton-photon (lepton-Z-boson) mechanism of single production of dimuonium
and ditauonium is determined by three amplitudes presented in Fig. 1(a,b,c).
The initial expression for the amplitude of such a process can be represented
as a convolution of the amplitude of the production of a photon and a lepton-
antilepton pair with the Bethe-Salpeter wave function, which describes, within
the framework of the four-dimensional formalism, a lepton-antilepton bound
state moving with 4-momentum $P$. Transforming the amplitude during the
transition to three-dimensional formalism and using the law of transformation
of the quasipotential wave function during the transition from the moving
frame of reference to the rest frame of the bound state of leptons, we obtain
two decay amplitudes (Fig. 1(a,b)) in the form:

$\displaystyle{\cal M}_{1}(k,P)=\frac{4\pi\alpha m_{l}}{\sin
2\theta_{W}M_{Z}}\int\frac{d{\bf p}}{(2\pi)^{3}}\frac{\Psi({\bf
p})}{\frac{\varepsilon}{m}\frac{(\varepsilon+m)}{2m}[(r-p_{2})^{2}-m^{2}]}\times$
$\displaystyle
Tr\bigl{\\{}[\frac{\hat{v}-1}{2}-\hat{v}\frac{p^{2}}{2m(\varepsilon+m)}-\frac{\hat{p}}{2m}]\hat{\varepsilon}_{l\bar{l}}\frac{(\hat{v}+1)}{2\sqrt{2}}[\frac{\hat{v}+1}{2}-\hat{v}\frac{p^{2}}{2m(\varepsilon+m)}+\frac{\hat{p}}{2m}]\hat{\varepsilon}_{\gamma}(\hat{r}-\hat{p}_{2}+m)\bigr{\\}},$
(2) $\displaystyle{\cal M}_{2}(k,P)=\frac{4\pi\alpha m_{l}}{\sin
2\theta_{W}M_{Z}}\int\frac{d{\bf p}}{(2\pi)^{3}}\frac{\Psi({\bf
p})}{\frac{\varepsilon}{m}\frac{(\varepsilon+m)}{2m}[(p_{1}-r)^{2}-m^{2}]}\times$
$\displaystyle
Tr\bigl{\\{}[\frac{\hat{v}-1}{2}-\hat{v}\frac{p^{2}}{2m(\varepsilon+m)}-\frac{\hat{p}}{2m}]\hat{\varepsilon}_{l\bar{l}}\frac{(\hat{v}+1)}{2\sqrt{2}}[\frac{\hat{v}+1}{2}-\hat{v}\frac{p^{2}}{2m(\varepsilon+m)}+\frac{\hat{p}}{2m}](\hat{p}_{1}-\hat{r}+m)\hat{\varepsilon}_{\gamma}\bigr{\\}},$
(3)

where $r$ is the Higgs boson four-momentum, $\varepsilon_{\gamma}$ is the
photon polarization vector, $\varepsilon_{l\bar{l}}$ is the polarization
vector of dimuonium (ditauonium), $k$ is the photon four-momentum,
$\varepsilon(p)=\sqrt{p^{2}+m^{2}}$, $v=P/M$, $M$ is the mass of $l\bar{l}$
state. $p$ is the lepton relative four-momentum. The expressions for the
amplitudes (II), (II) are obtained as a result of the Lorentz transformations
in the same way as in our works [13, 14, 15]. They make it possible to take
into account relativistic corrections connected with the momentum of relative
motion of leptons $p$.

The simplifications that can be made in the amplitudes are primarily related
to the denominators of the lepton propagators. The magnitude of the relative
momentum is small compared to the mass of the Higgs boson $M_{H}$, so we have

$(p_{1}-r)^{2}-m^{2}\approx(r-p_{2})^{2}-m^{2}\approx\frac{1}{2}(M_{H}^{2}-M^{2}),$
(4)

where the lepton mass can be approximately replaced here by $m=\frac{M}{2}$
neglecting bound state effects.

After calculating the trace and extracting relativistic corrections in $p/m$
of second and fourth order in the numerator, we obtain the following
expression for the total decay amplitude:

$\displaystyle{\cal M}=\frac{64\pi\alpha r_{1}}{\sin
2\theta_{W}M_{Z}(r_{2}^{2}-1)}\tilde{\psi}(0)\Biggl{\\{}(\varepsilon_{\gamma}\varepsilon_{l\bar{l}})\Bigl{[}r_{2}^{2}\bigl{(}\frac{1}{2}-\frac{1}{6}\omega_{1}\bigr{)}-r_{1}+r_{1}\omega_{1}+\frac{4}{3}r_{1}\omega_{2}\Bigr{]}-$
$\displaystyle(v_{\gamma}\varepsilon)(v\varepsilon_{\gamma})\Bigl{[}1-\frac{1}{3}\omega_{1}\Bigr{]}\Biggr{\\}},$
(5)

where for convenience the particle mass ratios are introduced:
$r_{1}=\frac{m}{M}$, $r_{2}=\frac{M_{H}}{M}$.

Relativistic corrections in (II) of the second and fourth orders in $p/m$ are
determined in our approach by the parameters $\omega_{1}$, $\omega_{2}$ and
the value of the wave function at zero $\tilde{\psi}(0)$. The definition of
these parameters and their numerical calculation is given in the next section.
Note that we have retained the 4th order $\alpha$ corrections here only to
demonstrate the overall amplitude structure. In what follows, when obtaining
numerical estimates, we neglect fourth-order corrections.

Using the decay amplitude (II), we can calculate the differential decay width
according to the formula:

$d\Gamma=\frac{|{\bf P}|}{32\pi^{2}M_{H}^{2}}\overline{|M|^{2}}d\Omega.$ (6)

Three-momentum of final dilepton state is expressed in terms of masses of
higgs boson $(M_{H})$, $Z$-boson $(M_{Z})$, and dilepton state $(M)$:

$|{\bf
P}|=\begin{cases}\sqrt{\frac{M^{4}+(M_{H}^{2}-M_{Z}^{2})^{2}-2M^{2}(M_{H}^{2}+M_{Z}^{2})}{4M_{H}^{2}}},~{}~{}~{}H\rightarrow(l{\bar{l}})+Z\\\
\frac{M_{H}^{2}-M^{2}}{2M_{H}},~{}~{}~{}H\rightarrow(l{\bar{l}})+\gamma\\\
\sqrt{\frac{M_{H}^{2}}{4}-M^{2}},~{}~{}~{}H\rightarrow(l{\bar{l}})+(l{\bar{l}})\end{cases}.$
(7)

In the numerator of the amplitude (II) we define two functions $g_{1,2}$ as
follows:

$\displaystyle
N=g_{1}(\varepsilon_{l\bar{l}}\varepsilon_{\gamma})-g_{2}(v\varepsilon_{\gamma})(v_{\gamma}\varepsilon_{l\bar{l}}),$
(8) $\displaystyle
g_{1}=\Bigl{[}r_{2}^{2}\bigl{(}\frac{1}{2}-\frac{1}{6}\omega_{1}\bigr{)}-r_{1}+r_{1}\omega_{1}+\frac{4}{3}r_{1}\omega_{2}\Bigr{]},~{}~{}~{}g_{2}=\Bigl{[}1-\frac{1}{3}\omega_{1}\Bigr{]}.$

After summing over the polarizations of the photon and ortho-dileptonium we
obtain total decay width $H\to\gamma+(l\bar{l})$ in the form:

$\displaystyle\Gamma_{\gamma(l\bar{l})}=\frac{512\pi\alpha^{2}r_{1}^{2}}{r_{2}^{3}\sin^{2}2\theta_{W}M_{Z}^{2}(r_{2}^{2}-1)}\Biggl{[}\frac{1}{2}r_{2}^{2}g_{2}^{2}+3g_{1}^{2}-\frac{1}{4}g_{2}^{2}-\frac{1}{4}r_{2}^{4}g_{2}^{2}\Biggr{]}|\tilde{\psi}(0)|^{2}.$
(9)

The decay width of the Higgs boson into $Z$-boson and $(l\bar{l})$ (see Fig.
1(a,b)) is calculated similarly and has the following form:

$\displaystyle\Gamma_{Z(l\bar{l})}=\frac{32\pi\alpha^{2}\sqrt{[(r_{2}+1)^{2}-r_{3}^{2}][(r_{2}-1)^{2}-r_{3}^{2}]}}{r_{2}^{3}\sin^{4}2\theta_{W}M_{Z}^{2}(r_{2}^{2}-1)^{2}}|\tilde{\psi}(0)|^{2}\Biggl{\\{}g_{2Z}^{2}(-\frac{1}{4}+\frac{1}{16}r_{3}^{-2}+\frac{3}{8}r_{3}^{2}-\frac{1}{4}r_{3}^{4}+$
$\displaystyle\frac{1}{16}r_{3}^{6}-\frac{1}{4}r_{2}^{2}r_{3}^{-2}+\frac{1}{4}r_{2}^{2}+\frac{1}{4}r_{2}^{2}r_{3}^{2}-\frac{1}{4}r_{2}^{2}r_{3}^{4}+\frac{3}{8}r_{2}^{4}r_{3}^{-2}+\frac{1}{4}r_{2}^{4}+\frac{3}{8}r_{2}^{4}r_{3}^{2}-\frac{1}{4}r_{2}^{6}r_{3}^{-2}-\frac{1}{4}r_{2}^{6}+\frac{1}{16}r_{2}^{8}r_{3}^{-2})+$
$\displaystyle
g_{1Z}g_{2Z}(-\frac{1}{4}+\frac{1}{4}r_{3}^{-2}-\frac{1}{4}r_{3}^{2}+\frac{1}{4}r_{3}^{4}\frac{3}{4}r_{2}^{2}r_{3}^{-2}-\frac{1}{2}r_{2}^{2}-3/4r_{2}^{2}r_{3}^{2}+\frac{3}{4}r_{2}^{4}r_{3}^{-}2+\frac{3}{4}r_{2}^{4}-\frac{1}{4}r_{2}^{6}r_{3}^{-2})+$
$\displaystyle
g_{1Z}^{2}(\frac{5}{2}+\frac{1}{4}r_{3}^{-}2+\frac{1}{4}r_{3}^{2}-\frac{1}{2}r_{2}^{2}r_{3}^{-2}-\frac{1}{2}r_{2}^{2}+\frac{1}{4}r_{2}^{4}r_{3}^{-2})\Biggr{\\}},$
(10) $\displaystyle
g_{1Z}=(r_{2}^{2}-r_{3}^{2})(\frac{1}{8}+\frac{1}{4}a_{z}-\frac{1}{24}\omega_{1}-\frac{1}{12}\omega_{1}a_{z})+r_{1}(-\frac{1}{4}-\frac{1}{2}a_{z}+\frac{1}{4}\omega_{1}+\frac{1}{2}\omega_{1}a_{z}),$
$\displaystyle
g_{2Z}=(\frac{1}{4}+\frac{1}{2}a_{z}-\frac{1}{12}\omega_{1}-\frac{1}{6}a_{z}\omega_{1}),$
(11)

where we introduce another mass coefficient $r_{3}=\frac{M_{Z}}{M}$,
$a_{z}=2\sin^{2}\theta_{W}$.

To obtain total width of the decay $H\to Z+(l\bar{l})$, including the
contribution of the $ZZ$-production mechanism (see, Fig. 1(c)), it is
necessary to make the following substitution in (II):

$\displaystyle g_{1Z}\to
g_{1Z}+\frac{r_{2}^{2}}{r_{1}}g_{1ZZ},~{}~{}~{}g_{1ZZ}=-\frac{1}{4}-\frac{1}{2}a_{z}-\frac{1}{12}\omega_{1}-\frac{1}{6}\omega_{1}a_{z}.$
(12)

The expressions (8), (9)-(II) clearly show that the main parameters that
determine numerical values of decay branchings include, together with the fine
structure constant, the particle mass ratios $r_{1}$, $r_{2}$, $r_{3}$.

Figure 2: Quark loop and $W$-boson loop mechanisms of single dilepton bound
state production in the Higgs boson decay.

The studies we carried out in [13, 14] show that it is necessary to take into
account the loop ($W$-boson and quark) production mechanism $Z+(l\bar{l})$
shown in Fig. 4. Although formally an additional degree of $\alpha$ appears in
it compared to the $ZZ$-boson mechanism, the fine structure constant is not a
decisive factor in determining the order of magnitude of the contribution.
Numerous particle mass ratio parameters begin to play a more significant role
here. The loops of heavy quarks and the $W$-boson make their contribution to
the structure functions $g_{1}$ and $g_{2}$, which can be represented as
follows:

$g^{(loop)}_{1~{}Z}=\frac{\pi\alpha
r_{2}^{2}r_{3}^{2}\cos^{2}\theta_{W}}{2r_{1}}A^{Z}_{W}(r_{2}^{2}-r_{3}^{2})(1+\frac{1}{3}\omega_{1})+\sum_{Q=c,b,t}\frac{3\alpha
q_{Q}r_{2}^{2}m_{Q}^{2}}{r_{1}M^{2}}A^{Z}_{Q}(r_{2}^{2}-r_{3}^{2})(1+\frac{1}{3}\omega_{1}),$
(13) $g^{(loop)}_{2~{}Z}=\frac{\pi\alpha
r_{2}^{2}r_{3}^{2}\cos^{2}\theta_{W}}{r_{1}}A^{Z}_{W}(1+\frac{1}{3}\omega_{1})+\sum_{Q=c,b,t}\frac{6\alpha
q_{Q}r_{2}^{2}m_{Q}^{2}}{r_{1}M^{2}}A^{Z}_{Q}(1+\frac{1}{3}\omega_{1}),$ (14)

where relative coefficients with respect to the lepton-Z-boson mechanism are
indicated. The determination of the loop parameters $A^{Z}_{W}$, $A^{Z}_{Q}$
is described below as in the case of the production of $\gamma+(l\bar{l})$.
Numerical results of calculating the decay widths with the production of
$Z+(l\bar{l})$ are shown in Table 1.

The lepton-photon decay mechanism in Fig. 1 gives a decay width of fifth order
in $\alpha$. But our analysis shows that the lepton-photon mechanism in Fig. 2
has a significantly larger numerical contribution, although it has a higher
order $\alpha^{7}$. This occurs due to simultaneously changing mass
coefficients in the amplitude of such a decay. We previously studied a similar
decay mechanism in connection with the production of heavy quarkonia in [13,
14, 15].

The general structure of the tensor corresponding to the $W$-boson or quark
loops in Fig. 2 in the case of two external virtual photons is the following:

$T^{\mu\nu}_{Q,W}=A_{Q,W}(t)(g^{\mu\nu}(v_{1}v_{2})-v_{1}^{\nu}v_{2}^{\mu})+B_{Q,W}(t)[v_{1}^{2}v_{2}^{\mu}-v_{1}^{\mu}(v_{1}v_{2})][v_{2}^{2}v_{1}^{\nu}-v_{2}^{\nu}(v_{1}v_{2})],$
(15)

where $t=\frac{M_{h}^{2}}{4m_{Q}^{2}}$ or $t=\frac{M_{h}^{2}}{4m_{W}^{2}}$,
$m_{W}$ is the mass of $W$-boson, $m_{Q}$ is the mass of heavy quark in the
quark loop. The structure functions $A_{Q,W}(t)$, $B_{Q,W}(t)$ can be obtained
using an explicit expression for a loop integrals [17, 18, 19]. In the case
where one photon is real, the only contribution comes from the structure
function $A(t)$ which satisfies the dispersion relation:

$A_{W}(t)=A_{W}(0)+\frac{t}{\pi}\int_{1}^{\infty}\frac{ImA(t^{\prime})dt^{\prime}}{t^{\prime}(t^{\prime}-t+i0)}.$
(16)

Imaginary part $ImA(t)$ can be calculated using the Mandelstam-Cutkosky rule
[20]:

$\displaystyle
ImA_{W}=\frac{r_{4}^{2}}{64\pi}\frac{1}{t(4t-r_{4}^{2})^{2}}\Biggl{[}r_{4}^{2}\sqrt{t(t-1)}(r_{4}^{2}(2t+1)-4t-6)+$
(17) $\displaystyle
4t(6-12t+r_{4}^{2}(2t+3)-r_{4}^{4})arcsh(\sqrt{t-1})\Biggr{]},~{}~{}~{}r_{4}=\frac{M}{M_{W}}.$

Since the mass parameter $r_{4}$ is very small, to calculate the function
$A(t)$ itself from its imaginary part, one can use the expansion in $r_{4}$.
In leading order by $r_{4}$ we get:

$\displaystyle
A_{W}(t)=\frac{r_{4}^{2}}{16\pi^{2}}\left[2+\frac{3}{t}+\frac{3}{t^{2}}(2t-1)f^{2}(t)\right],$
(18) $\displaystyle f(t)=\begin{cases}\arcsin\sqrt{t},~{}~{}~{}t\leq 1\\\
\frac{i}{2}\left[\ln\frac{1-\sqrt{1-t^{-1}}}{1+\sqrt{1-t^{-1}}}-i\pi\right],~{}~{}~{}t>1\end{cases}.$
(19)

Since the contribution of the $W$-boson loop mechanism to the decay width
significantly exceeds the others, it is convenient to separately present the
corresponding formula for the decay width:

$\Gamma=\frac{8\pi^{3}\alpha^{4}(r_{2}^{2}-1)ctg^{2}\theta_{W}M_{Z}^{2}}{r_{2}^{3}M^{4}}|\tilde{\psi}(0)|^{2}A_{W}^{2}\bigl{[}3|g_{1W}|^{2}-\frac{1}{4}(r_{2}^{2}-1)^{2}|g_{2W}|^{2}\bigr{]},$
(20) $\displaystyle
g_{1W}=(r_{2}^{2}-1)(1+\frac{7}{3}\omega_{1}+\frac{11}{3}\omega_{2}),~{}~{}g_{2W}=2(1+\frac{7}{3}\omega_{1}+\frac{11}{3}\omega_{2}).$
(21)

We emphasize that this formula contains mass factors constructed from the
observed particle masses. Thus, the coupling effects in the lepton system are
taken into account despite their smallness.

Similarly to the $W$-loop mechanism, we can calculate the contribution of the
quark loop mechanism shown in Fig. 2. The imaginary part of the structure
function $A_{Q}(t)$ has the form:

$\displaystyle
ImA_{Q}=\frac{r_{5}^{2}}{32\pi}\frac{1}{t(4t-r_{5}^{2})^{3}}\Biggl{[}3r_{5}^{2}\sqrt{t(t-1)}(4t-r_{5}^{2})+$
(22) $\displaystyle
4t(r_{5}^{4}+2r_{5}^{2}(1-2t)+8(t-1)t)arcsh(\sqrt{t-1})\Biggr{]},~{}~{}~{}r_{5}=\frac{M}{m_{Q}}.$

After calculating the dispersion integral, the following result is obtained
for the function $A_{Q}(t)$ in leading order in $r_{5}^{2}$ in the form:

$\displaystyle
A_{Q}(t)=\frac{r_{5}^{2}}{16\pi^{2}}\left[\frac{1}{t}+\frac{(t-1)}{t^{2}}f^{2}(t)\right].$
(23)

Total decay width, taking into account quark loops, can then be obtained from
formula (20), making the following substitution:

$g_{1W,2W}\to
g_{1W,2W}\left(1+\sum_{Q=c,b,t}\frac{24q_{Q}^{2}m_{Q}^{2}}{\cos^{2}{\theta_{W}}M_{Z}^{2}}\frac{A_{Q}}{A_{W}}\right).$
(24)

Numerical values for the width of single dilepton state production are
discussed below along with paired production processes.

##  III Pair production of leptonic bound states

In addition to the processes of single production of bound states of leptons,
we also consider processes of pair production of orto-positronium
($e^{+}e^{-}$), orto-dimuonium ($\mu^{+}\mu^{-}$), orto-ditauonium
($\tau^{+}\tau^{-}$) in the Higgs boson decay. There are five different
production mechanisms presented in Figs. 3, 4, 5, 6, 7, which we consider
below.

###  III.1 Lepton - photon production mechanism A

The first lepton-photon production mechanism can be separated into two stages.
On the first stage two lepton - antilepton pair are produced. The first pair
is produced in direct Higgs boson decay. After that lepton or antilepton can
emit a photon or Z-boson, that produces another lepton - antilepton pair. On
the second stage free leptons are combined in bound states.

The width of the decay process $H\to(l\bar{l})_{S=1}+(l\bar{l})_{S=1}$ is
determined by an expression (6). Four-momenta $p_{1}$, $p_{2}$ of one lepton-
antilepton pair and $q_{1}$, $q_{2}$ of second pair can be expressed in terms
of total $P$, $Q$ and relative $p$, $q$ 4-momenta in the form:

$p_{1,2}=\frac{1}{2}P\pm p,~{}~{}~{}q_{1,2}=\frac{1}{2}Q\pm
q,~{}~{}~{}(pP)=(qQ)=0.$ (25)

Total amplitude of pair production ${\cal M}(P,Q)$ now contains two
convolutions of the Higgs boson decay amplitude with the wave functions of two
lepton bound states:

$\displaystyle{\cal
M}(P,Q)^{(A)}=4\pi\alpha\Gamma_{Hl{\bar{l}}}\int\frac{d{\bf
p}}{(2\pi)^{3}}\int\frac{d{\bf
q}}{(2\pi)^{3}}Tr\Bigl{\\{}\Psi(p,P)\frac{-{\hat{r}}+{\hat{p}}_{1}+m}{(-r+p_{1})^{2}-m^{2}}\gamma^{\mu}\Psi(q,Q)\gamma^{\nu}+$
$\displaystyle\Psi(p,P)\gamma^{\mu}\frac{{\hat{r}}-{\hat{q}}_{1}+m}{(r-q_{1})^{2}-m^{2}}\Psi(q,Q)\gamma^{\nu}\Bigr{\\}}D^{\mu\nu}_{\gamma}(k_{1}),$
(26)

where $\Gamma_{Hl{\bar{l}}}=m(\sqrt{2}G_{F})^{1/2}$ is a vertex of the Higgs-
lepton-antilepton interaction, $D^{\mu\nu}_{\gamma}(k_{1})$ is a photon
propagator, $k_{1}=p_{2}+q_{2}$.

Figure 3: Lepton - photon mechanism A of pair dilepton production in Higgs
boson decay.

Taking into account the law of transformation of the wave function of a bound
state upon transition to a rest frame, we obtain wave functions in (III.1) in
the form [13, 14]:

$\displaystyle\Psi(p,P)=\frac{\psi_{0}({\bf
p})}{\frac{\varepsilon(p)}{m}\frac{\varepsilon(p)+m}{2m}}\Bigl{[}\frac{{\hat{v}}_{1}-1}{2}-{\hat{v}}_{1}\frac{p^{2}}{2m(\varepsilon(p)+m)}-\frac{{\hat{p}}}{2m}\Bigr{]}\times$
$\displaystyle\frac{1}{2\sqrt{2}}{\hat{\varepsilon}}_{1}(P,S_{z})(1+{\hat{v}}_{1})\Bigl{[}\frac{{\hat{v}}_{1}+1}{2}-{\hat{v}}_{1}\frac{p^{2}}{2m(\varepsilon(p)+m)}+\frac{{\hat{p}}}{2m}\Bigr{]},$
(27) $\displaystyle\Psi(q,Q)=\frac{\psi_{0}({\bf
q})}{\frac{\varepsilon(q)}{m}\frac{\varepsilon(q)+m}{2m}}\Bigl{[}\frac{{\hat{v}}_{2}-1}{2}-{\hat{v}}_{2}\frac{q^{2}}{2m(\varepsilon(q)+m)}+\frac{{\hat{q}}}{2m}\Bigr{]}\times$
$\displaystyle\frac{1}{2\sqrt{2}}{\hat{\varepsilon}}_{2}(Q,S_{z})(1+{\hat{v}}_{2})\Bigl{[}\frac{{\hat{v}}_{2}+1}{2}-{\hat{v}}_{2}\frac{q^{2}}{2m(\varepsilon(q)+m)}-\frac{{\hat{q}}}{2m}\Bigr{]},$
(28)

where we introduce two auxiliary 4-vectors $v_{1}=\frac{P}{M}$,
$v_{2}=\frac{Q}{M}$, $\varepsilon(p)=\sqrt{{\bf p}^{2}+m^{2}}$, $m$ is the
lepton mass. In equations (III.1), (III.1) the projection operators on
$(^{3}S_{1})$ dilepton states are introduced. Projection operators are
constructed in terms of the Dirac bispinors in the rest frame as follows:

$\Pi_{S=1}=\frac{(1+\hat{v}_{1})}{2\sqrt{2}}{\hat{\varepsilon}}_{1}(P,S_{z}).$
(29)

The denominators of propagators in amplitude (III.1) can be simplified due to
the condition $M_{h}\gg m$ as follows:

$\displaystyle\frac{1}{k_{1}^{2}}=\frac{1}{(p_{2}+q_{2})^{2}}\approx\frac{4}{M_{h}^{2}},~{}~{}~{}\frac{1}{(-r+p_{1})^{2}-m^{2}}\approx\frac{1}{(r-q_{1})^{2}-m^{2}}\approx\frac{2}{M_{h}^{2}}.$
(30)

The calculation of trace in (III.1) of the product of $\gamma$-factors and
various simplifications are performed in the Form package [21]. Considering
small factors $\frac{{\bf p}}{m}$, $\frac{{\bf q}}{m}$, we perform an
expansion on them in the numerator of the amplitude, while preserving the
second-order terms. Then the numerator of this amplitude takes the form:

$\displaystyle
N^{(A)}=\left(g_{1}^{(A)}g^{\alpha\beta}-g_{2}^{(A)}v_{1}^{\beta}v_{2}^{\alpha}\right)\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}),$
(31)

where the introduced functions are equal:

$\displaystyle
g_{2}^{(A)}=1-\frac{1}{9}\omega_{1}^{2},~{}~{}g_{1}^{(A)}=-\frac{1}{2}+\frac{1}{2}r_{2}^{2}-r_{1}+\frac{2}{3}\omega_{1}r_{1}+\frac{1}{18}\omega_{1}^{2}-\frac{1}{18}\omega_{1}^{2}r_{2}^{2}+\frac{1}{3}\omega_{1}^{2}r_{1}.$
(32)

Special relativistic parameters $\omega_{1}$, $\tilde{\psi}(0)$ take into
account relativistic corrections connected with the relative motion of
leptons. The calculation of such parameters is discussed in details below in
section IV.1. After all transformations, the production amplitude (III.1)
takes the form:

$\displaystyle{\cal M}(P,Q)^{(A)}=\frac{64\pi\alpha
mM(\sqrt{2}G_{F})^{1/2}{\tilde{\psi}}^{2}(0)}{M_{H}^{4}}(g_{1}^{(A)}g^{\alpha\beta}-g_{2}^{(A)}v_{1}^{\beta}v_{2}^{\alpha}\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}).~{}~{}~{}$
(33)

###  III.2 Lepton - Z-boson production mechanism B

Figure 4: Lepton - $Z$-boson mechanism B of pair dilepton production in the
Higgs boson decay.

The second mechanism of pair production differs from the previous one by
replacing the photon with a $Z$-boson. The decay amplitudes of the Higgs boson
are shown in Fig. 4. The total amplitude of the process is determined in this
case by the following expression:

$\displaystyle{\cal M}(P,Q)^{(B)}=-i\Gamma_{Hl{\bar{l}}}\int\frac{d{\bf
p}}{(2\pi)^{3}}\int\frac{d{\bf
q}}{(2\pi)^{3}}Tr\Bigl{\\{}\Psi(p,P)\frac{-{\hat{r}}+{\hat{p}}_{1}+m}{(-r+p_{1})^{2}-m^{2}}\Gamma^{\mu}_{Zl{\bar{l}}}\Psi(q,Q)\Gamma^{\nu}_{\gamma
l{\bar{l}}}+$ $\displaystyle\Psi(p,P)\Gamma^{\mu}_{\gamma
l{\bar{l}}}\frac{{\hat{r}}-{\hat{q}}_{1}+m}{(r-q_{1})^{2}-m^{2}}\Psi(q,Q)\Gamma^{\nu}_{{Zl{\bar{l}}}}\Bigr{\\}}D^{\mu\nu}_{Z}(k_{2}),$
(34)

where $D^{\mu\nu}_{Z}(k_{2})$ is the $Z$-boson propagator,
$k_{2}=p_{2}+q_{2}$. The $Z$-boson - lepton-antilepton interaction vertex is
determined by:

$\displaystyle\Gamma^{\mu}_{Zl{\bar{l}}}=\frac{e}{\sin{2\theta_{W}}}\gamma_{\mu}\left[\frac{1}{2}(1-\gamma^{5})+2\sin^{2}{\theta_{W}}\right].$
(35)

When transforming the denominators of the $Z$-boson and lepton propagators, we
also neglect the relative momenta $p$ and $q$ and write them in the form:

$\displaystyle\frac{1}{k_{2}^{2}-M_{Z}^{2}}=\frac{1}{(p_{2}+q_{2})^{2}}\approx\frac{4}{M_{H}^{2}-4M_{Z}^{2}},~{}\frac{1}{(-r+p_{1})^{2}-m^{2}}\approx\frac{1}{(r-q_{1})^{2}-m^{2}}\approx\frac{2}{M_{H}^{2}}.$
(36)

After the trace calculation the numerator of total amplitude has the same
general structure as (31):

$\displaystyle
N^{(B)}=\left(g_{1}^{(B)}g^{\alpha\beta}-g_{2}^{(B)}v_{1}^{\beta}v_{2}^{\alpha}\right)\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}),$
(37)

and the functions $g_{1}^{(B)}$, $g_{2}^{(B)}$ in the considered approximation
are equal:

$\displaystyle g_{2}^{(B)}=2(1-2a_{z}+2a_{z}^{2})(9-\omega_{1}^{2}),$ (38)
$\displaystyle
g_{1}^{(B)}=(1-2a_{z}+2a_{z}^{2})(3+\omega_{1})(-3+\omega_{1}-6r_{1}+6\omega_{1}r_{1}+3r_{2}^{2}-\omega_{1}r_{2}^{2}),~{}~{}~{}~{}~{}~{}~{}$
(39)

where $a_{z}=2\sin^{2}\theta_{W}$. The resulting expression for total
amplitude of this decay mechanism has the form:

$\displaystyle{\cal M}(P,Q)^{(B)}=\frac{64\pi\alpha
mM(\sqrt{2}G_{F})^{1/2}{\tilde{\psi}}^{2}(0)}{M_{H}^{2}(M_{H}^{2}-4M_{Z}^{2})\sin^{2}2\theta_{W}}(g_{1}^{(B)}g^{\alpha\beta}-g_{2}^{(B)}v_{1}^{\beta}v_{2}^{\alpha})\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}).~{}~{}~{}$
(40)

###  III.3 Lepton - photon production mechanism C

There are two more amplitudes in the lepton-photon decay mechanism, which are
shown in Fig. 5. Their structure is different from the previous ones in
subsection III.1, so we consider them separately.

Figure 5: Lepton - photon mechanism C of pair dilepton production in the Higgs
boson decay.

The first difference is that here there is a product of two different traces.
The second difference is that the denominator of the photon propagator does
not contain the large mass of the Higgs boson or $Z$-boson. The denominator of
the photon propagator $k_{3}^{2}$ ($k_{3}=q_{1}+q_{2}$), which is equal to the
mass of the bound state of leptons, can be considered as a factor that
increases the contribution of such a mechanism to the overall decay amplitude.
Let us present here total amplitude of the process, which is obtained
according to the same rules as the previous one in Fig. 3:

$\displaystyle{\cal M}(P,Q)^{(C)}=-i\Gamma_{Hl{\bar{l}}}\int\frac{d{\bf
p}}{(2\pi)^{3}}\int\frac{d{\bf
q}}{(2\pi)^{3}}Tr\Bigl{\\{}\Psi(p,P)\frac{(-{\hat{r}}+{\hat{p}}_{1}+m)}{(-r+p_{1})^{2}-m^{2}}\Gamma^{\mu}_{\gamma
l{\bar{l}}}+$ (41)
$\displaystyle\Psi(p,P)\gamma^{\mu}\frac{({\hat{r}}-{\hat{q}}_{1}+m)}{(r-p_{2})^{2}-m^{2}}\Bigr{\\}}Tr\Bigl{\\{}\Psi(q,Q)\gamma^{\nu}\Bigr{\\}}D^{\mu\nu}_{\gamma}(k_{3}).$

Now the mass of the Higgs boson is contained in only one of the denominators
in the amplitude (41) in the form:

$\displaystyle\frac{1}{k_{3}^{2}}=\frac{1}{q_{1}+q_{2}}\approx\frac{1}{M^{2}},~{}~{}~{}\frac{1}{(-r+p_{1})^{2}-m^{2}}\approx\frac{1}{(r-p_{2})^{2}-m^{2}}\approx\frac{2}{M_{H}^{2}}.$
(42)

The numerator of this amplitude is also a convolution of the ortho-dileptonic
polarization vectors with a tensor, which is expressed through two new
functions in the form:

$\displaystyle
N^{(C)}=\left(g_{1}^{(C)}g^{\alpha\beta}-g_{2}^{(C)}v_{1}^{\beta}v_{2}^{\alpha}\right)\varepsilon^{\alpha}(P,S_{z})\varepsilon^{\beta}(Q,S_{z}),$
(43) $\displaystyle
g_{1}^{(C)}=r_{2}^{2}-1-2r_{1}+\frac{4}{3}\omega_{1}r_{1}+\frac{1}{9}\omega_{1}^{2}-\frac{1}{9}\omega_{1}^{2}r_{2}^{2}+\frac{2}{3}\omega_{1}^{2}r_{1},~{}~{}~{}g_{2}^{(C)}=2-\frac{2}{9}\omega_{1}^{2}.$
(44)

Total amplitude of the production of a pair of dileptonia in this mechanism is
represented as follows:

$\displaystyle{\cal
M}(P,Q)^{(C)}=-\frac{8(4\pi\alpha)^{3/2}m{\tilde{\psi}}^{2}(0)}{M_{H}^{2}M_{Z}\sin
2\theta_{W}}(g_{1}^{(C)}g^{\alpha\beta}-g_{2}^{(C)}v_{1}^{\beta}v_{2}^{\alpha}\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}).$
(45)

###  III.4 Lepton - Z-boson production mechanism D

Replacing the photon with $Z$-boson, we obtain two more pair production
amplitudes shown in Fig. 6 But unlike the previous photon production of a
bound state of leptons, now the denominator of the $Z$-boson propagator still
contains a large mass of the $Z$-boson, which may mean the suppression of this
mechanism compared to the lepton-photon one. However, we have also included
this amplitude for completeness.

The amplitude for the production of a leptonium pair takes the form:

$\displaystyle{\cal M}(P,Q)^{(D)}=-i\Gamma_{Hl{\bar{l}}}\int\frac{d{\bf
p}}{(2\pi)^{3}}\int\frac{d{\bf
q}}{(2\pi)^{3}}Tr\Bigl{\\{}\Psi(p,P)\frac{-{\hat{r}}+{\hat{p}}_{1}+m}{(-r+p_{1})^{2}-m^{2}}\Gamma^{\mu}_{Zl{\bar{l}}}+$
$\displaystyle\Psi(p,P)\Gamma^{\mu}_{Zl{\bar{l}}}\frac{{\hat{r}}-{\hat{q}}_{1}+m}{(r-p_{2})^{2}-m^{2}}\Bigr{\\}}Tr\Bigl{\\{}\Psi(q,Q)\Gamma^{\nu}_{{Zl{\bar{l}}}}\Bigr{\\}}D^{\mu\nu}_{Z}(k_{4}),$
(46)

where $k_{4}=q_{1}+q_{2}$.

Figure 6: Lepton - $Z$-boson mechanism D of pair dilepton production in the
Higgs boson decay.

The propagator denominators contain two large masses of the $Z$-boson and the
Higgs boson:

$\displaystyle\frac{1}{k_{4}^{2}-M_{Z}^{2}}=\frac{1}{(q_{1}+q_{2})^{2}-M_{Z}^{2}}\approx\frac{1}{M^{2}-M_{Z}^{2}},~{}~{}\frac{1}{(-r+p_{1})^{2}-m^{2}}\approx\frac{1}{(r-p_{2})^{2}-m^{2}}\approx\frac{2}{M_{H}^{2}}.$
(47)

After calculating two traces in the numerator (III.4) we obtain:

$\displaystyle
N^{(D)}=\left(g_{1}^{(D)}g^{\alpha\beta}-g_{2}^{(D)}v_{1}^{\beta}v_{2}^{\alpha}\right)\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}),$
(48) $\displaystyle g_{2}^{(D)}=2(-1+2a_{z})(9-\omega_{1}^{2}),$ (49)
$\displaystyle
g_{1}^{(D)}=(-1+2a_{z})(3+\omega_{1})(-3+\omega_{1}-6r_{1}+6\omega_{1}r_{1}+3r_{2}^{2}-\omega_{1}r_{2}^{2}).$

Then the production amplitude (III.4) has the same general structure as
previous amplitudes:

$\displaystyle{\cal M}(P,Q)^{(D)}=\frac{16\pi\alpha
mM(\sqrt{2}G_{F})^{1/2}{\tilde{\psi}}^{2}(0)}{M_{H}^{2}(M^{2}-M_{Z}^{2})\sin^{2}(2\theta_{W})}(g_{1}^{(D)}g^{\alpha\beta}-g_{2}^{(D)}v_{1}^{\beta}v_{2}^{\alpha})\varepsilon_{1}^{\alpha}(P,S_{z})\varepsilon_{2}^{\beta}(Q,S_{z}).$
(50)

###  III.5 ZZ-boson production mechanism E

The production of a pair of bound states of particles (quarks, leptons) can
occur due to $ZZ$-mechanism shown in Fig. 7. A pair of $Z$-bosons arises
directly from the decay of the Higgs boson. The required vertex functions have
the form:

$\Gamma^{\alpha\beta}=\frac{2e}{\sin
2\theta_{W}}M_{Z}g^{\alpha\beta},~{}~{}~{}\Gamma^{\mu}_{Zl\bar{l}}=\frac{e}{\sin
2\theta_{W}}\gamma^{\mu}\left[\frac{1}{2}(1-\gamma_{5})+2\sin^{2}\theta_{W}\right].$
(51)

In this case, we have direct (a) and crossed (b) amplitudes (see Fig.7):

$\displaystyle{\cal
M}(P,Q)^{(E1)}=\Gamma^{\mu^{\prime}\nu^{\prime}}\int\frac{d{\bf
p}}{(2\pi)^{3}}\int\frac{d{\bf
q}}{(2\pi)^{3}}Tr\Bigl{\\{}\Psi(p,P)\Gamma^{\mu}_{Zl{\bar{l}}}\Bigr{\\}}Tr\Bigl{\\{}\Psi(q,Q)\Gamma^{\nu}_{Zl{\bar{l}}}\Bigr{\\}}D^{\mu\mu^{\prime}}_{Z}(k_{5})D^{\nu\nu^{\prime}}_{Z}(k_{6}),$
(52)

where $k_{5}=p_{1}+p_{2}$, $k_{6}=q_{1}+q_{2}$.

$\displaystyle{\cal
M}(P,Q)^{(E2)}=\Gamma^{\mu^{\prime}\nu^{\prime}}\int\frac{d{\bf
p}}{(2\pi)^{3}}\int\frac{d{\bf
q}}{(2\pi)^{3}}Tr\Bigl{\\{}\Psi(p,P)\Gamma^{\mu}_{Zl{\bar{l}}}\Psi(q,Q)\Gamma^{\nu}_{Zl{\bar{l}}}\Bigr{\\}}D^{\mu\mu^{\prime}}_{Z}(k_{7})D^{\nu\nu^{\prime}}_{Z}(k_{8}),$
(53)

where $k_{7}=p_{1}+q_{1}$, $k_{8}=p_{2}+q_{2}$.

Figure 7: Pair dilepton ZZ-boson production mechanism in the Higgs boson
decay.

By choosing the $Z$-boson propagator in the unitary gauge, we can represent
the total amplitude as follows:

$\displaystyle{\cal
M}(P,Q)^{(E)}=\frac{8(4\pi\alpha)^{3/2}M}{M_{z}^{3}\sin^{3}\theta_{W}}\tilde{\psi}^{2}(0)\Biggl{\\{}(\varepsilon_{1}\varepsilon_{2})\left[g_{1}^{E,dir}+\frac{r_{3}^{4}}{(r_{3}^{2}-\frac{r^{2}_{2}}{4})^{2}}g_{1}^{E,cr}\right]-$
(54)
$\displaystyle(\varepsilon_{1}v_{2})(\varepsilon_{2}v_{1})\left[g_{2}^{E,dir}+\frac{r_{3}^{4}}{(r_{3}^{2}-\frac{r^{2}_{2}}{4})^{2}}g_{2}^{E,cr}\right]\Biggr{\\}},$
$\displaystyle
g_{1}^{E,dir}=\frac{1}{4}+\frac{1}{6}+a_{z}+\frac{2}{3}a_{z}\omega_{1}+a_{z}^{2}+\frac{2}{3}a_{z}^{2}\omega_{1},~{}~{}~{}g_{2}^{E,dir}=0,$
(55) $\displaystyle
g_{1}^{E,cr}=-\frac{1}{4}-\frac{r_{2}^{4}}{128r_{3}^{4}}-\frac{r_{2}^{4}\omega_{1}}{192r_{3}^{4}}-\frac{r_{2}^{4}a_{z}\omega_{1}}{48r_{3}^{4}}-\frac{r_{2}^{4}a^{2}_{z}\omega_{1}}{48r_{3}^{4}}+\frac{r_{2}^{2}}{16r_{3}^{2}}+\frac{r_{2}^{2}\omega_{1}}{24r_{3}^{2}}+\frac{r_{2}^{2}a_{z}\omega_{1}}{6r_{3}^{2}}+\frac{r_{2}^{2}a^{2}_{z}\omega_{1}}{6r_{3}^{2}}-$
$\displaystyle\frac{1}{6}\omega_{1}-\frac{1}{2}a_{z}-\frac{1}{3}\omega_{1}a_{z}-\frac{1}{2}a_{z}^{2}-\frac{1}{3}\omega_{1}a_{z}^{2},$
(56) $\displaystyle
g_{2}^{E,cr}=-\frac{r_{2}^{2}}{64r_{3}^{4}}-\frac{r_{2}^{2}\omega_{1}}{96r_{3}^{4}}-\frac{r_{2}^{2}\omega_{1}a_{z}}{24r_{3}^{4}}-\frac{r_{2}^{2}\omega_{1}a^{2}_{z}}{24r_{3}^{4}}+\frac{1}{8r_{3}^{2}}+\frac{\omega_{1}}{12r_{3}^{2}}+\frac{\omega_{1}a_{z}}{3r_{3}^{2}}+\frac{\omega_{1}a_{z}^{2}}{3r_{3}^{2}}.$
(57)

###  III.6 Quark loop and W-boson loop production mechanism

From studies of paired meson production in the decay of the Higgs boson, it is
known that the $W$-boson and quark loop mechanisms can make a significant
contribution to the overall decay width. The corresponding amplitudes are
presented in Fig. 8. In each column here we have direct and crossed
amplitudes. In this case, the contribution of the crossed amplitude is
suppressed by the degree of small mass parameter, which arises from the
denominators of virtual photons.

Figure 8: Pair dilepton production mechanism in the Higgs boson decay with
quark and $W$-boson loop.

Since in the case of boson and quark loops the amplitude includes tensor (15),
we can immediately obtain total contribution of all loops in the decay
amplitude:

$\displaystyle{\cal
M}(P,Q)^{(loop)}=\frac{64\sqrt{4\pi\alpha}\pi^{2}\alpha^{2}M_{Z}ctg\theta_{W}}{M^{3}r_{2}^{4}}\tilde{\psi}^{2}(0)\Biggl{\\{}(\varepsilon_{1}\varepsilon_{2})\Bigl{(}g_{1}^{W}+\sum_{Q=c,b,t}\frac{24m_{Q}^{2}}{M_{Z}^{2}\cos^{2}\theta_{W}}g_{1}^{Q}\Bigr{)}-$
$\displaystyle(\varepsilon_{1}v_{2})(\varepsilon_{2}v_{1})\Bigl{(}g_{2}^{W}+\sum_{Q=c,b,t}\frac{24m_{Q}^{2}}{M_{Z}^{2}\cos^{2}\theta_{W}}g_{2}^{Q}\Bigr{)}\Biggr{\\}},$
(58) $\displaystyle
g_{1}^{W}=A^{W}(\frac{1}{3}\omega_{1}-\frac{1}{8}r_{2}^{2}-\frac{1}{12}\omega_{1}r_{2}^{2}-r_{2}^{4}-\frac{2}{3}\omega_{1}r_{2}^{4}+\frac{1}{2}r_{2}^{6}+\frac{1}{3}\omega_{1}r_{2}^{6})+$
(59) $\displaystyle
B^{W}(-\frac{1}{3}\omega_{1}r_{2}^{2}+\frac{1}{8}r_{2}^{4}+\frac{1}{6}\omega_{1}r_{2}^{4}-\frac{1}{32}r_{2}^{6}-\frac{1}{48}\omega_{1}r_{2}^{6}),$
(60) $\displaystyle
g_{2}^{W}=A^{W}(\frac{1}{4}+\frac{1}{6}\omega_{1}+r_{2}^{4}+\frac{2}{3}\omega_{1}r_{2}^{4})+B^{W}(-\frac{2}{3}\omega_{1}+\frac{1}{3}\omega_{1}r_{2}^{2}-\frac{17}{16}r_{2}^{4}-\frac{17}{24}\omega_{1}r_{2}^{4}).$
(61)

The coefficients $g_{1}^{Q}$, $g_{2}^{Q}$ are obtained from $g_{1}^{W}$,
$g_{2}^{W}$ by replacing $A^{W}\to A^{Q}$, $B^{W}\to B^{Q}$.

All decay amplitudes presented above, obtained within the framework of various
mechanisms, have the same structure and are expressed in terms of the
functions $g_{1}$ and $g_{2}$. Accordingly, total decay amplitude has the same
structure and contains total functions $g_{1}^{tot}$ and $g_{2}^{tot}$ with
corresponding coefficients.

Total width of the decay of the Higgs boson into a pair of bound states of
leptons can be obtained from (22) in terms of the functions $g_{1}^{tot}$,
$g_{2}^{tot}$ in the form:

$\displaystyle\Gamma=\frac{2048\pi^{4}\alpha^{5}r_{3}^{2}\sqrt{(\frac{r_{2}^{2}}{4}-1)}ctg^{2}\theta_{W}}{M^{5}r_{2}^{10}}|\tilde{\psi}(0)|^{4}\Bigl{[}|g_{1}^{tot}|^{2}(3-r_{2}^{2}+\frac{1}{4}r_{2}^{4})+$
$\displaystyle|g_{2}^{tot}|^{2}(r_{2}^{4}-\frac{1}{2}r_{2}^{6}+\frac{1}{16}r_{2}^{8})+(g_{2}^{tot}g_{1}^{\ast~{}tot}+g_{1}^{tot}g_{2}^{\ast~{}tot})(-r_{2}^{2}+\frac{3}{4}r_{2}^{4}-\frac{1}{8}r_{2}^{6})\Bigr{]},$
(62)

where the general factor corresponds to amplitude (III.6) so that the
functions $g^{W}_{1}$ and $g^{W}_{2}$ are included in the complete functions
$g_{1}^{tot}$ and $g_{2}^{tot}$ with unit coefficients.

##  IV Relativistic corrections

Relativistic corrections represent an important class of corrections for the
production of bound quark states. In the case of the production of bound
states of leptons, such corrections are small, since they are determined by
the value of the fine structure constant. Taking them into account is
connected with increasing the accuracy of calculating the decay width.

When calculating decay widths, we take into account relativistic corrections
of several types. The first group of relativistic corrections is connected
with the law of transformation of the wave function of a bound state into a
rest frame and the momenta of the relative motion of leptons in the decay
amplitude. The second group is determined by the Breit potential.

###  IV.1 Relativistic parameters in production amplitude

To calculate relativistic corrections to the $H$-boson decay amplitude, we use
the same approach as in our previous works [13, 14]. After calculating the
trace in the numerator of the amplitudes, corrections of a certain order of
smallness in the relative momentum of leptons in the bound state are
identified. We take into account second-order corrections in ${\bf
p}^{2}/m^{2}$, ${\bf q}^{2}/m^{2}$ and introduce an expansion of the form:

$\displaystyle\frac{|{\bf
p}|}{2m}=\sum_{n=1}^{\infty}\left(\frac{\varepsilon(p)-m}{\varepsilon(p)+m}\right)^{n+\frac{1}{2}},~{}~{}~{}\frac{|{\bf
q}|}{2m}=\sum_{n=1}^{\infty}\left(\frac{\varepsilon(q)-m}{\varepsilon(q)+m}\right)^{n+\frac{1}{2}}.$
(63)

This expansion emphasizes that the relative momenta of leptons are small
compared to their mass. In quantum electrodynamics, the form of the wave
function of a bound state in the initial approximation is known analytically.
It follows from this that $|{\bf p}|\sim W=\mu\alpha$.

Using expansion (63), relativistic corrections in $p$ can be represented as
the following characteristic integrals in momentum space:

$\displaystyle I^{(i)}=\int\frac{d{\bf
p}}{(2\pi)^{3}}\frac{\varepsilon(p)+m}{2\varepsilon(p)}\psi^{C}_{1S}({\bf
p})\left(\frac{\varepsilon(p)-m}{\varepsilon(p)+m}\right)^{i},$ (64)
$\psi^{C}_{1S}(p)=\frac{8\sqrt{\pi}W^{5/2}}{(p^{2}+W^{2})^{2}},~{}~{}~{}W=\frac{m}{2}\alpha.$
(65)

In the amplitudes presented earlier, relativistic effects are expressed in
terms of special parameters, which are determined by the integral ratios (64):

$\omega_{1}=\frac{I^{(1)}}{I^{(0)}},~{}~{}~{}\omega_{2}=\frac{I^{(2)}}{I^{(0)}},~{}~{}~{}I^{(0)}=\tilde{\psi}(0)=\int\frac{d{\bf
p}}{(2\pi)^{3}}\frac{\varepsilon(p)+m}{2\varepsilon(p)}\psi^{C}_{1S}({\bf
p}).$ (66)

After this, when calculating values (66), we perform an expansion in $p/m$
taking into account terms of the second degree. All nonrelativistic expansions
in $|{\bf p}|/m$ that we perform assume that $|{\bf p}|$ is small compared to
$m$. Therefore, although integrals (66) are formally convergent, when
calculating them it is necessary to take the value of the momentum at the pole
of the wave function. Proceeding in this way, we obtain the following results
when calculating relativistic corrections:

$\displaystyle\tilde{\psi}(0)=\int\frac{d{\bf
p}}{(2\pi)^{3}}\psi^{C}_{1S}({\bf p})\left[1-\frac{{\bf
p}^{2}}{4m^{2}}\right]=\psi^{C}_{1S}(0)\left[1+\frac{3}{4}\alpha^{2}\right],$
(67) $\displaystyle I^{(1)}=\int\frac{d{\bf p}}{(2\pi)^{3}}\psi^{C}_{1S}({\bf
p})\frac{{\bf
p}^{2}}{4m^{2}}=-\frac{3}{4}\alpha^{2}\psi^{C}_{1S}(0),~{}~{}~{}\omega_{1}=-\frac{3}{4}\alpha^{2}.$
(68)

The parameters $\omega_{i}$ for $i>1$ are proportional to higher powers of
$\alpha$ than (67) and are omitted from now on. The results (67) are obtained
using the following expression [22, 23]:

$\displaystyle\int\frac{d{\bf p}}{(2\pi)^{3}}\frac{{\bf
p}^{2}}{4m^{2}}\psi^{C}_{n}({\bf
p})=-\frac{[3+2(n-1)(n+1)]}{4n^{2}}\alpha^{2}\psi^{C}_{n}(0).$ (69)

###  IV.2 Relativistic corrections to the wave function

Wave function (65) was obtained from Schrödinger equation solution with
Coulomb potential. In our previous works we had estimated relativistic
correction to wave function by accounting corresponding terms in quark-
antiquark potential [13, 14, 15]. In present research we also can estimate
leading order relativistic correction to wave function as follows [9, 24, 25,
26, 27]:

$\psi^{(1)}_{1S}(0)=\int{\tilde{G}}_{1S}(0,{\bf r})\Delta V({\bf
r})\psi^{(0)}_{1S}({\bf r})d{\bf r},$ (70)

where ${\tilde{G}}_{1S}$ is reduced Coulomb Green function, $\psi^{(0)}_{1S}$
is defined by (65), $\Delta V$ is perturbative potential.

When accurately calculating the corresponding corrections to the wave
function, it is necessary to use the Green’s function with two non-zero
arguments:

$\displaystyle{\tilde{G}}_{1S}({\bf r}_{1},{\bf r}_{2})=-\frac{\mu
W}{\pi}e^{-(x_{>}+x_{<})}\Bigl{(}\frac{1}{2x_{>}}-\ln(2x_{>})-\ln(2x_{<})+Ei(2x_{<})+$
(71)
$\displaystyle\frac{7}{2}-2\gamma_{E}-(x_{>}+x_{<})+\frac{1-e^{2x_{<}}}{2x_{<}}\Bigr{)},$

where $x_{>}=max(x_{1},x_{2})$, $x_{<}=min(x_{1},x_{2})$, $x_{1.2}=Wr_{1,2}$,
$\gamma_{E}=0.577216\dots$ is Euler constant. When calculating corrections to
the wave function at zero, they are determined by the Green’s function with
one zero argument:

$\displaystyle{\tilde{G}}_{1S}({\bf r},0)=\frac{\mu
W}{4\pi}\frac{e^{-x}}{x}(4x(\ln(2x)+\gamma_{E})+4x^{2}-10x-2),$ (72)

where $x=Wr$.

The Breit potential contains several terms, the contribution of which must be
calculated:

$\displaystyle\Delta V_{1}=\frac{\pi\alpha}{m^{2}}\delta({\bf r})-\frac{{\bf
p}^{4}}{4m^{3}},$ (73) $\displaystyle\Delta
V_{2}=-\frac{\alpha}{2m^{2}}\frac{1}{r}\left({\bf p}^{2}+\frac{{\bf r}({\bf
r}{\bf p}){\bf p}}{r^{2}}\right),$ (74) $\displaystyle\Delta
V_{3}=\frac{\pi\alpha}{m^{2}}\left[\frac{7}{3}{\bf S}^{2}-2\right]\delta({\bf
r}).$ (75)

It is convenient to start calculating corrections to the wave function at zero
with the $\delta$ term of the potential:

$\displaystyle\psi^{(1)}_{1S}(0)={\tilde{G}}_{1S}(0,0)\frac{\pi\alpha}{m^{2}}\psi^{C}_{1S}(0),$
(76)

where $\tilde{G}_{1S}(0,0)$ represents a divergent expression. Along with the
coordinate representation, it is useful to use the momentum representation
when calculating corrections. In momentum representation this correction is
determined by the divergent integral in the form:

$\displaystyle\psi^{(1)}_{1}(0)=\frac{\pi\alpha}{m^{2}}\psi^{C}_{1S}(0)\int\tilde{G}_{1S}({\bf
q},{\bf p})\frac{d{\bf q}}{(2\pi)^{3}}\frac{d{\bf p}}{(2\pi)^{3}},$ (77)
$\displaystyle\tilde{G}_{1S}({\bf p},{\bf q})=-\frac{64\pi}{\alpha
W^{4}}\left[\frac{\pi^{2}W^{5}\delta({\bf p}-{\bf q})}{4({\bf
p}^{2}+W^{2})}+\frac{W^{6}}{4({\bf p}^{2}+W^{2})({\bf p}-{\bf q})^{2}({\bf
q}^{2}+W^{2})}+R({\bf p},{\bf q})\right]~{}~{}~{}.$ (78)

In parallel with this correction, let us also consider the contribution from
the second term (73):

$\displaystyle\psi^{(1)}_{2}(0)=-\frac{1}{4m^{3}}\int\tilde{G}_{1S}(0,{\bf
r})\hat{\bf p}^{4}\psi(r)d{\bf r})=~{}$
$\displaystyle-\frac{\pi\alpha}{m^{2}}\psi^{C}(0)\int\tilde{G}({\bf q},{\bf
p})\frac{d{\bf p}d{\bf
q}}{(2\pi)^{6}}\left[1-\frac{W^{2}(2p^{2}+W^{2})}{(p^{2}+W^{2})^{2}}\right].$
(79)

From expressions (77) , (IV.2) it follows that the divergent terms in the
matrix elements cancel each other so that the contribution of the operator
from (73) will be finite. After calculating three integrals with functions
from (78) , we obtain the following summary result:

$\displaystyle\Delta\psi^{(1)}(0)=\psi^{C}(0)[-\frac{63}{128}\alpha^{2}].$
(80)

Let us now consider one more divergent term in $\psi(0)$, which is determined
by the first part of operator (74):

$\displaystyle\Delta V^{(1)}_{2}=-\frac{\alpha}{2m^{2}r}\hat{\bf p}^{2}.$ (81)

Moving to the momentum representation, we can present the matrix element
$\Delta V^{(1)}_{2}$ in the form:

$\displaystyle\Delta\psi^{(2)}(0)=-\frac{\alpha}{2m^{2}}\int\tilde{G}_{1S}(0,r)\frac{1}{r}\hat{\bf
p}^{2}\psi(r)d{\bf r}=~{}~{}~{}~{}~{}~{}~{}~{}$ (82)
$\displaystyle=\psi^{C}(0)\frac{-2\mu\alpha^{2}}{\pi m^{2}}\int\tilde{G}({\bf
q},{\bf p})\frac{d{\bf q}}{(2\pi)^{3}}\frac{d{\bf p}}{(2\pi)^{3}}\frac{d{\bf
p^{\prime}}}{({\bf p}-{\bf p^{\prime}})^{2}}\left[\frac{1}{(p^{\prime
2}+W^{2})}-\frac{W^{2}}{(p^{\prime 2}+W^{2})^{2}}\right].~{}~{}~{}$

Isolating the divergent term in this expression, which is determined by the
first term in square brackets and the free Green’s function, we calculate it
taking into account the leading logarithmic term:

$\displaystyle\Delta\psi^{(2)}(0)=\psi^{C}(0)\frac{4W^{2}}{\pi
m^{2}}\int^{\sim m}\frac{d{\bf q}d{\bf p}}{(2\pi)^{3}(p^{2}+W^{2})({\bf
p}-{\bf p^{\prime}})^{2}(p^{\prime
2}+W^{2})}=\psi^{C}(0)\frac{1}{2}\alpha^{2}\ln\alpha^{-1}.$ (83)

It is useful to note that the same expression for the correction can be
obtained using the relativistic correction to the wave function from [24]:

$\displaystyle\Delta\psi_{rel}(k)=\psi^{C}(k)\frac{\alpha^{2}}{4}\left[\ln\frac{(k^{2}+W^{2})}{W^{2}}+\frac{(k^{2}-W^{2})}{kW}\arctan(\frac{k}{W})\right].$
(84)

Calculating the momentum integral, we obtain in leading logarithmic order:

$\displaystyle\Delta\psi_{rel}(0)=\psi^{C}(0)\frac{1}{2}\alpha^{2}\ln\alpha^{-1}.$
(85)

Other finite contributions from expression (74), corresponding to the terms of
the Coulomb Green’s function (78) have the form:

$\displaystyle\Delta\psi_{rel}^{(1)}=\psi^{C}(0)\left[-\frac{\alpha^{2}}{8}\right],~{}~{}~{}$
(86)
$\displaystyle\Delta\psi_{rel}^{(2)}=\psi^{C}(0)\frac{W^{2}}{2m^{2}\pi^{6}}\int\frac{p^{\prime
2}d{\bf p^{\prime}}d{\bf p}d{\bf q}}{({\bf p}-{\bf p^{\prime}})^{2}(p^{\prime
2}+1)^{2}(q^{2}+1)({\bf q}-{\bf p})^{2}(p^{2}+1)}=$ (87)
$\displaystyle\psi^{C}(0)\left[\frac{\alpha^{2}}{12}(\pi^{2}-\frac{3}{2})\right],$
$\displaystyle\Delta\psi_{rel}^{(3)}=\psi^{C}(0)\frac{2W^{2}}{m^{2}\pi^{4}}\int
R({\bf p},{\bf q})d{\bf p}d{\bf
q}\left[\frac{2}{p}\arctan(p)-\frac{1}{p^{2}+1}\right]=\psi^{C}(0)\left[\frac{\alpha^{2}}{12}(\frac{51}{4}-\pi^{2})\right].$
(88)

The second part of the operator $\Delta V_{2}$ also gives a finite
contribution so that we represent the full contribution of this operator in
the form:

$\displaystyle\Delta\psi_{2}=\psi^{C}(0)\left[\frac{1}{2}\alpha^{2}\ln\alpha^{-1}+\frac{5}{8}\alpha^{2}\right].$
(89)

Since we are studying the production of dimuonium and ditauonium in states
with a certain spin, we also take into account in $\psi(0)$ the spin-dependent
correction from the Breit Hamiltonian (75). Calculation of a matrix element
with a delta function taking into account the leading logarithm (the second
term in (78) gives the following result:

$\displaystyle\Delta\psi^{(3)}(0)=\psi^{C}(0)\left[2-\frac{7}{3}S(S+1)\right]\frac{1}{4}\alpha^{2}\ln\alpha^{-1}.$
(90)

The leading order correction in $\alpha$ is the correction for electron vacuum
polarization, which is determined by the potential:

$\displaystyle\Delta
V_{vp}(r)=-\frac{\alpha^{2}}{3\pi}\int_{1}^{\infty}\rho(\xi)d\xi\frac{1}{r}e^{-2m_{e}\xi
r},~{}~{}~{}\rho(\xi)=\frac{\sqrt{\xi^{2}-1}(2\xi^{2}+1)}{\xi^{4}}.$ (91)

After analytical integration over coordinates in (70), we present this
correction in integral form:

$\displaystyle\psi(0)=\psi^{C}(0)\left[a_{vp}\frac{\alpha}{\pi}\right],~{}~{}~{}r_{6}=\frac{m_{e}}{m\alpha},$
(92) $\displaystyle
a_{vp}=\int_{1}^{\infty}\frac{\rho(\xi)d\xi}{6(1+r_{6}\xi)}\left[2r_{6}^{2}\xi^{2}+7r_{6}\xi+2(1+r_{6}\xi)\ln(1+r_{6}\xi)+3\right].$

The numerical estimate of this correction for dimuonium is determined by the
electron, muon vacuum polarization, and for ditauonium by the electron, muon,
$\tau$ lepton vacuum polarization in the form:

$\displaystyle
a_{vp}(e^{+}e^{-})=0.004383,~{}~{}~{}a_{vp}(\mu^{+}\mu^{-})=0.942790,~{}~{}~{}a_{vp}(\tau^{+}\tau^{-})=3.613251.$
(93)

Summing contributions (86), (87), (88), (89), (90) and (93), we obtain a
correction to the wave function at the zero of the bound state of leptons as
follows:

$\displaystyle\psi(0)=\psi^{C}(0)\left\\{1+a_{vp}\frac{\alpha}{\pi}+\left(2-\frac{7}{6}S(S+1)\right)\frac{1}{2}\alpha^{2}\ln\alpha^{-1}-\frac{3}{128}\alpha^{2}\right\\}.$
(94)

##  V Numerical results and Conclusions

Theoretical study of rare exclusive decay processes of the Higgs boson
represents an important task in the general direction of studying the Higgs
sector in the Standard Model. In this work, the calculation of the observed
Higgs boson decay widths is carried out within the framework of the
relativistic quasipotential model, which we previously used in studying the
pair production of mesons and baryons. The creation of bound states of
particles in various reactions is also important for testing the theory of
bound states in quantum field theory.

In Table 1 we present the numerical results of calculating the decay widths of
the Higgs boson with the production of single lepton bound states and a photon
or Z-boson. The obtained values for $\Gamma$ of these rare decays are small,
and their observation is possible only at high-luminosity colliders being
designed in the future, when the production of a significant number of Higgs
bosons ($\sim 10^{10})$ will be possible [2, 3]. Taking into account the main
channel of the decay of lepton ortho-states into three photons, one can expect
their reliable detection.

In Table 1 we also present the results of pair production of bound states of
leptons, which are numerically very strongly suppressed in relation to single
production processes or similar processes with the production of charmoniums
or bottommoniums. Despite the fact that there is a mass factor $1/M^{2}$ (M is
the leptonium mass), which increases the contribution of a number of decay
amplitudes (see, for example, the amplitudes in Fig. 5, 8), nevertheless, the
second square of the modulus of the wave function of the bound state contains
the factor $M^{4}\alpha^{3}$ so that the pair production of lepton bound
states becomes completely unlikely.

Table 1: Numerical results for branching fractions Final state | Result of work [8] | Our result
---|---|---
$(e^{+}e^{-})+\gamma$ | $3.5\times 10^{-12}$ | $0.88\times 10^{-10}$
$(\mu^{+}\mu^{-})+\gamma$ | $3.5\times 10^{-12}$ | $1.12\times 10^{-11}$
$(\tau^{+}\tau^{-})+\gamma$ | $2.2\times 10^{-12}$ | $3.48\times 10^{-12}$
$(e^{+}e^{-})+Z$ | $5.2\times 10^{-13}$ | $7.87\times 10^{-13}$
$(\mu^{+}\mu^{-})+Z$ | $5.7\times 10^{-13}$ | $9.85\times 10^{-13}$
$(\tau^{+}\tau^{-})+Z$ | $1.4\times 10^{-11}$ | $5.68\times 10^{-11}$
$(e^{+}e^{-})+(e^{+}e^{-})$ | — | $2.05\times 10^{-19}$
$(\mu^{+}\mu^{-})+(\mu^{+}\mu^{-})$ | — | $1.13\times 10^{-20}$
$(\tau^{+}\tau^{-})+(\tau^{+}\tau^{-})$ | — | $1.09\times 10^{-18}$

In this work, we investigated the role of relativistic effects caused by the
law of transformation of wave functions of bound states of particles during
the transition to the rest frame, relativistic corrections in the decay
amplitude itself and corrections to the wave function of positronium,
dimuonium and ditauonium at zero. The leading corrections are of order
$\alpha$ and $\alpha^{2}\ln\alpha$. In addition, along with relativistic
corrections, we also took into account the effects of particle coupling in our
numerical results in Table 1, understanding by a purely non-relativistic
calculation result one in which the mass of the bound state of leptons is
equal to the sum of the lepton masses. The bound state effects occur because
various mass coefficients enter into the final result in large degrees. In the
case of dimuonium and ditauonium, these corrections are of order of a percent,
which is mainly due to the effect of vacuum polarization of order $O(\alpha)$
in the bound state wave function. For positronium, the leading corrections are
those of order $O(\alpha^{2})$, which leads to their significant reduction to
0.01 percent.

In calculating decay widths, we examine the contribution of various decay
mechanisms, the relative contribution of which changes in the transition from
positronium to dimuonium and ditauonium. Numerous mass parameters $r_{i}$ in
the presented analytical formulas for decay widths, which determine the final
result, change greatly during the transition between lepton bound states. But
we can definitely say that the contribution from bosonic and quark loops (see
Fig. 7,8) has the greatest values, along with the contribution of the ZZ-
mechanism in the case of pair leptonium production. The widths of rare Higgs
boson decays were calculated earlier in [8]. In general, we can say that our
results and [8] are consistent, although there are some differences in
individual decay widths.

###### Acknowledgements.

The work is supported by the Foundation for the Advancement of Theoretical
Physics and Mathematics “BASIS” (Grant No. 22-1-1-23-1).

## References

  * Aad  _et al._ [2022] G. Aad, B. Abbott, D. C. Abbott, K. Abeling, S. H. Abidi, A. Aboulhorma, and et al. [The Atlas Collaboration], A detailed map of higgs boson interactions by the atlas experiment ten years after the discovery, Nature 607, 52 (2022).
  * De Blas  _et al._ [2018] J. De Blas, R. Franceschini, F. Riva, P. Roloff, U. Schnoor, M. Spannowsky, J. Wells, A. Wulzer, J. Zupan, S. Alipour-Fard, and et al., The clic potential for new physics, CERN Yellow Rep. Monogr. 3, 1 (2018).
  * Contino  _et al._ [2017] R. Contino, D. Curtin, A. Katz, M. Mangano, G. Panico, M. Ramsey-Musolf, G. Zanderighi, C. Anastasiou, W. Astill, G. Bambhaniya, and et al., Physics at a 100 tev pp collider: Higgs and EW symmetry breaking studies, CERN Yellow Rep. 3, 255 (2017).
  * Kartvelishvili  _et al._ [2009] V. Kartvelishvili, A. V. Luchinsky, and A. A. Novoselov, Double vector quarkonia production in exclusive higgs boson decays, Phys. Rev. D 79, 114015 (2009).
  * Gao and Gong [2022] D.-N. Gao and X. Gong, Higgs boson decays into a pair of heavy vector quarkonia, Phys. Lett. B 832, 137 (2022).
  * Adam  _et al._ [2022] W. Adam, J. W. Andrejkovic, T. Bergauer, and et al. [The CMS Collaboration], Search for higgs boson decays into z and $j/\psi$ and for higgs and z boson decays into $j/\psi$ or $\upsilon$ pairs in pp collisions at $\sqrt{13}$ tev, Phys. Lett. B 842, 137534 (2022).
  * Martynenko and Martynenko [2023a] A. P. Martynenko and F. A. Martynenko, Relativistic theory of paired heavy meson and baryon production, Phys. Atom. Nucl. 86, 878 (2023a).
  * d’Enterria and Le [2024] D. d’Enterria and V. D. Le, Rare and exclusive few-body decays of the higgs, z, w bosons, and the top quark, J. Phys. G: Nucl. Part. Ph. (2024).
  * Karshenboim  _et al._ [1998] S. G. Karshenboim, V. G. Ivanov, U. D. Jentschura, and G. Soff, Bound states of the muon antimuon system: Lifetimes and hyperfine splitting, J. Exp. Theor. Phys. 86, 226 (1998).
  * Ginzburg  _et al._ [1998] I. F. Ginzburg, U. D. Jentschura, S. G. Karshenboim, F. Krauss, V. G. Serbo, and G. Soff, Production of bound mu+ mu- systems in relativistic heavy ion collisions, Phys. Rev. C 58, 3565 (1998).
  * Bogomyagkov  _et al._ [2018] A. Bogomyagkov, V. P. Druzhinin, E. Levichev, A. Milstein, and S. Sinyatkin, Low-energy electron-positron collider to search and study $(\mu^{+}\mu^{-})$ bound state, EPJ Web Conf. 181, 01032 (2018).
  * d’Enterria and Shao [2023] D. d’Enterria and H.-S. Shao, Prospects for ditauonium discovery at colliders, Phys. Lett. B 842, 137960 (2023).
  * Belov  _et al._ [2021] I. N. Belov, A. V. Berezhnoy, A. E. Dorokhov, A. K. Likhoded, A. P. Martynenko, and F. A. Martynenko, Higgs boson decay to paired bc: Relativistic and one-loop corrections, Nucl. Phys. A 1015, 122285 (2021).
  * Martynenko and Martynenko [2023b] A. P. Martynenko and F. A. Martynenko, Paired double heavy baryons production in decays of the higgs boson, Symmetry 15, 1944 (2023b).
  * Faustov  _et al._ [2023] R. N. Faustov, A. P. Martynenko, and F. A. Martynenko, Relativistic corrections to paired production of charmonium and bottomonium in decays of the higgs boson, Phys. Rev. D 107, 056002 (2023).
  * Borodylin  _et al._ [2022] V. I. Borodylin, R. N. Rogalyov, and S. R. Slabospitskii, CORE 3.2 (COmpendium of RElations, version 3.2), arXiv:1702.08246v3 [hep-ph] 1702.08246v3, 1 (2022).
  * Boradjiev  _et al._ [2018] I. Boradjiev, E. Christova, and H. Eberl, Dispersion theoretic calculation of the $h\to z+\gamma$ amplitude, Phys. Rev. D 97, 073008 (2018).
  * Melnikov  _et al._ [1994] K. Melnikov, M. Spira, and O. I. Yakovlev, Threshold effects in two-photon decays of higgs particles, Z. Phys. 64, 401 (1994).
  * Marciano  _et al._ [2012] W. J. Marciano, C. Zhang, and S. Willenbrock, Higgs decay to two photons, Phys. Rev. D 85, 013002 (2012).
  * Cutkosky [1960] R. E. Cutkosky, Singularities and discontinuities of feynman amplitudes, J. Math. Phys. 1, 429 (1960).
  * Kuipers  _et al._ [2013] J. Kuipers, T. Ueda, J. A. M. Vermaseren, and J. Vollinga, Form version 4.0., Comp. Phys. Comm. 184, 1453 (2013).
  * Khriplovich and Milstein [1994] I. B. Khriplovich and A. I. Milstein, Large relativistic corrections to the positronium decay probability, Zh. Eksp. Teor. Fiz. 106, 689 (1994).
  * Faustov and Martynenko [1999] R. N. Faustov and A. P. Martynenko, Corrections of order $(z\alpha)^{6}m_{e}^{2}/m_{\mu}$ in the muonium fine structure, J. Exp. Theor. Phys. 88, 672 (1999).
  * Lepage [1978] G. P. Lepage, Two-body bound states in quantum electrodynamics, SLAC Report: SLAC-R-212 212, 1 (1978).
  * Caswell and Lepage [1979] W. E. Caswell and G. P. Lepage, $0(\alpha\ln(\alpha^{-1}))$ corrections in positroninm: Hyperfine splitting and decay rate, Phys. Rev. A 20, 36 (1979).
  * Khriplovich and Yelkhovsky [1990] I. B. Khriplovich and A. S. Yelkhovsky, On the radiative corrections $\alpha^{2}\ln\alpha$ to the positronium decay rate, Phys. Lett. B 246, 520 (1990).
  * Kniehl and Penin [2000] B. A. Kniehl and A. A. Penin, Order $\alpha^{3}\ln(1/\alpha)$ corrections to positronium decays, Phys. Rev. Lett. 85, 1210 (2000).