\section{Introduction}
Advancements in modern deep learning have resulted in substantial progress in our ability to model, reconstruct, and forecast complex time-series data, most notably through recurrent neural networks (RNNs) and neural ordinary differential equations (neural ODEs) \cite{chen2019neuralODE,rubanova2019latentODE,rangapuram2018deep,li2021learning}. Many of these techniques have been successfully
adapted to model the spatiotemporal dynamics of complex physiological systems, such as in reconstructing \cite{jiang2021label,ghimire2018generative} or forecasting the electrical activity of the cardiac system \cite{EPNET2021}.
While strong in modeling general system dynamics, however, these methods are not designed to uncover local activity with minute signal strength compared to the global dynamics. Such local activities are common and often signify important abnormal events in physiological systems, such as an extra foci that triggers an abnormal propagation of electrical waves amid normal propagation patterns (Figure \ref{fig:normalProp}).
Due to negligible optimization loss contribution given weak signal strength, reconstructing these local activities remains challenging despite deep learning advances in time-series modelling.

We present a fundamentally novel perspective of this problem: instead of general reconstruction, we stress that these local activities, while small in signal strength, are the cause of the subsequent global activities that have larger signal strength. Therefore,
we may be able to disentangle the \textit{cause} and \textit{effect} of the latent dynamics
by explicitly modeling and inferring how the latent state of a system is influenced by potential hidden internal \textit{interventions}.
Modeling the effect on system dynamics by exogenous inputs is well-studied in classic state-space systems (SSMs), and has seen recent successes in neural-network modeling of SSMs \cite{krishnan2015deep,de2019gruODEBayes,karl2016deep,gwak2020imode}.
These works, however, focus on external interventions that are either known or directly observed, whereas we consider hidden and internal causes of system's abnormal events.
As these internal causes are not directly observable, we propose to disentangle them from the observed effects on the system's dynamics. Rather than learning the complex intervened dynamics of a system from scratch as in existing works \cite{gwak2020imode}, we incorporate knowledge of the native \textit{intervention-free} dynamics of a system and focus only on learning how it is influenced by potential hidden interventions.

To this end we present a novel neural SSM with two major innovations. First, drawing inspirations from related work in intervention modeling \cite{gwak2020imode}, we introduce causal-effect modeling of the latent dynamics via a system of interacting neural ODEs that separately describe 1) the continuous-time dynamics of the internal intervention, and 2) its effect on the trajectory of the system's native state. Second, to disentangle the cause and effect from their collective observations, we 1) leverage data of intervention-free systems to pre-train a neural ODE of \textit{native} dynamics and  integrate it into the intervention ODEs, and 2) infer the hidden intervention by assuming it to be responsible for differences observed between the actual and hypothetical \textit{intervention-free} dynamics at each time frame.

We demonstrated a proof-of-concept of the presented framework on reconstructing ectopic foci that disrupt normal electrical propagation using remote observations. We compared the presented method to the two most relevant time-series modeling approaches:
1) a global neural ODE to describe latent system dynamics,
and 2) a neural SSM with a latent neural ODE lacking intervention modeling, developed for reconstructing cardiac electrical propagation \cite{jiang2021label}.
Experiments were conducted on synthetic data with controlled \textit{internal interventions}, simulated in 2D settings.
Our results indicated that the presented method delivers more accurate inverse estimations in terms of localizing the triggering events.

\section{Methodology}

We formulate a neural SSM where the system dynamics is described on a lower-dimensional latent manifold in separation from its emission to the data space as illustrated in Fig.~\ref{fig:overview}.
Consider cardiac electrical propagation $\mathbf{x_t}$ and its body-surface measurement $\mathbf{Y_t}$ with a physics-based relation $\mathbf{Y_t} = \mathbf{H}\mathbf{x}_t$.
We enable causal-effect intervention modelling by sequentially learning two separate neural ODE functions, $\mathcal{F}_z$ and $\mathcal{F}_a$, that describe \textit{native} system dynamics and causal dynamics, respectively. The effect of intervention $\mathbf{a}$ on $\mathbf{z}$ is modelled through a coupled neural ODE function, $\mathcal{F}_{z, a}$.
Because intervention $\mathbf{a}$ is not directly observable,
it is not possible to separately learn $\mathcal{F}_z$, $\mathcal{F}_a$ and $\mathcal{F}_{z, a}$ all from scratch from only their collective observations.
Instead, we poise that we must leverage knowledge of native dynamics $\mathcal{F}_z$, in order to
be able to generate \textit{hypothetical observation}
$\hat{\mathbf{Y}}_t = \mathbf{H} \hat{\mathbf{x}}_t$
assuming intervention-free dynamics since last observation $\mathbf{Y}_{t-1}$. Then using the residual between the actual observation $\mathbf{Y}_t$ and hypothetical observation
$\hat{\mathbf{Y}}_t$ allows us to model and estimate intervention $\mathbf{a}_t$.

A two-stage optimization process is leveraged: First, the \textit{native} dynamics of the system are learned on an \textit{intervention-free} subset of data (termed \textit{ODE-VAE}) which is then statically integrated into a system of intervention ODEs optimized on the intervention set (termed \textit{ODE-VAE-IM}).

\begin{figure}[!tb]
\label{fig:overview}
    \begin{center}
        \includegraphics[width=\textwidth]{figures/overviewUpdated.png}
        \caption{A) Graphical overview of the proposed network. Inputs \textbf{H} and $\mathbf{Y}_{0:T}$ are given to the intervention dynamics $\mathcal{F}_{a}$ and influence the propagation of $\mathcal{F}_{(z,a)}$. B) Schematic of a single intervention step, showing the estimation of the intervention latent variable $\mathbf{a}_k$ from data-space observations and \textit{hypothetical intervention-free observations}.
        }
        \label{fig:overview}
    \end{center}
\end{figure}

In sections \ref{sec:naiveDynamics}-\ref{sec:dynamics} we describe the native dynamics model, the predictive model enabling causal-effect learning, the estimation of the latent intervention variable through global observations, and the model optimization loop.

\subsection{Modeling and Learning Intervention-Free Native Dynamics} \label{sec:naiveDynamics}
An ODE function $\mathcal{F}_z$ is designed to learn the \textit{native} dynamics function that handles the propagation of the system in the absence of interventions. Since we leverage this function to hypothesize the intervention-free native state of a system given at any previous system state, it must have a strong long-term forecasting ability. Therefore, inspired by \cite{yildiz2019ode2vae}, 1.) an encoding network $Enc_z$ initializes the latent ODEs vector field $\mathbf{z}_0$ using the first few frames of sequential input $\mathbf{Y}_{0:k}$ (Eq.\ref{eqn:z_ode_1}), 2.) an initial value problem is solved up to time $\mathbf{T}$ using the ODE function $\mathcal{F}_z$ (Eq.\ref{eqn:z_ode_2}), and 3.) a decoding network $Dec_z$ converts the latent trajectory $\mathbf{z}_{0:T}$ to the output space $\mathbf{X}_{0:T}$ (Eq.\ref{eqn:z_ode_3}):
\begin{align}
    \mathbf{z}_0 &= Enc_z(\mathbf{Y}_{0:k})
    \label{eqn:z_ode_1} \\
    \mathbf{z}_i &= \mathbf{z}_{i-1} + \int_{t_{i-1}}^{t_i} \mathcal{F}_{z}(\mathbf{z}_\tau) d\tau
    \label{eqn:z_ode_2} \\
    \hat{\mathbf{x}}_i &= Dec_z(\mathbf{z}_i),
    \label{eqn:z_ode_3}
\end{align}
where $Enc_z$ and $Dec_z$ are represented by CNNs and the sequence input length $k$ is empirically tuned.

We train this native dynamics model on datasets of \textit{intervention-free} systems.
After initializing the latent state, the full trajectory is predicted using $\mathcal{F}_{z}$, from which the reconstructed $\hat{\mathbf{x}}$
is compared with the ground truth in the loss:
\begin{align}
    \mathcal{L}_{native}(\mathbf{x}_{0:T}, \hat{\mathbf{x}}_{0:T}) &=
    \beta [(\mathbf{x}_0 * log(\hat{\mathbf{x}}_0)) + (1 - \hat{\mathbf{x}}_0) * log(1 - \hat{\mathbf{x}}_0)] \\
    &+ \frac{1}{T} \sum_{i=1}^T (\mathbf{x}_i * log(\hat{\mathbf{x}}_i) + (1 - \hat{\mathbf{x}_i}) * log(1 - \hat{\mathbf{x}}_i),
\end{align}
where binary cross entropy (BCE) is applied to each frame and $\beta$ represents a loss weighting coefficient. We explicitly separate out the reconstruction of $\mathbf{x}_0$ to emphasize the accuracy of the initial conditions, $\mathbf{z}_0$, independent of the accuracy of the dynamic transition models. Because the network is asked to reconstruct the full sequence of $\mathbf{x}_{0:T}$ from only an estimated initial condition $\mathbf{z}_0$, it promotes learning strong, continuous dynamics models $\mathcal{F}_{z}$ that can simulate the system given any initial condition of the space.
This function, once trained, will be incorporated into the intervention-effects function $\mathcal{F}_{(z, a)}$ as a static dynamics function whose trajectory is influenced by the intervention state.

\subsection{Modelling and Learning Intervention Dynamics} \label{sec:dynamics}

The system dynamics predictive model is described by two ODE functions, $\mathcal{F}_a$ and $\mathcal{F}_{(z, a)}$, that aim to separately estimate the cause and effect of interventions in the latent space, respectively. To achieve this, we use the state-space Bayesian filtering setting where the system propagation undergoes a series of predict-from-dynamics and update-from-observations steps \cite{karl2016deep,fraccaro2017disentangled,de2019gruODEBayes}. The former comprises of the predictive model which, given a current state $\mathbf{s}_{t-1}$, predicts the next temporal state, $\hat{\mathbf{s}}_{t}$, using the transition function $\mathcal{F}$. The latter corrects the temporal prediction through an update function, $\mathcal{G}$, using an inferred state $\mathbf{s}_{enc}(t)$ from the current observation $\mathbf{Y}_i$.

\textbf{Predictive System Dynamics Model.}
$\mathcal{F}_a$ handles the prediction of the intervention state, $a_t$, while $\mathcal{F}_{(z, a)}$ handles the intervention-influenced prediction of the latent state $z_t$, denoted by the equations in Eq.~\ref{eqn:combined}, respectively.
\begin{align}
    \mbox{\textbf{Prediction:} }
    \begin{bmatrix}
        \hat{\mathbf{a}}_i \\
        \hat{\mathbf{z}}_i
    \end{bmatrix}
    &=
    \begin{bmatrix}
        \mathbf{a}_{i-1} \\
        \mathbf{z}_{i-1}
    \end{bmatrix}
    + \int_{t_{i-1}}^{t_i}
    \begin{bmatrix}
        \mathcal{F}_{a}(\mathbf{a}_\tau) \\
        \mathcal{F}_{z}(\mathbf{z}_\tau) + \mathcal{F}_{a}(\mathbf{a}_\tau)
    \end{bmatrix}
    d\tau \label{eqn:combined}
\end{align}

\textbf{Estimating the Intervention Latent Variable.}
Different from Section 2.1 wherein our goal is to learn a native dynamic model that has strong standalone forecasting ability,
here our interest is to extract information from available observations at each time frame to uncover hidden interventions.
A dynamical function describing the temporal cause of intervention alone is not enough to produce a meaningful causal-effect model. Rather, an observational source from which the latest intervention information can be inferred is required \cite{gwak2020imode,yin2021augmentingaphyn}. We hypothesize that the latent variable associated with the interventions rises from differences observed between data space observations and \textit{hypothetical observations} continuing under \textit{intervention-free} dynamics at each timestep. For inverse image reconstruction problems, this can be formulated as the element-wise difference in the inverse solution difference of $H\mathbf{\hat{x}}^\mathbf{z}_t$ and $\mathbf{Y}_t$. This difference is then used as input to an encoder to output the estimated latent intervention state.
\begin{align}
    \mbox{\textbf{Update:} } \mathbf{a}_{enc} &= Enc_a(|| \mathbf{H}\hat{\mathbf{x}_i} - \mathbf{Y}_i ||_{i:(i+k)}) \label{eqn:pienc} \\
    \mathbf{a}_i &= \mathcal{G}_a(\hat{\mathbf{a}}_{i}, \mathbf{a}_{enc}) \label{eqn:z_ode}
\end{align}
We make use of the Gated Recurrent Unit (GRU) cell \cite{cho2014learningGRU}, $\mathcal{G}_a$ to learn a weighted combination of the encoded intervention variable and the prediction coming from the intervention dynamics function.
A number of methods may be applied to improve the forward intervention estimation. One such choice, and the one utilized in this work, is to predict a number of timesteps ahead and get the inverse solution difference over a sequence of frames rather than just the current predicted frame. This provides increased temporal information to the intervention encoder by exposing the longer-term differences that would arise from the current state.

\textbf{Optimizing \textit{Intervention} Models.}
 The loss function is the  data reconstruction loss averaged over each time frame:
\begin{align}
    \mathcal{L}_{intv}(\mathbf{X}_{0:T}, \hat{\mathbf{X}}_{0:T}) &= \frac{1}{N * T} \sum_{n=0}^N \sum_{i=0}^T (\mathbf{x}_{n,i} * log(\hat{\mathbf{x}}_{n,i}) + (1 - \hat{\mathbf{x}_{n,i}}) * log(1 - \hat{\mathbf{x}}_{n,i}),
\end{align}
where BCE is applied to each time frame. As we are no  longer care about the solution of an initial value problem, the emphasis on $\mathbf{x}_0$ reconstruction is lifted.

\section{Experiments}
For the task of reconstructing ectopic foci that disrupt normal electrical propagation using remote observations, we include three baseline methods: \textit{ECGI}.) a first order Tikhonov regularization for reconstruction, \textit{ODE-VAE}.) the direct application of the native dynamics network to the given intervention set, and \textit{ODE-VAE-GRU}.) a Bayesian filtering ablation where a latent update mechanism (GRU) is applied at every timestep using the next 3 observations straight as input. This represents a model optimized on native and intervention dynamics together without clear disentanglement or ability to preserve \textit{intervention-free} dynamics. We denote the proposed method as \textit{ODE-VAE-IM} throughout the experiments. The observation data $\textbf{X}_{0:T}$ is passed through a Bernoulli filter between activated and deactivated nodes following min-max normalization and all network decoders have a Sigmoid activation on the output images. The intervention variable encoder (\textit{$Enc_a$}) has non-activated outputs.

Initial learning rates were found using the cyclical learning rate estimation technique \cite{smith2017cyclical} and are decayed by $0.5$ at set points throughout training. We used AdamW optimizer \cite{loshchilov2017decoupled} with a weight decay of $1e-2$, a batch size of 16, and a latent dimension of 12 across all methods. All experiments were run on NVIDIA Tesla T4s with 16GB memory, taking $\sim$8 hours to train. Our implementation and saved models are available at \url{https://github.com/qu-gg/causal-effect-neural-ssm}, along with more examples in Supplementary Material.

\textbf{Data Generation.} Following the Fitzhugh-Nagumo model \cite{izhikevich2006fitzhugh}, we simulated the transmembrane potentials on a 100*100mm 2D grid. This synthetic dataset includes two subsets 1.) native transmembrane potential set containing 1000 voltage maps, where the initial excitation location is chosen randomly across the grid and 2.) transmembrane potentials in the presence of an extra focci. Simulating the setting where the extra Foci intervenes in the normal dynamics, we generated a dataset of 705 samples with varying initial excitation location and time and location of the extra Foci.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{figures/PacingSet.png}
    \caption{Synthetic TMP per-step recons in which additional excitations occur.}
    \label{fig:dualProp}
\end{figure}

\subsection{Dynamics Reconstruction}
In this section, we showcase the reconstructive and temporal capabilities of the method on normal and intervention dynamics.

\textbf{Normal Dynamics.}
We first highlight the \textit{naive} dynamics model's performance on the non-intervention datasets used for pre-training. To ensure that a robust vector field is trained, our base dataset combines both single and double excitation propagations that occur within the first few frames. We posit that by introducing more propagation patterns in the base dynamics, even if they are not used in intervention, helps the ODE function and emission model to handle unstable vector state jumps. Figure~\ref{fig:normalProp} highlights two samples of per-step reconstructions, showcasing that strong dynamics are established.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{figures/normalRecons.png}
    \caption{Reconstruction of electrical propagation in which no interventions (foci) occur.}
    \label{fig:normalProp}
\end{figure}

\textbf{Ectopic Foci.}
Figure~\ref{fig:dualProp} highlights per-timestep reconstructions of a single sample across all methods. \textit{ODE-VAE}, the base dynamics model with no means of intervention, manages to reconstruct a propagated waveform of the \textit{native} dynamics. Leveraging this fixed native dynamics function, the presented \textit{ODE-VAE-IM} is able to uncover the hidden intervention dynamics.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{figures/pacingRecons.png}
    \caption{Reconstruction of electrical propagation in which ectopic foci occurs.}
    \label{fig:dualProp}
\end{figure}

\subsection{Localization Results}
To gauge the proposed method's numeric performance, we performed a localization test in which the time and location of interventional foci were identified from the reconstruction results across a subset of 100 samples. The mean absolute error (MAE) in time frames, the Euclidean distance in location, and the percentage of identified foci are used for quantitative metrics. Table~\ref{tab:metrics} showcases these results, in which the ability of \textit{ODE-VAE-IM} in identifying the presence of abnormal foci is highlighted.

\begin{table}[!tb]
\centering
\caption{Comparison in identifying the intervention foci's location and activation timestep. Identification percentage represents cases with clear foci reconstruction.}
\label{tab:metrics}
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Model} &  \textbf{\% Foci Identified} & \textbf{Timestep Activation MAE} & \textbf{Foci Location Error} \\
\hline
ODE-VAE-IM  &  0.85 & 3.56 & 54.90 \\
ODE-VAE-GRU &  0.79 & 3.92 & 53.71 \\
ODE-VAE     &  0.39 & 3.07 & 57.28 \\
ECGI        &  0.52 & 4.10 & 52.35 \\
\hline
\end{tabular}
\end{table}

\subsection{Latent Norm Ablation}
We performed an ablation study on the latent dynamics of the cause and effect by visualizing the $L^2$-Norm of their respective latent vector states over time for two examples in Figure~\ref{fig:latentNorms}. Temporal windows are highlighted in which the intervention dynamics correctly capture the windows of starting and ending intervention activity. To facilitate comparisons between each latent space, we first perform min-max normalization over the $L^2$-Norms of each latent space to bring them within the same data range.
\begin{figure}[!tb]
    \centering
    \includegraphics[width=\textwidth]{figures/latentNorms.png}
    \caption{Visualizations of the L2-Norm of system and intervention states over time.}
    \label{fig:latentNorms}
\end{figure}

\section{Conclusion}
In this work, we propose a proof-of-concept interventional modelling framework to tackle the problem of low-strength local activity that signifies important abnormal triggering events in a dynamic system. We introduce a two-ODE system that separately models the cause and effect of system's latent state under the influence of hidden internal interventions. To this end we leverage pre-trained dynamical functions describing \textit{intervention-free} native dynamics of a system. We demonstrated the frameworks performance on reconstruction of the ectopic foci causing abnormality in cardiac electrical activity.

\textbf{Future Work.} Future work aims to extend this framework into 3-D inverse reconstruction and pacing localization using graph neural networks to enable experimentation on clinical data, following \cite{jiang2020learning,jiang2021label}. One promising direction is learning multi-interventional dynamics simultaneously via attention mechanisms.

\textbf{Limitations.} The transition between training $\mathcal{F}_{z}$ and $\mathcal{F}_{a}$ has a period of re-training up to the original performance, which identity-based initialization or loss-influence annealing $\mathcal{F}_{a}$ may alleviate. $Enc_z$ has potential generalization problems when the initial frame distributions shift between native and intervention sets and requires experimentation on its general training setup.

\bibliographystyle{splncs04}


\section{Introduction}

After the discovery of the Higgs boson \cite{atlas,cms}, a period of detailed study of the processes connected
with this particle was begun. We can say that the study of the Higgs sector has become one of the most important areas
in particle physics \cite{cms1}.
The discovery of the Higgs boson confirmed the electroweak mechanism of symmetry breaking, but the nature
of this particle remains to be explored. To do this, it is necessary to create particle colliders,
on which Higgs bosons could be produced in significant quantities, which would make it possible to proceed
to a precise study of the parameters of the Higgs sector. It is possible that interaction processes involving
the Higgs boson can provide a transition to a New Physics that lies beyond the Standard Model.

Among the parameters of the Higgs sector, the coupling constants of the Higgs boson with various bosons g(HZZ),
g(HWW), leptons $g(H\tau\tau)$, $g(H\mu\mu)$ and quarks g(Hcc), g(Hbb) stand out. They determine the decay processes of the
Higgs boson into various particles \cite{pdg}. Due to its large mass and the presence of coupling constants with different
particles, the Higgs boson has numerous decay channels. The decay channel of the Higgs boson into a pair
of heavy quarks is interesting because it creates the possibility of the production of bound states
of heavy quarks. Thus, rare exclusive decay processes of the Higgs boson into a pair of charmonium or bottomonium
are of obvious interest both for studying the decay mechanisms and for studying the properties
of bound states of quarks.

The CMS collaboration began the search for rare Higgs decays into a pair of heavy vector quarkonia in 2019 \cite{cms2}.
The results of new upper limits on the branching fractions are obtained in \cite{cms3}:
\begin{equation}
B(H\to J/\Psi,J/\Psi)<3.8\cdot 10^{-4},
\end{equation}
\begin{equation}
B(H\to \Upsilon(1S),\Upsilon(1S))<1.7\cdot 10^{-3}.
\end{equation}

Theoretical studies of the production of a pair of heavy quarkonia in the decays of the Higgs boson began
about 40 years ago in \cite{bander,keung}. In these papers, theoretical formulas were obtained for estimating
the decay widths in the nonrelativistic approximation for some decay mechanisms. Relatively more recent theoretical
studies of these processes have been carried out in \cite{luchinsky,apm2021,apm2022,neng}.

In this work we continue the study of relativistic effects in the exclusive
paired charmonium and bottomonium production in the Higgs boson decay which was begun in \cite{apm2021,apm2022}
for $B_c$ mesons.
Our calculation of the decay widths
is performed on the basis of relativistic quark model used previously for the investigation
of relativistic corrections in different reactions in \cite{apm1,apm2,apm3}.
Despite the rare nature of the Higgs boson decays being studied, one can hope that such processes can be
investigated at the new Higgs boson factories.
The aim of the present work is to give new analysis of pair quarkonium production in Higgs boson decays in the
Standard Model:

1. In comparison with previous works \cite{luchinsky,neng} we consider different production mechanisms
(quark-gluon, quark-photon, photon-photon and Z-boson mechanisms) of the pair charmonium and
bottomonium production.

2. We take into account relativistic corrections connected with the relative motion
of heavy quarks both in the pair production amplitude itself and in the wave functions
of the bound states.

As a result, new numerical estimation for the Higgs boson decay rates is obtained.
An experimental study of these rare decays of the Higgs boson at future high-luminosity accelerators could
be useful in testing the Standard Model with greater accuracy.

\section{General formalism}

From the very beginning, it should be emphasized that there are several groups of amplitudes for the decay
of the Higgs boson with the formation of a pair of charmoniums $J/\Psi$ (bottomonium $\Upsilon$), which contribute
the same order of magnitude to the decay width. We study the amplitudes which are shown in Fig.~\ref{pic1}-Fig.~\ref{pic5}
and representing different decay mechanisms with pair vector meson production:
quark-gluon, quark-photon, quark-loop, W-boson loop and Z-boson.
The first group includes quark-gluon amplitudes shown in Fig.~\ref{pic1} (quark-gluon mechanism).
The factor determining the order of contributions can be represented as $\alpha_s/M_H^4$, where $M_H$ is
the mass of the Higgs boson. Such a factor can be distinguished from the very beginning due to the structure
of the interaction vertices and the denominators of the particle propagators.
The second group is formed by quark-photon amplitudes shown in Fig.~\ref{pic2} (quark-photon mechanism).
The order of contribution is determined here
by the following factor $\alpha/M_H^2 M_{{\cal Q\bar Q}}^2$. The third group is formed by the amplitudes
in Fig.~\ref{pic3}-Fig.~\ref{pic4} which contain the quark or W-boson loop with two photons
that create a pair of $J/\Psi$ (or $\Upsilon$) mesons (quark-loop and W-boson loop mechanisms).
The primary common factor in this case takes the form $\alpha^2/M_{{\cal Q\bar Q}}^4$.
Finally, the last group of amplitudes in Fig.~\ref{pic5} is determined by the interaction of the Higgs boson
with a pair of Z-bosons (Z-boson mechanism).
To achieve good calculation accuracy, it is necessary to take into account the contribution
of all amplitudes from these groups in Fig.~\ref{pic1}-Fig.~\ref{pic5}.

Let consider firstly the Higgs boson decay amplitudes shown in Fig.~\ref{pic1}, \ref{pic2}. For pair production
of quarkonia in the leading order of perturbation theory, it is necessary to obtain at the first stage two
free quarks and two free antiquarks. Then they can form bound states with some probability at the next stage.
In the quasipotential approach the decay amplitude
can be presented as a convolution of a perturbative production amplitude
of two $c$-quark and $\bar c$-antiquark pairs and the quasipotential
wave functions of the final mesons \cite{apm2021,apm2022}:
\begin{equation}
\label{eq1}
{\cal M}(P,Q)=-i(\sqrt{2}G_F)^{\frac{1}{2}}\frac{2\pi}{3}M_{{\cal{Q\bar Q}}}
\int\frac{d{\bf p}}{(2\pi)^3}\int\frac{d{\bf q}}{(2\pi)^3}\times
\end{equation}
\begin{displaymath}
\times Tr\left\{\Psi^{\cal V}(p,P)\Gamma_1^{\nu}(p,q,P,Q)\Psi^{\cal V}(q,Q)
\gamma_\nu+
\Psi^{\cal V}(q,Q)\Gamma_2^{\nu}(p,q,P,Q)\Psi^{\cal V}(p,P)
\gamma_\nu\right\},
\end{displaymath}
where $M_{{\cal{Q\bar Q}}}$ is the mass of quarkonium,
Four-momenta $p_1$ and $p_2$ of $c$-quark and $\bar c$-antiquark in the pair
forming the first ${\cal (Q\bar Q)}$ meson, and four-momenta $q_2$ and $q_1$ for
quark and antiquark in the second meson ${\cal (Q\bar Q)}$ are expressed in terms of relative and total
four momenta as follows:
\begin{equation}
\label{eq2}
p_{1,2}=\frac12 P \pm p,\quad (pP)=0; \qquad q_{1,2}=\frac12 Q \pm q,\quad (qQ)=0,
\end{equation}
A superscript ${\cal P}$ indicates a pseudoscalar meson ${\cal (Q\bar Q)}$, a superscript ${\cal V}$
indicates a vector meson ${\cal (Q\bar Q)}$.
The vertex functions $\Gamma_{1,2}$ are presented below in leading order and in Fig.~\ref{pic1}.
Heavy quarks $c$, $b$ and antiquarks $\bar c$, $ \bar b$ are outside
the mass shell in the intermediate state:
$p_{1,2}^2\not= m^2$, so that $p_1^2-m^2=p_2^2-m^2$,
which means that there is a symmetrical exit of particles from the mass shell.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{pic1c.eps}
\caption{Quark-gluon mechanism of the pair charmonium production
in the Higgs boson decay.
Dashed line shows the Higgs boson and wavy line corresponds to the gluon.}
\label{pic1}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{pic2.eps}
\caption{Quark-photon mechanism of the pair charmonium production
in the Higgs boson decay.
Dashed line shows the Higgs boson and wavy line corresponds to the photon.}
\label{pic2}
\end{figure}

In Eq.~\eqref{eq1} we integrate over the relative three-momenta of quarks and antiquarks
in the final state.
The systematic account of all terms depending on the relative quark momenta $p$ and $q$
in the decay amplitude is
important for increasing the accuracy of the calculation.
$p=L_P(0,\mathbf p)$ and $q=L_Q(0,\mathbf q)$
are the relative four-momenta obtained by the Lorentz transformation of four-vectors
$(0,\mathbf p)$ and $(0,\mathbf q)$ to the reference frames moving with the four-momenta $P$ and $Q$.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{pic3.eps}
\caption{Quark loop mechanism of the pair charmonium production
in the Higgs boson decay.
Dashed line shows the Higgs boson and wavy line corresponds to the photon.}
\label{pic3}
\end{figure}

The relativistic wave functions of the bound quarks accounting for the transformation from the
rest frame to the moving one with four momenta $P$, and $Q$, are
\begin{equation}
\label{eq3}
\Psi(p,P)=\frac{\Psi_0(\mathbf p)}{\bigl[\frac{\epsilon(p)}{m}
\frac{\epsilon(p)+m}{2m}\bigr]}
\left[
\frac{\hat v_1-1}{2}+\hat v_1\frac{\mathbf p^2}{2m(\epsilon(p)+m)}-\frac{\hat p}{2m}
\right] \times \\
\hat\varepsilon(P,S_z)
\end{equation}
\begin{displaymath}
(1+\hat v_1)
\left[
\frac{\hat v_1+1}{2}+\hat v_1\frac{\mathbf p^2}{2m(\epsilon(p)+m)}+\frac{\hat p}{2m}
\right],
\end{displaymath}
\begin{equation}
\label{eq4}
\Psi(q,Q)=\frac{\Psi_0(\mathbf q)}{\bigl[\frac{\epsilon(q)}{m}
\frac{\epsilon(q)+m}{2m}\bigr]}
\left[\frac{\hat v_2-1}{2}+\hat v_2\frac{\mathbf q^2}{2m(\epsilon(q)+m)}+\frac{\hat q}{2m}
\right] \times \\
\hat\varepsilon(Q,S_z)
\end{equation}
\begin{displaymath}
(1+\hat v_2)
\left[
\frac{\hat v_2+1}{2}+\hat v_2\frac{\mathbf q^2}{2m(\epsilon(q)+m)}-\frac{\hat q}{2m}
\right],
\end{displaymath}
where the hat symbol means a contraction of the four--vector with the Dirac gamma matrices;
$v_1=P/M_{{\cal{Q\bar Q}}}$, $v_2=Q/M_{{\cal{Q\bar Q}}}$; $\epsilon(p)=\sqrt{m^2+\mathbf p^2}$,
$m$ is $c(b)$-quark mass,
and $M_{{\cal{Q\bar Q}}}$ is the mass of charmonium (bottomonium) state.
$\varepsilon^\lambda(P,S_z)$ is the polarization vector of the $J/\Psi(\Upsilon)$ meson.

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{pic41.eps}
\includegraphics[scale=0.3]{pic42.eps}
\caption{W-boson loop mechanism of the pair charmonium production
in the Higgs boson decay.
Dashed line shows the Higgs boson and wavy line corresponds to the photon.}
\label{pic4}
\end{figure}

Expressions \eqref{eq3} and \eqref{eq4} represent complicated functions depending on relative
momenta ${\bf p}$, ${\bf q}$ including the bound state wave function in the rest frame $\Psi_0({\bf p}) $.
The color part of the meson wave function in the amplitudes \eqref{eq3}-\eqref{eq4} is taken as
$\delta_{ij}/\sqrt{3}$ (color indexes $i, j, k=1, 2, 3$).

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{pic5.eps}
\caption{Z-boson mechanism of the pair charmonium production
in the Higgs boson decay.
Dashed line shows the Higgs boson and wavy line corresponds to the Z-boson.}
\label{pic5}
\end{figure}

The general structure of expressions \eqref{eq3}-\eqref{eq4} allows us to say that they are the product of the wave
functions of mesons in the rest frame and special projection operators resulting from the transformation
from the moving reference frame to the reference frame in which the meson is at rest.
Expressions \eqref{eq3}-\eqref{eq4} make it possible to correctly take into account the relativistic
corrections connected with the relative momenta of quarks in the final states.
It is useful to note that the expression for the projection operators was obtained in the framework
of nonrelativistic quantum chromodynamics in \cite{bodwin2002} in a slightly different form for the
case when the quark momenta lie on the mass shell.
The transformation law of the bound state wave functions of quarks, which is used in the derivation
of equations \eqref{eq3}-\eqref{eq4} , was obtained in the Bethe-Salpeter
approach in \cite{brodsky} and in quasipotential method in \cite{faustov}.

As follows from \eqref{eq3}-\eqref{eq4}, when constructing the decay amplitudes of the Higgs boson, projection
operators are introduced for quark-antiquark pairs onto the spin states with total spin $S=1$
of the following form:
\begin{equation}
\label{eq5}
\hat\Pi^{\cal V}=[v_2(0)\bar u_1(0)]_{S=1}=\hat\varepsilon\frac{1+\gamma^0}{2\sqrt{2}}.
\end{equation}

Total amplitude of the Higgs boson decay to paired vector quarkonium in the case of quark-gluon mechanism
in the leading order in strong coupling constant $\alpha_s$ can be presented in the form:
\begin{equation}
\label{eq6}
{\mathcal M}=\frac{4\pi}{3}M_{{\cal{Q\bar Q}}}\alpha_s \Gamma_Q
\int\!\frac{d\mathbf{p}}{(2\pi)^3}\int\!\frac{d\mathbf{q}}{(2\pi)^3}
\mathrm{Tr}\bigl\{\mathcal T_{12}+\mathcal T_{34}\bigr\},
\end{equation}
\begin{equation}
\label{eq7}
{\mathcal T_{12}}=\Psi^{\cal V}(p,P)
\left[\frac{\hat p_1-\hat r+m}{(r-p_1)^2-m^2}\,\gamma_\mu+\gamma_\mu\,
\frac{\hat r-\hat q_1+m}{(r-q_1)^2-m^2}\right]D^{\mu\nu}(k_2)
\Psi^{\cal V}(q,Q)\gamma_\nu,
\end{equation}
\begin{equation}
\label{eq8}
{\mathcal T_{34}}=\Psi^{\cal V}(q,Q)
\left[\frac{\hat p_2-\hat r+m}{(r-p_2)^2-m^2}\,\gamma_\mu+
\gamma_\mu\frac{\hat r-\hat q_2+m}{(r-q_2)^2-m^2}\,\right]D^{\mu\nu}(k_1)
\Psi^{\cal V}(p,P)
\gamma_\nu,
\end{equation}
where $\alpha_{s}=\alpha_s\left(\frac{M_H^2}{4\Lambda^2}\right)$,
$\Gamma_Q=m(\sqrt{2} G_F)^{\frac{1}{2}}$.

Four-momentum of Higgs boson squared $r^2=M_H^2=(P+Q)^2=2M_{{\cal{Q\bar Q}}}^2+2PQ$,
the gluon four-momenta are $k_1=p_1+q_1$, $k_2=p_2+q_2$.
Relative momenta $p$, $q$ of heavy quarks enter in the gluon propagators $D_{\mu\nu}(k_{1,2})$
and quark propagators as well as in relativistic wave functions~\eqref{eq3} and \eqref{eq4}.
Accounting for the small ratio of relative quark momenta $p$ and $q$ to the mass of the Higgs boson $M_H$,
we can simplify the inverse denominators of quark and gluon propagators as follows:
\begin{equation}
\label{eq9}
\frac{1}{(p_1+q_1)^2}\approx \frac{1}{(p_2+q_2)^2}=\frac{4}{M_H^2},
\end{equation}
\begin{equation}
\label{eq10}
\frac{1}{(r-q_1)^2-m_1^2}=\frac{1}{(-r-p_1)^2-m_1^2}=\frac{1}{(r-p_2)^2-m_1^2}=\frac{1}{(-r-q_2)^2-m_1^2}
=\frac{2}{M_H^2}.
\end{equation}
In \eqref{eq9}-\eqref{eq10} we completely neglect corrections of the form $|{\bf p}|/M_H$, $|{\bf q}|/M_H$.
At the same time, we keep in the amplitudes \eqref{eq7}, \eqref{eq8} the second-order correction
for small ratios $|{\bf p}|/m$, $|{\bf q}|/m$ relative to the leading order result.
Calculating the trace in obtained expression
in the package FORM \cite{form}, we find relativistic amplitudes of the paired meson production
in the form:
\begin{equation}
\label{eq12}
{\cal M_{VV}}^{(1)}=\frac{256\pi}{3M_H^4}(\sqrt{2}G_F)^{\frac{1}{2}}m M_{{\cal{Q\bar Q}}} \alpha_s
\varepsilon_1^{\lambda}\varepsilon_2^{\sigma}
F_{1,{\cal VV}}^{\lambda\sigma}|\tilde\Psi_{\cal V}(0)|^2,
\end{equation}
where $\varepsilon_1^\lambda$, $\varepsilon_2^\sigma$ are the polarization vectors of spin 1
mesons. The superscript in amplitude designation and subscript in tensor function $F_{{\cal VV}}$ designation
denotes the contribution of the quark-gluon mechanism in Fig.~\ref{pic1}.

The contribution of the amplitudes in Fig.~\ref{pic2} must also be taken into account, since these
amplitudes have the same order as the previous ones, despite the replacement $\alpha_s\to \alpha$.
This is due to the presence in the denominator of the mass of the meson instead of the mass
of the Higgs boson \cite{bodwin2014}.
The expression for the production amplitudes of the pair $J/\Psi$ has a similar structure with slight
changes in the common factors:
\begin{equation}
\label{eq13}
{\cal M_{VV}}^{(2)}=\frac{288\pi}{M_H^2M_{\cal{Q\bar Q}}}(\sqrt{2}G_F)^{\frac{1}{2}}m e_Q^2\alpha
\varepsilon_1^{\lambda}\varepsilon_2^{\sigma}
F_{2,{\cal VV}}^{\lambda\sigma}|\tilde\Psi_{\cal V}(0)|^2.
\end{equation}

The tensor corresponding to the quark or W-boson loops in Fig.~\ref{pic3}-Fig.~\ref{pic4}
has the structure (the subscript denotes the contribution of the quark or bosonic loop):
\begin{equation}
\label{eq13d}
T^{\mu\nu}_{Q,W}=A_{Q,W}(t)(g^{\mu\nu}(v_1v_2)-v_1^\nu v_2^\mu)+B_{Q,W}(t)[v_2^\mu-v_1^\mu (v_1v_2)][v_1^\nu-v_2^\nu(v_1v_2)],
\end{equation}
where  $t=\frac{M_h^2}{4m_Q^2}$ or $t=\frac{M_h^2}{4m_W^2}$. The structure functions $A_{Q,W}(t)$, $B_{Q,W}(t)$ can be
obtained using an explicit expression for a loop integrals (see Appendix A,B).

The amplitudes in Fig.~\ref{pic3}-\ref{pic5} with quark, W-boson loops and $ZZ$ in
an intermediate state contain the contributions of direct and crossed diagrams.
In this case, direct diagrams, in which virtual photons or Z-bosons give vector quarkonia in the final
state, are dominant in terms of the mass factor.
However, the structure of the numerators of the direct and cross amplitudes is different, which can
eventually lead to numerically close contributions.
In what follows, we present the expressions
for these amplitudes only in the leading order:
\begin{equation}
\label{eq13a}
{\cal M_{VV}}^{(3)}=\frac{2052\pi^2}{m_Q M_{\cal{Q\bar Q}}}(\sqrt{2}G_F)^{\frac{1}{2}}e_q^2 e_Q^2\alpha^2
\varepsilon_1^{\lambda}\varepsilon_2^{\sigma}
F_{3,{\cal VV}}^{\lambda\sigma}|\tilde\Psi_{\cal V}(0)|^2,
\end{equation}
\begin{equation}
\label{eq13b}
{\cal M_{VV}}^{(4)}=\frac{48\pi^2 M_Z M_W}{M^4_{\cal{Q\bar Q}}}(\sqrt{2}G_F)^{\frac{1}{2}} e_q^2\alpha^2\cos\theta_W
\varepsilon_1^{\lambda}\varepsilon_2^{\sigma}
F_{4,{\cal VV}}^{\lambda\sigma}|\tilde\Psi_{\cal V}(0)|^2,
\end{equation}
\begin{equation}
\label{eq13c}
{\cal M_{VV}}^{(5)}=\frac{48\pi\alpha}{M_Z^2\sin^2 2\theta_W}(\sqrt{2}G_F)^{\frac{1}{2}}
\varepsilon_1^{\lambda}\varepsilon_2^{\sigma}
F_{5,{\cal VV}}^{\lambda\sigma}|\tilde\Psi_{\cal V}(0)|^2,
\end{equation}
where $m$ is the mass of heavy quark c, b, produced in the vertex of Higgs boson decay, $m_Q$ is
the heavy quark mass in quark loop, $e_q$ is the charge (in units e) of heavy quark (c or b) entering
in final mesons, $e_Q$ is the charge (in units e) of heavy quark (c or b) in the quark loop.

The tensor functions in amplitudes \eqref{eq12}-\eqref{eq13c} have the following form:
\begin{equation}
\label{eq15}
F^{\alpha\beta}_{i,{\cal VV}}=g^{(i)}_1 v_1^\alpha v_2^\beta+g^{(i)}_2g^{\alpha\beta},~~~
g^{(1)}_1=-2+\frac{2}{9}\omega_{1}^2,
\end{equation}
\begin{displaymath}
g^{(1)}_2=-1-2r_2+r_1^2+\frac{4}{3}r_2\omega_{1}+\frac{1}{9}\omega_{1}^2+\frac{2}{3}r_2\omega_1^2
-\frac{1}{9}r_1^2\omega_1^2,
\end{displaymath}
\begin{equation}
\label{eq16}
g^{(2)}_1=4-\frac{4}{9}\omega_{1}^2,~
g^{(2)}_2=2+4r_2-2r_1^2-\frac{8}{3}r_2\omega_1-\frac{2}{9}\omega_1^2-\frac{4}{3}r\omega_1^2+\frac{2}{9}r_1^2\omega_{1}^2,
\end{equation}
\begin{equation}
\label{eq16a}
g^{(3)}_{1,Q}= -A_Q(t)(1 + \frac{2}{3}\omega_1 + \frac{1}{9}\omega_1^2)+
B_Q(t)(1 + \frac{2}{3}\omega_1 + \frac{1}{9}\omega_1^2),~
\end{equation}
\begin{displaymath}
g^{(3)}_{2,Q} =A_Q(t)\left(-1 - \frac{2}{3}\omega_1 - \frac{1}{9}\omega_1^2 +\frac{1}{2}r_1^2+\frac{1}{3}\omega_1r_1^2+
\frac{1}{18}\omega_1^2r_1^2\right),
\end{displaymath}
\begin{equation}
g^{(4)}_1=-A_W(t)\left(1+\frac{2}{3}\omega_1+\frac{1}{9}\omega_1^2\right)+B_W(t)\left(1+\frac{2}{3}\omega_1+
\frac{1}{9}\omega_1^2\right),
\end{equation}
\begin{equation}
\label{eq16b}
g^{(4)}_2=A_W(t)\left(-1-\frac{2}{3}\omega_1-\frac{1}{9}\omega_1^2+\frac{1}{2}r_1^2+\frac{1}{3}\omega_1r_1^2+
\frac{1}{18}\omega_1^2r_1^2\right),~
\end{equation}
\begin{equation}
\label{eq16c}
g^{(5)}_{2}=\left(1+\frac{1}{3}\omega_1\right)^2
\left(\frac{1}{2}-a_z\right)^2-\frac{M_z^4}{3\left(\frac{M_h^2}{4}-M_Z^2
\right)^2}\Bigl(-\frac{1}{4}-\frac{1}{6}\omega_1-\frac{1}{36}\omega_1^2+
\end{equation}
\begin{displaymath}
+\frac{1}{2}a_z+\frac{1}{3}\omega_1a_z+\frac{1}{18}\omega_1^2a_z-\frac{1}{2}a_z^2-\frac{1}{3}\omega_1a_z^2-\frac{1}{18}
\omega_1^2a_z^2\Bigr),~~~g^{(5)}_{1}=0,
\end{displaymath}
where the parameter $r_1=\frac{M_H}{M_{{\cal Q\bar Q}}}$, $r_2=\frac{m}{M_{{\cal Q\bar Q}}}$,
$a_z=2|e_Q|\sin^2\theta_W$.

The decay widths of the Higgs boson into a pair of vector quarkonia states are
determined by the following expressions (see also \cite{apm2021,apm2022}):
\begin{equation}
\label{eq18}
\Gamma_{\cal{VV}}=\frac{2^{14}\sqrt{2}\pi\alpha_s^2 m^2G_F  |\tilde\Psi_{\cal V}(0)|^4
\sqrt{\frac{r_1^2}{4}-1}}{9M_H^5 r_1^5}\sum_{\lambda,\sigma}
\vert\varepsilon_1^\lambda\varepsilon_2^\sigma F^{\lambda\sigma}_{{\cal VV}}\vert^2,
\end{equation}
\begin{equation}
\label{eq18a}
F^{\lambda\sigma}_{{\cal VV}}=\Bigl[g^{(1)}_1+\frac{9}{16}r_1^2\frac{e_q^2\alpha}{\alpha_s}g^{(2)}_1+
\sum_Q \frac{27\pi}{8}r_1^4\frac{e_Q^2e_q^2\alpha^2m^2_Q}{\alpha_s mM_{{\cal Q\bar Q}}}g^{(3)}_{1,Q}+
\frac{9\pi e_q^2\alpha^2 r_1^4 M_Z M_W}{64\alpha_s m M_{\cal Q\bar Q}}g^{(4)}_1+
\end{equation}
\begin{displaymath}
\frac{9M_H^4\alpha}{16 M_Z^2m M_{\cal Q\bar Q}\alpha_s}\frac{\left(\frac{1}{2}-2|e_q|\sin^2\theta_W\right)^2}{
\sin^22\theta_W}g^{(5)}_{1}
\Bigr]v_1^\sigma v_2^\lambda+
\Bigl[g^{(1)}_2+\frac{9}{16}r_1^2\frac{e_q^2\alpha}{\alpha_s}g^{(2)}_2+
\sum_Q\frac{27\pi}{8}r_1^4\frac{e_q^2e_Q^2\alpha^2 m^2_Q}{\alpha_s m M_{{\cal Q\bar Q}}}g^{(3)}_{2,Q}+
\end{displaymath}
\begin{displaymath}
\frac{9\pi e_q^2\alpha^2 r_1^4 M_Z M_W}{64\alpha_s m M_{\cal Q\bar Q}}\cos\theta_W g^{(4)}_2+
\frac{9M_H^4\alpha}{16 M_Z^2m M_{\cal Q\bar Q}\alpha_s}\frac{1}{
\sin^22\theta_W}g^{(5)}_{2}\Bigr]g^{\lambda\sigma}.
\end{displaymath}
We found it convenient to separate in square brackets in \eqref{eq18a} the coefficients denoting the relative
contribution of different decay mechanisms with respect to the quark-gluon mechanism in Fig.~\ref{pic1}.
The common factor in \eqref{eq18} corresponds to the amplitudes of the quark-gluon decay mechanism.

\begin{table}[htbp]
\caption{Numerical values of the relativistic parameters~\eqref{eq18}.}
\bigskip
\label{tb1}
\begin{ruledtabular}
\begin{tabular}{|c|c|c|c|c|c|}
$n^{2S+1}L_J$ & Meson &$M^{exp}_{{\cal Q\bar Q}}$, GeV &$\tilde R(0),~~~GeV^{3/2}$ & $\omega_{1}$ & $\omega_{2}$ \\ \hline
$1^1S_0$    & $\eta_c$     &  2.9839  & 0.92 &  0.20    &   0.0087      \\    \hline
$1^3S_1$    & $J/\Psi$     &  3.0969  & 0.81 &  0.20    &   0.0078      \\    \hline
$1^1S_0$    & $\eta_b$     &  9.3987  & 1.95 &  0.05    &   0.0044      \\    \hline
$1^3S_1$    & $\Upsilon$   &  9.4603  & 1.88 &  0.05    &   0.0044      \\    \hline
\end{tabular}
\end{ruledtabular}
\end{table}

General expression for the decay rate \eqref{eq18} contains numerous parameters.
One part of the parameters, such as quark masses, the masses of mesons are determined
within the framework of quark models as a result of calculating the observed quantities.
The parameters of quark models are found from the condition of the best agreement with
experimental data \cite{pdg}.
Another part of the relativistic parameters can also be found in the quark model as
a result of calculating integrals with wave functions of quark bound states in the
momentum representation.
Thus, we can say that the approach to calculation based on the relativistic quark model is closed, since
it allows calculating all the necessary parameters without resorting to additional hypotheses.

The amplitudes of the pair production of vector charmonium and bottomonium in the decay of the Higgs boson
are expressed in terms of function $F_{VV}$,
that is presented in the form of an expansion in ${|\bf p}|/m$,
${|\bf q}|/m$ up to terms of the second order.
As a result of algebraic transformations, it turns out to be convenient to express
relativistic corrections in terms of relativistic parameters
$\omega_{n}$. In the case of S-states $\omega_{n}$ are determined by
the momentum integrals $I_{n}$ in the form \cite{apm4,apm5}:
\begin{equation}
\label{eq19}
I_{n}^{P,V}=\int_0^\infty p^2R^{P,V}(p)\frac{(\epsilon(p)+m)}
{2\epsilon(p)}
\left(\frac{\epsilon(p)-m}{\epsilon(p)+m}\right)^n dp,
\end{equation}
\begin{equation}
\label{eq20}
\tilde R(0)=\frac{\sqrt{2}}{\sqrt{\pi}}\int_0^\infty \frac{(\epsilon(p)+m)}
{2\epsilon(p)}p^2R(p)dp,~~
\omega^{P,V}_{1}=\frac{I^{P,V}_{1}}{I^{P,V}_{0}},~~
\omega^{P,V}_{2}=\frac{I^{P,V}_{2}}{I^{P,V}_{0}},
\end{equation}
where superscripts $P,V$ denote pseudoscalar and vector states.

\begin{table}[htbp]
\caption{Numerical results for the decay widths in the nonrelativistic approximation
and with the account for relativistic corrections.}
\bigskip
\label{tb2}
\begin{ruledtabular}
\begin{tabular}{|c|c|c|}
Final state & Nonrelativistic decay width   & Relativistic decay width  \\
$(Q\bar Q)(Q\bar Q)$    & $\Gamma_{nr}$    in GeV  &$\Gamma_{rel}$  in GeV \\  \hline
$J/\Psi$ +$J/\Psi$      & $3.29 \cdot 10^{-12}$ &  $0.69 \cdot 10^{-12}$ \\
$1^3S_1+1^3S_1$         &          &           \\  \hline
$\Upsilon$ +$\Upsilon$  & $0.63 \cdot 10^{-12}$ & $0.74\cdot 10^{-12}$ \\
$1^3S_1+1^3S_1$         &          &           \\  \hline
\end{tabular}
\end{ruledtabular}
\end{table}

Another source of relativistic corrections is related with the Hamiltonian of the heavy quark
bound states which allows to calculate
the bound state wave functions of pseudoscalar and vector mesons (S-states).
The exact form of the bound state wave functions $\Psi_0({\bf q})$
is important to obtain more reliable predictions for the decay widths.
In the nonrelativistic approximation the Higgs boson decay width with a production of a pair
of quarkonium contains the fourth power of the nonrelativistic wave function at the origin
for S-states.
The value of the decay width is very sensitive
to small changes of $\tilde R(0)$. In the nonrelativistic QCD there
exists corresponding problem of determining the magnitude of the
color-singlet matrix elements \cite{bbl}. To account for relativistic corrections to the meson
wave functions we describe the dynamics of heavy quarks by the QCD generalization of
the standard Breit Hamiltonian in the center-of-mass reference frame
\cite{repko1,pot1,godfrey,lucha1995,rqm2,rqm3,repko2}.

\begin{table}[htbp]
\caption{Numerical values of parameters A and B for W-boson and quark loop mechanisms.}
\bigskip
\label{tb22}
\begin{ruledtabular}
\begin{tabular}{|c|c|c|}
Parameter  & $H\to J/\psi\:J/\psi$ &$H\to \Upsilon\:\Upsilon$    \\ \hline
$A_c$  & $-2.74\cdot 10^{-4}+1.07\cdot 10^{-4}i$ &  $-3.39\cdot 10^{-3}+0.99\cdot 10^{-3}i$  \\ \hline
$A_b$  & $-5.68\cdot 10^{-5}+7.87\cdot 10^{-4}i$ &  $-1.53\cdot 10^{-3}+0.74\cdot 10^{-3}i$  \\ \hline
$A_t$  & $7.02\cdot 10^{-7}$ &  $6.54\cdot 10^{-6}$  \\ \hline
$A_W$ & $7.94\cdot 10^{-5}$ &   $7.40\cdot 10^{-4}$  \\ \hline
$B_c$  & $0.31\cdot 10^{-4}+0.62\cdot 10^{-10}i$ &  $0.13\cdot 10^{-3}+0.51\cdot 10^{-7}i$  \\ \hline
$B_b$  & $-0.42\cdot 10^{-7}+0.42\cdot 10^{-10}i$ &  $0.37\cdot 10^{-3}+0.35\cdot 10^{-7}i$   \\ \hline
$B_t$ & $0.12\cdot 10^{-14}$ &  $-0.13\cdot 10^{-4}$   \\ \hline
$B_W$ & $1.42\cdot 10^{-11}$ &   $2.76\cdot 10^{-4}$  \\ \hline
\end{tabular}
\end{ruledtabular}
\end{table}

\begin{table}[htbp]
\caption{The contributions of different mechanisms to the Higgs boson decay widths in GeV.}
\bigskip
\label{tb3}
\begin{ruledtabular}
\begin{tabular}{|c|c|c|}
\multicolumn{3}{|c|}{The contribution accounting for relativistic corrections} \\ \hline
Contribution  & $H\to J/\psi\:J/\psi$ &$H\to \Upsilon\:\Upsilon$    \\ \hline
Fig.\ref{pic1}  & $0.36\cdot 10^{-15}$ &  $0.10\cdot 10^{-12}$  \\ \hline
Fig.\ref{pic2}  & $0.80\cdot 10^{-12}$ &  $0.16\cdot 10^{-12}$   \\ \hline
Fig.\ref{pic3} & $0.70\cdot 10^{-13}$ &  $0.37\cdot 10^{-12}$   \\ \hline
Fig.\ref{pic4} & $0.74\cdot 10^{-13}$ &   $0.68\cdot 10^{-13}$  \\ \hline
Fig.\ref{pic5} & $0.22\cdot 10^{-12}$ &  $1.45\cdot 10^{-12}$   \\ \hline
Total contribution & $0.69\cdot 10^{-12}$   &  $0.74\cdot 10^{-12}$    \\ \hline
\end{tabular}
\end{ruledtabular}
\end{table}

Using the effective Hamiltonian of quarkonia constructed on the basis of perturbative quantum chromodynamics
with allowance for nonperturbative effects, a model of heavy quarkonium for S-states is obtained. Its main
elements are presented in our previous works \cite{apm2022,apm1,apm2,apm4,apm5}. It was used to calculate
those relativistic parameters that
are included in the decay width of the Higgs boson \eqref{eq18}. The numerical values
of  the parameters $\omega_1$, $\omega_2$ are presented in the Table~\ref{tb1}. The parameter $\omega_{2}$
is not included in the decay width, since corrections of the order $O({\bf q}^4)$, $O({\bf p}^4)$ connected with it,
are omitted.

The results of numerical calculations of the decay widths and parameters A and B
are presented in Tables~\ref{tb2}, \ref{tb22}, \ref{tb3}.
Comparing the contributions of different meson pair production mechanisms, it is important
to emphasize that their relative magnitude in the total result is determined by such important parameters
as $\alpha$ and $\alpha_s$ and the particle mass ratio.
So, for example, in the case of the
production of a pair of vector mesons, the decay amplitude has two contributions from the gluon (Fig.~\ref{pic1})
and photon (Fig.~\ref{pic2}) diagrams. On the one hand, the contribution of photon amplitudes is proportional to $\alpha$,
which leads to its decrease in comparison with the contribution from gluon amplitudes. But on the other hand,
the photon contribution contains an additional factor $e_Q^2r_1^2$, which leads
to an increase in the decay widths.

\section{Numerical results and conclusion}

The study of rare exclusive decay processes of the Higgs boson is an important task, which makes it possible
to refine the values of the interaction parameters of particles in the Higgs sector. In this paper, we have attempted
to present a complete study of decay processes $H\to J/\Psi J/\Psi$, $H\to \Upsilon \Upsilon$
in the Standard Model by considering the various decay mechanisms shown
in Fig.\ref{pic1}-Fig.\ref{pic5}. To improve the calculation accuracy, we take into account
relativistic corrections, which were not previously considered in \cite{luchinsky,neng}.
The numerical contribution of different decay mechanisms
is analyzed and it is shown that in the case of double production of charmonium, the main decay mechanism
is the quark-photon one in Fig.\ref{pic2}, ZZ-mechanism in Fig.\ref{pic5}
although taking into account other mechanisms is also necessary to achieve high calculation accuracy
(see also \cite{bodwin2014}). For bottomonium pair production
all mechanisms give close contributions to the width, but the ZZ-mechanism is dominant.
Their sum ultimately determines the full numerical result.
The parameters A and B, which determine the numerical values of the contributions
of the W-boson loop mechanism and the quark loop mechanism, are obtained in analytical and numerical form.
They are presented separately in Table~\ref{tb22}.

In our approach, we distinguish two types of relativistic corrections. The corrections of the first type are
determined by the relative momenta p and q in the production amplitude of two quarks and two antiquarks.
The corrections of the second type are determined by the transformation law of the meson wave functions,
which results in expressions \eqref{eq3}-\eqref{eq4}. The $\psi({\bf p})$ wave functions entering into
\eqref{eq3}-\eqref{eq4} in the
rest frame of the bound state are found from the solution of the Schr\"odinger equation with a potential that includes relativistic corrections. Having an exact expression for the decay amplitude (see the example in the Appendix C),
we carry out a series of transformations with it, extracting second-order corrections in p and q. The model used
is described in more detail in \cite{apm2022}. In our approach, all arising relativistic parameters $\omega_1$,
$\omega_2$, $\tilde R(0)$ are determined within the framework of the relativistic quark model (see Table~\ref{tb1}).
Accounting for relativistic corrections in this work shows that such contributions lead to
a significant change in nonrelativistic results. Here it must be emphasized that we call nonrelativistic
results such results that are obtained at $p=0$, $q=0$ and neglecting relativistic corrections in the quark
interaction potential. The main parameter that greatly reduces the nonrelativistic results is $\tilde R(0)$,
which enters the fourth power of the decay width.
Therefore, the difference between the relativistic and nonrelativistic results in Table~\ref{tb2} turns out
to be more significant in the case of pair production of charmonium.

The paper considers various mechanisms for the production of a pair of vector mesons in the decay of the
Higgs boson, which are presented in Fig.~\ref{pic1}-Fig.~\ref{pic5}. There are other amplitudes
for the production of a pair of $J/\Psi$, $\Upsilon$ mesons,
which are calculated but not included in detail in the work. So, for example, there is a pair production mechanism
(HH mechanism), when the original Higgs boson turns into a HH pair, which then gives a pair of vector mesons
in the final state. In a sense, it is similar to the ZZ mechanism, but the amplitude structure is different,
which results in a contribution to the decay width that is several orders of magnitude smaller than those
given in Table~\ref{tb1}.
Despite the obvious difference in the amplitudes in Fig.~\ref{pic1} and Fig.~\ref{pic2} in terms of the mass factor
$M^2/M_h^2$,
we include the amplitudes of Fig.~\ref{pic1} in the consideration, in contrast to work \cite{neng}.
We do not consider amplitudes like those shown in Fig.~\ref{pic3}, but with two gluons. In this case, the direct
amplitude (Fig.~\ref{pic3} (left)) vanishes due to the color factor, and the cross amplitude
(Fig.~\ref{pic3} (right)) is suppressed by the mass factor $M^2/M_h^2$.
A distinctive feature of our calculations is that when studying the contributions of the amplitudes
in Fig.~\ref{pic3}-\ref{pic4}, we take into account two structure functions A and B \eqref{eq13d} in the tensor
function of the triangular loop. They are calculated analytically and numerically, and the corresponding
results are presented in Table~\ref{tb3} and Appendixes A and B.
Although the function B initially contains higher powers of the factor $r_4$ (see Appendix A, B),
which should lead to a decrease in the numerical values for $r_4 \ll 1$, nevertheless, the structure
of the coefficient functions in \eqref{eq18} is such that the function B cannot be neglected.

In Table~\ref{tb3}, we present separately the numerical values of the contributions from different decay mechanisms
with an accuracy of two significant figures after the decimal point. First of all, the problem was to obtain
the main contribution to the decay width. This was the contribution from the quark-photon mechanism in
Fig.~\ref{pic2} and the ZZ mechanism in Fig.~\ref{pic5}, which are one order of magnitude
greater than the other contributions in the case of charmonium production. Its value is due to the structure
of the decay amplitudes, in which small denominators $1/M^2$ appear from the photon propagators
in contrast to other amplitudes in which there is a factor $1/M_H^2$. In addition,
the numerator in these amplitudes contains amplifying factors in powers of a large parameter $r_1$. In each
considered decay mechanism, there are radiative corrections $O(\alpha_s)$ that were not taken into account. They represent
a major source of theoretical uncertainty, which we estimate to be about $30\%$ in $\alpha_s$. Other available
errors connected with the parameters of the Standard Model and higher order relativistic corrections
do not exceed $10\%$.
On the whole, our complete numerical estimates of the decay widths agree in order of magnitude
with the results \cite{neng}.

\acknowledgments
The work is supported by the Foundation for the Advancement
of Theoretical Physics and Mathematics "BASIS" (grant No. 22-1-1-23-1).
The authors are grateful to A.~V.~Berezhnoy, I.~N.~Belov, D.~Ebert, V.~O.~Galkin, A.~K.~Likhoded
for useful discussions.