# Application of deep learning to the estimation of normalization coefficients
in diffusion-based covariance models

###### Abstract

Variational data assimilation in ocean models depends on the ability to model
general correlation operators in the presence of coastlines. Grid-point
filters based on diffusion operators are widely used for this purpose, but
come with a computational bottleneck – the costly estimation of normalization
factors for every model grid point. In this paper, we show that a simple
convolutional neural network can effectively learn these normalization factors
with better accuracy than the current operational methods. Our network is
tested with a two-dimensional diffusion operator from the NEMOVAR ocean data
assimilation system, applied to a global ocean grid with approximately one
degree horizontal resolution. The network is trained on exact normalization
factors estimated by a brute-force method. Knowing that convolutional networks
can only model translation-equivariant functions, we ensure that the
normalization estimation problem is indeed translation-equivariant.
Specifically, we show how the number of inputs of this problem can be reduced
while preserving translation equivariance. Adding the distance to the
coastline as an input channel is found to improve the performance of the
network around coastlines. Extensions to three-dimensional diffusion and to
higher horizontal resolutions are discussed. Removing the computational
bottleneck associated with normalization opens the way to using adaptive
correlation models for operational ocean data assimilation. The code for this
work is publicly available at https://github.com/FolkeKS/DL-
normalization/tree/core-features

##  1 Introduction

Error covariance matrices are essential for determining how different pieces
of information are weighted in a data assimilation (DA) problem. For large
problems, such as those typically encountered in atmospheric and ocean DA
applications, it is impossible to represent them as explicit full-rank
matrices. Instead, they must be modeled, either as (localized) reduced-rank
matrices or in terms of parameterized covariance functions that have a
tractable number of tunable parameters and that are computationally efficient
when applied as the kernel of a covariance operator (Bannister, 2008).
Covariance operators are commonly used in variational DA where only matrix-
vector products are required for the iterative algorithms used to obtain the
solution. In this article, we consider covariance matrices belonging to this
latter (operator) class. Specifically, we focus on covariance operators that
can be formulated as a sequence of differential operators and interpreted as
the solution of a diffusion equation (Weaver and Courtier, 2001; Mirouze and
Weaver, 2010; Weaver and Mirouze, 2013). When the diffusion equation is solved
implicitly with an Euler-backward step, diffusion operators represent
covariance operators whose smoothing kernels are from the Matérn class
(Guttorp and Gneiting, 2006; Lindgren et al., 2011).

Diffusion-based covariance operators are convenient in DA applications
involving domains with complex boundary geometry and requiring inhomogeneous
and anisotropic smoothing tensors. This is the case for background-error
covariances in ocean DA (Weaver and Courtier, 2001; Weaver et al., 2020).
Specifically, and especially within the context of this article, these
operators are needed to model correlation operators, which, by definition,
correspond to covariance matrices whose diagonal elements (variances) are all
equal to one. The symmetric property of the covariance matrix is maintained by
pre- and post-multiplying the diffusion-based correlation matrix by the
diagonal matrix of standard deviations. By factoring the covariance operator
this way, the correlations and variances can be controlled independently.

In order to isolate the correlation operator, diffusion-based covariance
matrices must be normalized by their intrinsic variances. We refer to the
diagonal matrix whose elements are the inverse of the intrinsic variances as
the matrix of normalization factors. The normalization factors act to
compensate the attenuation of the covariance amplitude that occurs from a
diffusive process. As with the actual variances, the normalization factors can
be applied symmetrically through pre- and post-multiplication by a diagonal
matrix that contains the square root of the normalization factors. The
normalization factors have units of length, area or volume depending on
whether the covariance operator is applied in one, two or three dimensions.

Except in idealized or specific cases, for which we can resort to well-known
analytical expressions (Guttorp and Gneiting, 2006), the normalization factors
are not known to high accuracy at all grid points. Their values depend on the
diffusion model parameters, such as the smoothing tensor; on the boundary
geometry, such as coastlines in an ocean model; on the type of boundary
conditions (Neumann, Dirichlet or mixed); and on the details of the numerical
discretization scheme and computational solver. Regardless of the domain and
covariance model details, it is always possible to diagnose the exact
normalization factor at any grid point using brute force, whereby the
covariance matrix is applied to a canonical vector, which has a value of one
at the grid point of interest and a value of zero at all other grid points.
However, this method is not practical for model configurations involving
millions of grid points, especially if normalization factors have to be
regularly recomputed, as is required in a cycled DA system when newly
available estimates of the diffusion model parameters are applied on each
cycle. A randomization or Monte Carlo method is much more practical than the
brute-force method as it can be used to generate sufficiently accurate
normalization factors for any number of grid points at the cost of several
thousand applications of the diffusion operator. Nevertheless, this is still
too costly to be applied routinely (i.e., on each DA cycle) to three-
dimensional (3D) variables in a real-time operational environment where
timeliness is critical.

Using a diffusion-based covariance model for a global ocean DA system, Weaver
et al. (2020) compared randomization with various other methods for
approximating the normalization factors. Their work showed that none of the
other methods is accurate and robust enough to provide a viable alternative to
randomization. As a compromise between allowing for adaptive covariances on
the one hand, and keeping computational costs acceptable on the other, Weaver
et al. (2020) proposed a quasi-adaptive correlation model for 3D variables
that relies on splitting the diffusion operator into a one-dimensional (1D)
vertical component and a two-dimensional (2D) horizontal component (referred
to by Weaver et al. (2020) as a 2D$\times$1D formulation) and only allowing
the 1D component to be adaptive. Furthermore, they proposed a method whereby
the normalization factors can be well approximated by a product of
normalization factors computed separately for the 1D and 2D diffusion-model
components. This way, the normalization factors for the adaptive and
inexpensive 1D component can be estimated “on-the-fly” to a high accuracy
using the randomization or even brute-force method. The normalization factors
for the static and expensive 2D component applied in each model level need to
be estimated only once, using the conventional randomization approach. In this
article, we show how 2D normalization factors in the 2D$\times$1D formulation
can be estimated accurately and cheaply using a machine-learning strategy,
thus opening the way for the use of 2D adaptive correlations.

The article is organized as follows. In Section 2, we show how the problem of
estimating normalization factors for the diffusion operator can be transformed
to make it suitable for a solution algorithm based on a Convolutional Neural
Network (CNN). Specific care is taken to ensure that the function to regress
does indeed lie in the class of functions that can be approximated with any
accuracy by CNNs. The method has been applied to a 2D diffusion operator that
forms the basis of a correlation model in a global ocean DA system. Results
are presented and discussed in Section 3 where they are compared to those
computed with a standard randomization method. A summary and outlook are given
in Section 4.

##  2 Proposed method

We present here the machine learning approach used to estimate normalization
factors for a diffusion-based covariance operator. For at least two reasons,
the natural choice for this problem is a Convolutional Neural Network (CNN).
First, the data are stored on a structured grid, where every point has
neighboring grid points in all four directions. This makes it very similar to
an image-to-image translation exercise, where CNNs are known to perform well.
The structured nature of the data allows us to use standard CNNs, instead of
more general but less efficient Graph Neural Networks for instance. Second,
the normalization problem is essentially local. While a normalization factor
at a given location is, in theory, affected by the diffusion tensor at all
points on the globe, in practice, it is mainly sensitive to the diffusion
tensor at points in a small neighborhood. The impact of distant diffusion
tensors decays roughly with the amplitude of the modeled correlation function.
In DA applications, the modeled correlations are approximately Gaussian, with
length scales that are short compared to the domain size. As a result, they
decay to zero quite rapidly and the problem can effectively be considered as
local. This makes CNNs with even relatively small receptive fields good
candidates for solving the normalization problem.

###  2.1 Reducing the dimension of the problem

We now explain how the problem can be simplified by reducing the number of
inputs required to train the network. No accuracy is lost in the process. Our
assumption is that reducing the number of inputs (a.k.a. channels) will make
the problem easier to learn for a neural network, and will reduce the memory
cost and execution time of both training and applying the network.

We want to approximate a non-linear function

$\displaystyle
f:(\bm{\kappa}_{11},\bm{\kappa}_{22},\mathbf{e}_{1v},\mathbf{e}_{2v},\mathbf{e}_{1u},\mathbf{e}_{2u},\mathbf{w},\mathbf{m})\;\to\bm{\gamma}$
(1)

where $\bm{\kappa}_{11}$ and $\bm{\kappa}_{22}$ are the diagonal elements of
the horizontal diffusion tensor, with physical units of squared length;
$\mathbf{e}_{1v},\mathbf{e}_{2v},\mathbf{e}_{1u}$ and $\mathbf{e}_{2u}$ are
scale factors that define the local directional grid lengths at the $v$\- and
$u$-points indicated in Figure 1; $\mathbf{w}:=\mathbf{e}_{1T}\mathbf{e}_{2T}$
is the grid cell area defined at the cell center or $T$-point in Figure 1;
$\mathbf{m}$ is the land-sea mask, also defined at the $T$-point; and
$\bm{\gamma}$ denotes the associated normalization factor. Here, we are using
conventional notation from the ocean model NEMO (Nucleus for European
Modelling of the Ocean; Madec et al., 2023). All of these inputs are vectors
of the size of the grid, some of them being located in the centers of the
cells, some on the edges. Note that we limit ourselves to the case of diagonal
diffusion tensors here, which implies that the principal directions of
anisotropy are assumed to be aligned with the grid edges. This assumption is
currently used in the operational ocean DA systems at ECMWF and Met Office.
The diagonal elements of the tensor, $\bm{\kappa}_{11}$ and
$\bm{\kappa}_{22}$, are hereafter called diffusivity fields or diffusivities.
The problem also depends on external parameters that we consider as constant
here. These parameters include the total number of diffusion steps $M$ and the
type of boundary condition. In our experiments, $M$ is set to 10 to represent
an approximate Gaussian function, and a Neumann boundary condition is imposed
where the normal derivative at coastlines is set to zero.

Figure 1: The diffusivities $\bm{\kappa}_{11}$ and $\bm{\kappa}_{22}$ are,
respectively, located at a $u$-point and $v$-point on a staggered Arakawa
C-grid. The normalization factors $\bm{\gamma}$ are located at $T$-points.
Here, we scale $\bm{\kappa}_{11}$ and $\bm{\kappa}_{22}$ by the ratios
$\mathbf{e}_{2u}/\mathbf{e}_{1u}$ and $\mathbf{e}_{1v}/\mathbf{e}_{2v}$,
respectively.

In the specific context of diffusion, we can simplify the problem by
performing a local re-scaling of the diffusivity fields $\bm{\kappa}_{11}$ and
$\bm{\kappa}_{22}$. This allows us to account for the land-sea mask and for
the grid structure implicitly. More precisely, problem (1) can be transformed,
without approximation, as:

$\displaystyle\widetilde{f}:(\bm{\alpha}_{1},\bm{\alpha}_{2},\mathbf{w})\;\to\;\bm{\gamma}$
(2)

where we have introduced the scaled diffusivity fields

$\displaystyle\bm{\alpha}_{1}=\bm{\epsilon}_{u}\frac{\bm{\kappa}_{22}\bm{e}_{1v}}{\bm{e}_{2v}}\hskip
14.22636pt\mbox{and}\hskip
14.22636pt\bm{\alpha}_{2}=\bm{\epsilon}_{v}\frac{\bm{\kappa}_{11}\bm{e}_{2u}}{\bm{e}_{1u}},$

$\bm{\epsilon}_{u}$ and $\bm{\epsilon}_{v}$ being defined from the land-sea
mask $\mathbf{m}$ and from the boundary condition type, as explained in
Mirouze and Weaver (2010, their section 5.2). This re-formulation is made
possible by noting that the discretized diffusion operator can be expressed in
terms of $\bm{\alpha}_{1}$, $\bm{\alpha}_{2}$ and $\mathbf{w}$ only, as shown
by Mirouze and Weaver (2010). In our case, we used (zero) Neumann boundary
conditions, which is equivalent to setting $\bm{\epsilon}_{u}$ (respectively,
$\bm{\epsilon}_{v}$) equal to zero on east-west (respectively, north-south)
coastlines and equal to one elsewhere. Other types of boundary conditions
would lead to different values of the $\bm{\epsilon}_{u}$ and
$\bm{\epsilon}_{v}$ fields on the coastlines, as detailed in Mirouze and
Weaver (2010).

We remark that there is no loss of information through this transformation, as
we have

$\displaystyle
f(\bm{\kappa}_{11},\bm{\kappa}_{22},\mathbf{e}_{1v},\mathbf{e}_{2v},\mathbf{e}_{1u},\mathbf{e}_{2u},\mathbf{w},\mathbf{m})=\widetilde{f}(\bm{\alpha}_{1},\bm{\alpha}_{2},\mathbf{w}).$
(3)

The problem is thus reduced from eight input fields to three. Importantly,
both formulations of the normalization problem are equivariant to translation,
a key property that we will explore further in the following subsection.

###  2.2 Equivariance to translation

A fully convolutional neural network represents functions that are equivariant
to translation (Goodfellow et al., 2016, Chapter 9.4). A function $f$ is
equivariant to a spatial translation $T$ if it satisfies the following
identity:

$\displaystyle f(T(\mathbf{x}))=T(f(\mathbf{x})),$ (4)

where $\mathbf{x}$ denotes the inputs of the function, which should match the
inputs used by the corresponding network. As the CNN is itself equivariant to
translation by construction, it is intrinsically limited on how well it can
fit a function that does not satisfy this property.

Both formulations (1) and (2) of the normalization problem are equivariant to
translation. The equivariance to translation of formulation (2) is
demonstrated in Appendix A. For the sake of simplicity, the analysis is
restricted to one-dimensional (1D) diffusion, although the same line of
reasoning can be applied to 2D diffusion. The equivariance to translation
property is illustrated numerically in Figure 2. This figure has been obtained
using a 2D diffusion model in Python that closely reproduces the 2D diffusion
on the sphere from NEMOVAR but uses a regular latitude/longitude grid.

Figure 2: Numerical test illustrating the equivariance relation from Equation
(4). Panel (a) shows $f(\mathbf{x})$, where $f$ is defined as in Equation (1)
and $\mathbf{x}$ corresponds to random diffusivity fields $\bm{\kappa}_{11}$
and $\bm{\kappa}_{22}$, and the coefficients
$\mathbf{e}_{1v},\mathbf{e}_{2v},\mathbf{e}_{1u},\mathbf{e}_{2u},\mathbf{w}$
and $\mathbf{m}$ are associated with the regular grid used with this test
problem. Panel (b) shows $T(f(\mathbf{x}))$, panel (c) shows
$f(T(\mathbf{x}))$, and panel (d) shows the relative difference between (b)
and (c) (i.e., $|1-f(T(\mathbf{x}))/T(f(\mathbf{x}))|$ ). The relative
difference (d) is equal to zero at machine precision.

It is important to remark that the equivariance to translation property would
not be verified if we neglected to provide the grid information to the neural
network, assuming that it could learn them, as in

$\displaystyle
f_{(\mathbf{e}_{1v},\mathbf{e}_{2v},\mathbf{e}_{1u},\mathbf{e}_{2u},\mathbf{w})}:(\bm{\kappa}_{11},\bm{\kappa}_{22},\mathbf{m})\;$
$\displaystyle\to\bm{\gamma}$ (5)
$\displaystyle\mbox{or}\;\;f_{(\mathbf{e}_{1v},\mathbf{e}_{2v},\mathbf{e}_{1u},\mathbf{e}_{2u})}:(\bm{\kappa}_{11},\bm{\kappa}_{22},\mathbf{m},\mathbf{w})\;$
$\displaystyle\to\bm{\gamma}.$ (6)

In this case, a pure CNN could not fit the function with an arbitrarily small
error. In formulation (1), the grid parameters are needed to preserve
equivariance to translation but mostly represent redundant information, making
formulation (2) the more natural choice.

###  2.3 Generation of the data set

There are two steps in the generation of the training, validation, and test
data sets: 1) design realistic input diffusivity fields; and 2) compute the
associated normalization fields. We use the so-called ORCA1 configuration of
NEMO (Madec et al., 2023), which refers to the standard tri-polar global grid
with approximately one-degree horizontal resolution. All fields are defined in
the surface level of ORCA1 in our 2D diffusion experiments.

For the first step, we generate reference diffusivity fields and add random,
spatially-correlated perturbations to them to produce a set of distinct input
samples, with complex spatial variations similar to those that would be
obtained from ensemble methods such as the one used by Weaver et al. (2020).
The reference field is generated following the method described in Mogensen et
al. (2012, their Section 4.6.2 and their Table 2). The reference field is then
perturbed by applying a multiplicative factor defined as a spatially
correlated Gaussian random field centered on 1.0 and with standard deviation
0.15. The correlation model used to define the perturbations is the same
correlation model associated with the reference diffusivity fields. Some
additional procedures are then applied to the perturbed diffusivity fields.
First, we reduce the value of a diffusivity field near coastlines when the
geometric mean of the implied directional length-scales at the diffusivity
grid point exceeds the distance to the coastline from that grid point. Second,
we increase the value of a diffusivity when the geometric mean of the implied
directional length-scales is less than the geometric mean of the grid scale
factors at the diffusivity point. We generate a total of 190 samples, which
differ from each other by the multiplicative noise factor only.

For the second step, we calculate the normalization factors for each of the
190 input samples using the brute-force method described in Section 1. This
method is numerically affordable for a 2D field in our relatively coarse
global ocean configuration. However, for a 3D field or with higher resolution,
the method would be extremely expensive and not practical for generating large
data sets. This issue is discussed in the last section of this article.

An example of an input sample of diffusivity fields along with their
corresponding true value of the normalization field is illustrated in Figure
3.

Figure 3: An example of a sample of randomly-perturbed diffusivity fields
$\bm{\kappa}_{11}$ and $\bm{\kappa}_{22}$ (left and middle panels) and the
associated normalization factors (right panel).

###  2.4 Pre-processing of the data

As explained in Section 22.1, the diffusivity fields are scaled to reduce the
number of inputs while preserving equivariance to translation. We also
observed that adding a distance-to-coast input channel reduces large errors
near coastlines. The results we present here are all obtained with this
additional input channel. The data are then separated into training,
validation, and test sets containing 10, 10 and 170 samples, respectively. We
deliberately limited ourselves to 10 input/output training pairs to reflect
the increased cost of data generation when considering the problem in higher
dimensions or with higher resolution. Using statistics from the training data
set, each of the four input channels and the output channel are centered with
respect to their global mean and normalized by their standard deviation.

Finally, the inputs are extended by periodicity in both the longitudinal and
latitudinal directions by a few grid points. This ensures that no assumptions
need to be made by the neural network near the boundaries of the “image” while
keeping the network output at the correct (original) dimensions. This
extension by periodicity accounts for the fact that the input channels are
defined on different points on the staggered grid, as illustrated in Figure 1.

###  2.5 Model architecture and training specifications

The final network, inspired by the ResNet of He et al. (2015), consists of
Residual Blocks (ResBlocks) chained together as illustrated in Figure 4. A
ResBlock uses skipped connections to add, or otherwise transfer, the output of
one layer to a layer deeper in the network. Our version of a ResBlock is
detailed in Figure 5. The network contains 4 ResBlocks along with a layer for
both the network input and output, totaling 10 convolutional layers and
approximately 3 $\times 10^{5}$ trainable parameters.

Figure 4: Architecture containing four ResBlocks along with convolutional
layers for the network input and output, marked in orange. The first layer
consists of 64 filters, each convolving over the 4 input channels and making
the input of a ResBlock contain 64 feature maps. Inversely, the final layer
consists of only one filter to obtain the correct output dimension of
$292\times 360\times 1$.
Figure 5: Our implementation of a ResBlock consists of two 3x3x64
convolutional layers with 64 filters, each followed by batch normalization and
the application of an Exponential Linear Unit (ELU) activation function. The
inputs of a block are added to the output of the second convolutional layer
before applying the activation function.

We use a relative loss function where the continent grid points are masked;
i.e., the model is only penalized for errors on grid points in the ocean:

$\displaystyle\ell=\frac{1}{N}\sum^{N}_{n=0}\left(\frac{\hat{\gamma}_{n}^{2}-\gamma_{n}^{2}}{\gamma_{n}^{2}}\right)^{2}.$

where $\gamma_{n}^{2}$ is the true normalization factor,
$\hat{\gamma}_{n}^{2}$ is the predicted normalization factor, and $N$ is the
number of ocean grid points.

The model is trained from scratch using Pytorch version 1.11.0 (Paszke et al.,
2019) default weight initialization and the Adam optimizer with default
learning rate of $10^{-3}$.

The hardware used is an Nvidia Tesla v100 16GB GPU. The best model is selected
based on performance on the validation data set after 3000 epochs, which
corresponds to approximately 9 hours of training. This high number of epochs
is chosen because we use a small training data set (less than 20 MB of data),
and because we observed a continued increase in accuracy both on the training
and validation data sets as the number of epochs increases.

##  3 Evaluation of the neural network

The neural network described in the previous section has been evaluated
against the randomization method with a large sample size of $10^{4}$ samples.
This sample size is a typical choice for this application, giving estimates
that are considered accurate enough for operational use (Weaver et al., 2020).

###  3.1 Error metrics

We define two global error metrics to evaluate and compare different
estimation methods. Following Weaver et al. (2020), we define the relative
error $\varepsilon_{n}$ for an ocean grid point $n$ as

$\displaystyle\varepsilon_{n}=\frac{\hat{\gamma}_{n}^{2}-\gamma_{n}^{2}}{\gamma_{n}^{2}},$

where $\gamma_{n}^{2}$ is the true factor computed with the brute-force method
and $\hat{\gamma}_{n}^{2}$ is the estimate obtained using either the
randomization method or the CNN.

From this local quantity, we define two global precision measures. The first
is the mean absolute relative error $\varepsilon$, which is defined as

$\displaystyle\varepsilon=\frac{1}{d}\sum^{d}_{l=0}\frac{1}{N}\sum^{N}_{n=0}|\varepsilon_{n}|,$

where $N$ is the number of ocean grid points and $d$ is the number of samples
available. For the neural network evaluation, we present results on the test
data set with $d=170$. For the randomization method, we only have $d=1$. This
is enough to have an accurate estimate, as we know the randomization errors
are “well-behaved” in the sense that the variance estimates are unbiased and
follow a $\chi^{2}$ distribution.

While we want the estimate to be accurate on average, we also want to control
large local errors. We therefore introduce the mean maximum absolute relative
error $\varepsilon^{\infty}$:

$\displaystyle\varepsilon^{\infty}=\frac{1}{d}\sum^{d}_{l=0}\max\limits_{0\leq
n\leq N}|\varepsilon_{n}|.$

Although the neural network is optimized to achieve a small mean square
relative error through the loss function $\ell$ defined in the previous
section, the final network is chosen based on the values of both $\varepsilon$
and $\varepsilon^{\infty}$ on the validation data.

###  3.2 Results

All neural network results shown in this section are from the test data set.
The performance metrics of the neural network and of a randomization method
with $10^{4}$ samples are compared in Table 1. The neural network consistently
beats the randomization method in terms of mean absolute relative error
(approximately 3 times smaller). It has a similar performance in terms of mean
maximum error.

Method | $\varepsilon^{\infty}$ | $\varepsilon$
---|---|---
Randomization with $10^{4}$ samples | 5.91% | 1.15%
CNN | 4.3% | 0.4%
Table 1: Relative errors of the two estimation methods.

To have more information on the distribution of errors, we show the histogram
of the relative errors across all grid points and all test samples in Figure
6. Similarly, Figure 7 shows the maximum errors of the different test samples.

Figure 6: Histogram of relative errors across all grid points and all test
samples. Each bar has a width of approximately $0.53\%$. The upper and lower
bounds are set to $\pm 5\%$, leaving some extreme values hidden. These values
are highlighted in Figure 7.  Figure 7: Histogram of the maximum absolute
relative error. Each bar has a width of approximately $1.73\%$. The maximum
values are taken from each test sample.

Figure 6 shows that the neural network is globally unbiased, with errors
following a Gaussian-like distribution centered on zero. A closer look at the
distribution of the maximum absolute error in Figure 7 reveals an interesting
outlier for one test sample, with a maximum error of $0.32$. After
investigation, this error was found to be caused by a single grid point with
coastlines on several edges. In this grid cell, the corresponding input
diffusivity field $\bm{\kappa}_{11}$ has an abnormally large perturbed value,
even after its original value is reduced by the distance-to-the-coast control
procedure outlined in Section 22.3. The anomalous scaled diffusivity has a
value of 10.2 after the standardization procedure explained in Section 22.4.
For comparison, the second largest scaled diffusivity in the test data set is
8.6, and the largest scaled diffusivity in the training data set is 7.5 (using
the same standardization). Importantly, this result highlights the fact that
the network can make large errors for inputs that are too far from the range
of values in the training data.

An example of the neural network outputs on a random sample of the test data
set is shown in Figures 8 and 9.

Figure 8: A sample from the test data set showing the true normalization
factors and the factors estimated by the network, together with their relative
error. Figure 9: Same as Figure 8 but with a zoom over coastal European seas.
Figure 10: Mean absolute relative error as a function of the Manhattan
distance to the coastline. These results were obtained from a network with
distance-to-coast as an additional input channel.

There is no visible spatial structure in the errors, except when close to the
coastline. This is evident in Figure 10, which shows the error as a function
of the distance to the coastline. It clearly shows that the error is larger
for points within 10 grid points of the coastlines than for points further
away. For points adjacent to a coastline, the error is up to 10 times larger.
This observation motivated the introduction of the distance-to-coast as a
fourth input of the neural network, which was in fact used to produce the
results in Figure 10. Although the errors are still relatively larger in the
vicinity of the coastlines, this extra input channel helps maintain the
amplitude of the error below that of the randomization method.

##  4 Summary and outlook

In this paper, we evaluated the potential of a convolutional neural network
(CNN) to estimate normalization factors for a covariance model based on a
diffusion operator. The normalization factors are required to transform the
implied diffusion-based covariance matrix into a covariance matrix with unit
diagonal (i.e., a correlation matrix). We trained a CNN on pairs of 2D
diffusivity fields and the associated exact 2D normalization factors. The data
were generated on a global grid of the NEMO ocean model (Madec et al., 2023)
using a 2D diffusion operator from the NEMOVAR ocean data assimilation system
(Mogensen et al., 2012), with randomly perturbed diffusivity fields. The
trained network yields more accurate estimates than those obtained by a
randomization method with 10000 samples, which is the current method of choice
for this problem. Although it has not been quantified in this article, once
trained, the CNN runs at a fraction of the computational cost of
randomization.

The use of machine learning for determining normalization factors is being
explored independently at the UK Met Office (James While, personal
communication) and at the European Centre for Medium-Range Weather Forecasts
(ECMWF; Marcin Chrust and Vincent Chabot, personal communication), using
variations of the approach described in this article. We would like to
highlight three aspects of our work that we believe to be important and
relevant for other machine learning approaches to this problem.

  1. 1.

When using a CNN for this problem, the input diffusivity fields should be
scaled (see Section 22.2) to reduce the number of input channels, while
ensuring the problem lies in the class that can be represented by the network.
The appropriate scaling involves the local grid geometry and the land-sea mask
used to account for the boundary conditions at coastlines. Equivariance to
translation can then be obtained by using the scaled diffusivity fields
together with the local grid-cell area elements as input channels.

  2. 2.

The normalization factors are harder to determine in the vicinity of
coastlines, particularly where the geometry is complex. In practice, adding
the distance to the nearest coastline as an extra input channel to the neural
network helps to reduce the estimation error. The distance-to-coast input
channel may become more important if the presence of a coastline outside the
receptive field of the network impacts the normalization factor at a given
grid point, which might occur, for instance, when moving to higher resolution
grids.

  3. 3.

A surprisingly small training data set (10 pairs of 2D fields in our case) is
enough to outperform the randomization method with 10000 samples.

The network presented here is trained on a global ocean grid with coarser
resolution than used operationally at ECMWF and the Met Office. Here, the
input and output fields are defined at 1 degree horizontal resolution compared
to the common operational resolutions of 1/4 or 1/12 degree. However, we
expect that a CNN with similar structure would perform equally well at higher
resolutions, although the required receptive field of the neural network would
possibly change due to larger diffusivities relative to the increased grid
resolution.

In principle, the network could be extended to account for 3D diffusion or
2D$\times$1D diffusion (see Section 1) by using 3D convolutional kernels and
by taking care to apply the same scaling procedure to the vertical
diffusivities. Within the specific 2D$\times$1D framework, which is the 3D
diffusion formulation used operationally at ECMWF and the Met Office, the
normalization approximation based on the separable assumption described in
Section 1 could be used with the existing 2D convolutional approach to reduce
computational costs. While the geometry of the coastlines would change in each
ocean model level, the analytical problem itself would not change, and the
approach we have presented should be able to account for any coastline
structure. Including data from different ocean model levels in the training
step would likely improve the accuracy, but would not be essential.

Even for the 2D diffusion problem, generating a training data set of exact
normalization factors would be prohibitively expensive at 1/4 and 1/12 degree
resolutions. Indeed, the cost of a single application of the diffusion
operator scales at least linearly with the number of grid points, and the
brute-force approach requires as many applications as grid points. As a
result, the brute-force approach scales at least quadratically with the number
of grid points. Conversely, the computational cost of the randomization method
only increases as the cost of a single application of the diffusion operator
when increasing resolution. The intersection point where the two methods have
comparable costs is found when there are as many grid points as random samples
(typically 10000); i.e., for grids even smaller than those considered in this
article. Therefore, for high-resolution grids, the randomization approach
seems to be the only viable method to generate a training data set. From a
training point of view, it is unclear whether it is possible to beat the
randomization method by training on data generated by it. It might be possible
if the network does not try to overfit the data but manages to generalize
beyond the random error in the training data. It might also be desirable to
train on a larger random sample than is typical, to reduce sampling error as
much as possible.

An alternative or complementary solution would be to train the network on a
small selection of grid points only, for which the associated normalization
factors could be estimated with the brute-force method. These points could be
sampled homogeneously in the model domain, or could be sampled more heavily
where the network performs more poorly; i.e., near the coastlines. This
approach could be useful at least for evaluating and testing the network with
an accurate data set. It could also be used to provide refined training of a
neural network that is pre-trained with randomization-generated data.

Estimating normalization factors for a diffusion-based covariance operator is
a computationally expensive operation, but well suited for a machine learning
approach. Being able to estimate these factors accurately and quickly opens
opportunities for using adaptive covariance parameters diffusivity fields) in
operational ocean data assimilation. In particular, for centers that can
estimate background-error correlations from an ensemble of forecasts, this
would allow the use of correlations-of-the-day in parametric background-error
covariance models based on the diffusion operator. Furthermore, this would
enable the application of localization-of-the-day with ensemble-variational
data assimilation (EnVar) schemes that use a diffusion operator for localizing
a sample estimate of the background-error covariance matrix.

###### Acknowledgements.

The authors wish to thank Andrea Piacentini for his help in setting up and
running the NEMOVAR experiments. The code for this work is publicly available
at https://github.com/FolkeKS/DL-normalization/tree/core-features and the
associated data set at
https://kaggle.com/datasets/e6e1e745161300b59a7ecf7695f02fedd5ab6f5853067004c127029d15fde68e

##  Appendix A Equivariance to translation of one-dimensional diffusion

As introduced in Section 22.2, as fully convolutional networks are equivariant
to translation, they are bound to be inaccurate if used to model a function
that is not itself equivariant to translation. Figure 2 illustrates
experimentally that, with the right choice of inputs, the function
$\widetilde{f}$ from Equation (2), is equivariant yo translation. This
property can also be demonstrated analytically, which is the purpose of this
appendix.

To lighten notation, we will only consider here the case of 1D diffusion.
Consequently, there is only one input feature ($\bm{\alpha}$) to consider
instead of two ($\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}$):

$\widetilde{f}:(\bm{\alpha},\mathbf{w})\;\to\;\bm{\gamma}.$ (7)

Our goal is to show that

$\widetilde{f}(\mathbf{S}\bm{\alpha},\mathbf{S}\mathbf{w})\,=\,\mathbf{S}\widetilde{f}(\bm{\alpha},\mathbf{w})\,=\,\mathbf{S}\bm{\gamma},$
(8)

where $\mathbf{S}$ is the matrix associated with 1D translation that
periodically shifts elements of a vector by one position ‘upwards’. For
$p,q\in\llbracket 0,n-1\rrbracket$, the coefficient $[\mathbf{S}]_{p,q}$ on
the $p$-th line and $q$-th column of $\mathbf{S}$ is defined as

$[\mathbf{S}]_{p,q}=\begin{cases}1&\text{ if }p-q=1;\\\ 1&\text{ if
}p=n-1\text{ and }q=0;\\\ 0&\text{ otherwise.}\\\ \end{cases}$ (9)

The inverse translation matrix $\mathbf{S}^{-1}=\mathbf{S}^{\rm T}$
periodically shifts elements of a vector by one position ‘downwards’.

A diffusion operator $\mathbf{D}$ can be decomposed into a symmetric sequence
of sub-operators (see Equation (4)) from Weaver et al., 2020):

$\mathbf{D}=\left[\mathbf{A}^{-1}\right]^{M/2}\mathbf{W}^{-1}\left[(\mathbf{A}^{-1})^{\rm
T}\right]^{M/2},$ (10)

where $\mathbf{A}$ is the self-adjoint, positive-definite matrix that must be
inverted to perform one step of implicit diffusion, and $\mathbf{W}$ is a
diagonal matrix whose diagonal contains the elements of the vector
$\mathbf{w}$ in Equation (1). These elements correspond to the surface of each
grid cell in 2D and the length of each grid cell in 1D. The matrix
$\mathbf{A}$ is self-adjoint with respect to the inner product with weighting
matrix $\mathbf{W}$ (Weaver et al., 2020, p4), which implies that the matrix
$\mathbf{W}\mathbf{A}$ is symmetric; i.e., self-adjoint with respect to the
canonical inner product. In terms of $\mathbf{W}\mathbf{A}$, Equation (10) can
then be rewritten as

$\mathbf{D}=\left[\left(\mathbf{W}\mathbf{A}\right)^{-1}\mathbf{W}\right]^{M/2}\mathbf{W}^{-1}\left[\mathbf{W}\left(\mathbf{W}\mathbf{A}\right)^{-1}\right]^{M/2}.$
(11)

We are interested in determining the matrix expression for the ‘translated’
diffusion operator, denoted $s(\mathbf{D})$, obtained by translating
$\bm{\alpha}$ and $\mathbf{w}$. To do so, we first look for the translated
counterparts of $\mathbf{W}$ and $\mathbf{WA}$, which we denote
$s(\mathbf{W})$ and $s(\mathbf{WA})$, respectively. For any vector
$\mathbf{v}$ of size $n$, we denote by $[\mathbf{v}]_{p}$ its $p$-th element
with $p\in\llbracket 0,n-1\rrbracket$. We also assume periodic boundary
conditions, which we impose by introducing cyclical indexes such that
$[\mathbf{v}]_{-1}=[\mathbf{v}]_{n-1}$ and
$[\mathbf{v}]_{n}=[\mathbf{v}]_{0}$.

For any vector $\mathbf{x}$, we have

$[\mathbf{Wx}]_{p}=[\mathbf{w}]_{p}[\mathbf{x}]_{p}.$ (12)

Applying the translation $\mathbf{S}$ to $\mathbf{w}$ shifts periodically its
elements:

$[s(\mathbf{W})\,\mathbf{x}]_{p}=[\mathbf{w}]_{p+1}[\mathbf{x}]_{p}.$ (13)

As this relation holds for all $p$, we can apply the inverse periodical shift
to all indexes:

$[s(\mathbf{W})\,\mathbf{x}]_{p-1}\,=\,[\mathbf{w}]_{p}[\mathbf{x}]_{p-1}.$
(14)

In matrix form, this system of equations becomes

$\mathbf{S}^{-1}s(\mathbf{W})\,\mathbf{x}=\mathbf{W}\,\mathbf{S}^{-1}\,\mathbf{x}.$
(15)

As this equation must hold for any vector $\mathbf{x}$, we have

$s(\mathbf{W})\,=\,\mathbf{S}\,\mathbf{W}\,\mathbf{S}^{-1}.$ (16)

The same property can be derived for $s(\mathbf{WA})$. A direct relation
between the elements of the matrix $\mathbf{WA}$, the scaled diffusivity
vector $\bm{\alpha}$ and the vector of weighting coefficients $\mathbf{w}$ can
be deduced from Equations (76) and (77) of Mirouze and Weaver (2010). In a 1D
domain, and with the notation used in this article, we have

$\displaystyle[(\mathbf{WA})\mathbf{x}]_{p}\,=\,\left([\mathbf{w}]_{p}-[\bm{\alpha}]_{p-1}-[\bm{\alpha}]_{p}\right)[\mathbf{x}]_{p}-[\bm{\alpha}]_{p-1}[\mathbf{x}]_{p-1}-[\bm{\alpha}]_{p}[\mathbf{x}]_{p+1},$
(17)

for $p\in\llbracket 0,n-1\rrbracket$. The scaled diffusivity
$[\bm{\alpha}]_{p}$ is equal to $\kappa_{p}/e_{p}$ where $\kappa_{p}$ is the
diffusivity on the edge between the $p$-th and $(p+1)$-th cells, and $e_{p}$
is the distance between the centres of the $p$-th and $(p+1)$-th cells. With
$\bm{\alpha}$ and $\mathbf{w}$ translated by $\mathbf{S}$, Equation (17)
becomes

$\displaystyle[s(\mathbf{WA})\mathbf{x}]_{p}\,=\,\left([\mathbf{w}]_{p+1}-[\bm{\alpha}]_{p}-[\bm{\alpha}]_{p+1}\right)[\mathbf{x}]_{p}-[\bm{\alpha}]_{p}[\mathbf{x}]_{p-1}-[\bm{\alpha}]_{p+1}[\mathbf{x}]_{p+1}.$
(18)

As Equation (18) holds for all values of $p$, we can replace $p$ by $p-1$:

$\displaystyle[s(\mathbf{WA})\mathbf{x}]_{p-1}\,=\,\left([\mathbf{w}]_{p}-[\bm{\alpha}]_{p-1}-[\bm{\alpha}]_{p}\right)[\mathbf{x}]_{p-1}-[\bm{\alpha}]_{p-1}[\mathbf{x}]_{p-2}-[\bm{\alpha}]_{p}[\mathbf{x}]_{p},$
(19)

which is equivalent to

$\displaystyle[\mathbf{S}^{-1}s(\mathbf{WA})\mathbf{x}]_{p}\,=\,\left([\mathbf{w}]_{p}-[\bm{\alpha}]_{p-1}-[\bm{\alpha}]_{p}\right)[\mathbf{S}^{-1}\mathbf{x}]_{p}-[\bm{\alpha}]_{p-1}[\mathbf{S}^{-1}\mathbf{x}]_{p-1}-[\bm{\alpha}]_{p}[\mathbf{S}^{-1}\mathbf{x}]_{p+1}.$
(20)

In matrix form, this system can be written as

$\mathbf{S}^{-1}s(\mathbf{WA})\,\mathbf{x}\,=\,\mathbf{WA}\,\mathbf{S}^{-1}\,\mathbf{x},$
(21)

which leads to

$s(\mathbf{WA})=\mathbf{S}\mathbf{WA}\mathbf{S}^{-1}.$ (22)

By inserting Equations (16) and (22) into Equation (10), we obtain

$s(\mathbf{D})\,=\,\mathbf{S}\,\mathbf{D}\,\mathbf{S}^{-1}.$ (23)

The elements on the diagonal of $\mathbf{D}$ are the square of the inverse of
the normalization factors $\gamma$:

$([\bm{\gamma}]_{p})^{-2}\,=\,\mathbf{e}_{p}^{\rm
T}\,\mathbf{D}\,\mathbf{e}_{p},$ (24)

where $\mathbf{e}_{p}$ is a vector whose elements are equal to zero except at
position $p$ where the element is equal to 1. We can estimate the diagonal
elements of $s(\mathbf{D})$ from the inner product

$\displaystyle\mathbf{e}_{p}^{\rm T}\,s(\mathbf{D})\,\mathbf{e}_{p}$
$\displaystyle\,=\,\mathbf{e}_{p}^{\rm
T}\,\mathbf{S}\,\mathbf{D}\,\mathbf{S}^{-1}\mathbf{e}_{p}$
$\displaystyle\,=\,(\mathbf{S}^{-1}\mathbf{e}_{p})^{\rm
T}\,\mathbf{D}\,(\mathbf{S}^{-1}\mathbf{e}_{p})$
$\displaystyle\,=\,(\mathbf{e}_{p+1})^{\rm T}\,\mathbf{D}\,\mathbf{e}_{p+1}$
$\displaystyle\,=\,([\bm{\gamma}]_{p+1})^{-2}$
$\displaystyle\,=\,([\mathbf{S}\bm{\gamma}]_{p})^{-2},$ (25)

where we have used the property that $\mathbf{S}$ is orthogonal
($\mathbf{S}^{-1}=\mathbf{S}^{\rm T}$).

Equation (25) shows that applying a translation $\mathbf{S}$ to the input
features $\bm{\alpha}$ and $\mathbf{w}$ is equivalent to applying the same
translation to $\bm{\gamma}$, and thus that $\widetilde{f}$ is equivariant to
translation:

$\widetilde{f}(\mathbf{S}\bm{\alpha},\mathbf{S}\mathbf{w})=\mathbf{S}\widetilde{f}(\bm{\alpha},\mathbf{w})=\mathbf{S}\bm{\gamma}.$
(26)

The same result can be obtained using the ‘downwards’ translation,
$\mathbf{S}^{-1}=\mathbf{S}^{\rm T}$, instead of the ‘upwards’ translation,
$\mathbf{S}$. By applying either of these translations multiple times, this
result can be extended to any periodic translation. Note that this property
does not automatically extend to other choices of input features. For example,
this result would not hold if the cell surfaces $\mathbf{w}$ were not included
in the input features, or if, instead of $\bm{\alpha}$, the diffusivity field
$\bm{\kappa}$ was included without the scaling factors representing the local
grid geometry ($\mathbf{e}_{1v}$, $\mathbf{e}_{2v}$, $\mathbf{e}_{1u}$ and
$\mathbf{e}_{2u}$ in Section 22.1 ).

## References

  * Bannister (2008) Bannister, R. N., 2008: A review of forecast error covariance statistics in atmospheric variational data assimilation. II: Modelling the forecast error covariance statistics. Quarterly Journal of the Royal Meteorological Society, 134 (637), 1971–1996, 10.1002/qj.340.
  * Goodfellow et al. (2016) Goodfellow, I. J., Y. Bengio, and A. Courville, 2016: Deep Learning. MIT Press, Cambridge, MA, USA, http://www.deeplearningbook.org.
  * Guttorp and Gneiting (2006) Guttorp, P., and T. Gneiting, 2006: Studies in the history of probability and statistics XLIX: On the Matérn correlation family. Biometrika, 93, 989–995, 10.1093/biomet/93.4.989.
  * He et al. (2015) He, K., X. Zhang, S. Ren, and J. Sun, 2015: Deep residual learning for image recognition. arXiv, URL https://arxiv.org/abs/1512.03385, 10.48550/ARXIV.1512.03385.
  * Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Lindström, 2011: An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73 (4), 423–498, 10.1111/j.1467-9868.2011.00777.x.
  * Madec et al. (2023) Madec, G., and Coauthors, 2023: NEMO Ocean Engine Reference Manual. Zenodo, URL https://doi.org/10.5281/zenodo.8167700, 10.5281/zenodo.8167700.
  * Mirouze and Weaver (2010) Mirouze, I., and A. T. Weaver, 2010: Representation of correlation functions in variational assimilation using an implicit diffusion operator. Quarterly Journal of the Royal Meteorological Society, 136 (651), 1421–1443, 10.1002/qj.643.
  * Mogensen et al. (2012) Mogensen, K., M. A. Balmaseda, and A. T. Weaver, 2012: The NEMOVAR ocean data assimilation system as implemented in the ECMWF ocean analysis for System 4\.  (668), 59, 10.21957/x5y9yrtm, URL https://www.ecmwf.int/node/11174.
  * Paszke et al. (2019) Paszke, A., and Coauthors, 2019: Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 8024–8035, URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Weaver and Courtier (2001) Weaver, A., and P. Courtier, 2001: Correlation modelling on the sphere using a generalized diffusion equation. Quarterly Journal of the Royal Meteorological Society, 127 (575), 1815–1846, 10.1002/qj.49712757518.
  * Weaver et al. (2020) Weaver, A. T., M. Chrust, B. Ménétrier, and A. Piacentini, 2020: An evaluation of methods for normalizing diffusion-based covariance operators in variational data assimilation. Quarterly Journal of the Royal Meteorological Society, 10.1002/qj.3918.
  * Weaver and Mirouze (2013) Weaver, A. T., and I. Mirouze, 2013: On the diffusion equation and its application to isotropic and anisotropic correlation modelling in variational assimilation. Quarterly Journal of the Royal Meteorological Society, 139 (670), 242–260, 10.1002/qj.1955.


# Magnetic dipole moments of the spin-$\frac{3}{2}$ doubly heavy baryons

Ulaş Özdem<EMAIL_ADDRESS>Health Services Vocational School of Higher
Education, Istanbul Aydin University, Sefakoy-Kucukcekmece, 34295 Istanbul,
Turkey

###### Abstract

The magnetic dipole moments of the spin-$\frac{3}{2}$ doubly charmed, bottom
and charmed-bottom baryons are obtained by means of the light-cone QCD sum
rule. The magnetic dipole moments of these baryons encode essential knowledge
of their inner structure and shape deformations. The numerical results are
given as, $\mu_{\Xi_{cc}^{*++}}=2.94\pm 0.95$, $\mu_{\Xi_{cc}^{*+}}=-0.67\pm
0.11$, $\mu_{\Omega_{cc}^{*+}}=-0.52\pm 0.07$, $\mu_{\Xi_{bb}^{*0}}=2.30\pm
0.55$, $\mu_{\Xi_{bb}^{*-}}=-1.39\pm 0.32$, $\mu_{\Omega_{bb}^{*-}}=-1.56\pm
0.33$, $\mu_{\Xi_{bc}^{*+}}=2.63\pm 0.82$, $\mu_{\Xi_{bc}^{*0}}=-0.96\pm 0.32$
and $\mu_{\Omega_{bc}^{*+}}=-1.11\pm 0.33$, respectively.

Electromagnetic form factors, Magnetic dipole moment, Doubly heavy baryons,
Light-cone QCD sum rules

##  I Motivation

The doubly heavy baryons (DHBs) presumably contain two heavy quark and one
light quark. One of them was firstly announced by the SELEX Collaboration in
the decay mode $\Xi_{cc}^{+}\rightarrow\Lambda_{c}^{+}K^{-}\pi^{+}$ with the
mass $M_{\Xi^{+}_{cc}}=3519\pm 1$ MeV Mattson:2002vu . However, neither Belle
Chistov:2006zj , nor FOCUS Ratti:2003ez , nor BABAR Aubert:2006qw  could
confirm the DHBs in $e^{-}\,e^{+}$ annihilations. It is worth pointing out
that the analysis of the SELEX experiment with other experimental groups is
achieved through different production mechanisms. Therefore, the results of
the SELEX Collaboration cannot be ruled out. In 2017, LHCb Collaboration
observed another doubly heavy baryon $\Xi_{cc}^{++}$ in the mass spectrum of
$\Lambda_{c}^{+}\,K^{-}\pi^{+}\,\pi^{+}$ with the mass
$M_{\Xi_{cc}^{++}}=3621.40\pm 0.72\pm 0.27\pm 0.14$ MeV Aaij:2017ueg .
Recently, the LHCb Collaboration reconstructed their analysis via decay modes
$\Xi_{cc}^{++}\rightarrow\Lambda_{c}^{+}\,K^{-}\pi^{+}\,\pi^{+}$ and
$\Xi_{cc}^{++}\rightarrow\Xi_{c}^{+}\,\pi^{+}$ and they reported their mass
value as $M_{\Xi_{cc}^{++}}=3621.55\pm 0.23\pm 0.30$ MeV Aaij:2019uaz . The
investigation for the DHBs may provide with valuable knowledge for
comprehension of the nonperturbative QCD effects. One of the several point of
views which makes the physics of DHBs charming is that the binding of two
charm quarks and a light quark provides a unique perspective for dynamics of
confinement. The research of the properties of DHBs is one of the active and
interesting branches of particle physics. Therefore, the weak Albertus:2006ya
; Li:2017ndo ; Wang:2017mqp ; Wang:2017azm ; Shi:2017dto ; Shi:2019hbf ;
Shi:2019fph , strong Hu:2005gf ; Xiao:2017udy  and radiative decays Li:2017pxa
; Yu:2017zst ; Lu:2017meb ; Cui:2017udv , the magnetic dipole moments (MDMs)
Can:2013zpa ; Branz:2010pq ; Bose:1980vy ; Patel:2008xs ; SilvestreBrac:1996bg
; Patel:2007gx ; Gadaria:2016omw ; JuliaDiaz:2004vh ; Faessler:2006ft ;
Can:2013tna ; Li:2017cfz ; Bernotas:2012nz ; Lichtenberg:1976fi ; Oh:1991ws ;
Simonis:2018rld ; Liu:2018euh ; Blin:2018pmj ; Meng:2017dni ; Dhir:2009ax ;
Ozdem:2018uue  and masses Bagan:1992za ; Roncaglia:1995az ; Ebert:1996ec ;
Tong:1999qs ; Itoh:2000um ; Gershtein:2000nx ; Kiselev:2001fw ; Kiselev:2002iy
; Narodetskii:2001bq ; Lewis:2001iz ; Ebert:2002ig ; Mathur:2002ce ;
Flynn:2003vz ; Vijande:2004at ; Chiu:2005zc ; Migura:2006ep ; Albertus:2006ya
; Martynenko:2007je ; Tang:2011fv ; Liu:2007fg ; Roberts:2007ni ;
Valcarce:2008dr ; Liu:2009jc ; Alexandrou:2012xk ; Aliev:2012ru ; Aliev:2012iv
; Namekawa:2013vu ; Karliner:2014gca ; Sun:2014aya ; Chen:2015kpa ;
Sun:2016wzh ; Shah:2016vmd ; Kiselev:2017eic ; Chen:2017sbg ; Hu:2005gf ;
Meng:2017fwb ; Narison:2010py ; Zhang:2008rt ; Guo:2017vcf ; Lu:2017meb ;
Xiao:2017udy ; Weng:2018mmf ; Can:2013zpa ; Branz:2010pq ; Bose:1980vy ;
Patel:2008xs ; SilvestreBrac:1996bg ; Patel:2007gx ; Gadaria:2016omw  of the
DHBs have been examined broadly in literature by the help of the quark models,
potential models, lattice QCD, Feynman-Hellmann theorem, extended
chromomagnetic model, chiral perturbation theory, heavy quark effective
theory, QCD sum rules, perturbative QCD, Faddeev approach, SU(3) flavor
symmetry, nonperturbative string approach, local diquark approach, light-cone
QCD sum rules, light front approach and extended on-mass-shell renormalization
scheme.

Electromagnetic properties one of the major and meaningful parameters of the
DHBs. As the electromagnetic properties symbolise necessary viewpoints of the
intrinsic properties of hadrons, it is quite significant to analyze the baryon
electromagnetic form factors, particularly the MDMs. The magnitude and sign of
the MDM ensure crucial data on size, structure and shape deformations of
baryons. Apparently, determining the MDM is an important step in our
comprehension of the baryon features with regards to quark-gluon degrees of
freedom. In this study, we are going to concentrate on the DHBs (from now on
we will represent these particles as $B^{*}_{QQ}$ ) with spin-parity
$J^{P}=\frac{3}{2}^{+}$, and calculate their MDMs by the help of the light-
cone QCD sum rule (LCSR) approach, which is one of the powerful
nonperturbative methods in hadron physics providing us to calculate properties
of the particles and processes. In LCSR, the hadronic properties are expressed
with regards to the vacuum condensates and the light-cone distribution
amplitudes of the on-shell particles Chernyak:1990ag ; Braun:1988qv ;
Balitsky:1989ry . Since the MDMs are quantities in terms of the properties of
the vacuum and distribution amplitudes of the particles, any uncertainties in
these parameters are reflected to the uncertainties of the predictions of the
MDMs. The first extraction of the MDMs of the spin-$\frac{3}{2}$ doubly heavy
baryons were obtained by Lichtenberg Lichtenberg:1976fi  by means of the naive
quark model. Then, many scientists have used different theoretical models to
compute the spin-$\frac{3}{2}$ doubly heavy baryon MDMs Albertus:2006ya ;
Shah:2016vmd ; Shah:2017liu ; Meng:2017dni ; Dhir:2009ax ; Dhir:2013nka ;
Bernotas:2012nz ; Simonis:2018rld ; Gadaria:2016omw . In Ref. Albertus:2006ya
, they have been calculated the static properties and semileptonic decays of
the DHBs with the help of nonrelativistic quark model. To examination the
dependence of their results on inter-quark interaction they use various quark
potential that contain hyperfine terms coming from one gluon exchange and
Coulomb term as well. In Refs. Shah:2016vmd ; Shah:2017liu , the masses of the
ground, orbital and radial states of the DHBs are evaluated in the framework
of hypercentral constituent quark model with Coulomb plus linear potential.
The MDMs of the ground state DHBs are also extracted. In Ref. Meng:2017dni ,
the MDMs of the spin-$\frac{3}{2}$ doubly charmed baryons are investigated up
to NNLO by means of the heavy baryon chiral perturbation theory. As a by
product, they obtained the numerical values of the doubly bottom and charmed-
bottom baryons up to NLO. In Refs. Dhir:2009ax ; Dhir:2013nka , they are
calculated the MDMs of singly and DHBs with spin-$\frac{3}{2}$ and
spin-$\frac{1}{2}$ in the framework of effective quark mass and screened
charge of quark. In Ref. Bernotas:2012nz , the MDMs of the DHBs are
investigated in the MIT bag model with center of mass motion corrections. In
Ref. Simonis:2018rld , the decay widths, MDMs and M1 transitions of the all
ground states heavy baryons are evaluated by means of the extended MIT bag
model. In Ref. Gadaria:2016omw , they are utilized the extended relativistic
harmonic confinement model to acquire masses and MDMs of the heavy flavored
ground state baryons.

The plan of the manuscript is as follows. In section II, the details of the
MDMs computations for the DHBs with spin-$\frac{3}{2}$ are presented.
Numerical analysis of the LCSR for the MDMs are given in section III. Section
IV is reserved for discussion and concluding remarks.

##  II Formalism

To obtain the MDM of the DHBs by using the LCSR approach, we begin with the
subsequent correlation function,

$\displaystyle\Pi_{\mu\nu\alpha}(p,q)$ $\displaystyle=$ $\displaystyle
i^{2}\int d^{4}x\,\int d^{4}y\,e^{ip\cdot x+iq\cdot y}$ (1)
$\displaystyle\times\langle
0|\mathcal{T}\\{J_{\mu}^{B^{*}_{QQ}}(x)J_{\alpha}(y)J_{\nu}^{B^{*\dagger}_{QQ}}(0)\\}|0\rangle.$

Here, $J_{\mu(\nu)}$ is the interpolating current of the $B^{*}_{QQ}$ baryons
and the electromagnetic current $J_{\alpha}$ is given as,

$J_{\alpha}=\sum_{q=u,d,s,c,b}e_{q}\bar{q}\gamma_{\alpha}q,$ (2)

where $e_{q}$ is the electric charge of the corresponding quark.

From a technical point of view, the Eq. (1) can be reconsidered in a more
appropriate form by the help of external background electromagnetic (EBGEM)
field,

$\displaystyle\Pi_{\mu\nu}(p,q)=i\int d^{4}x\,e^{ip\cdot x}\langle
0|\mathcal{T}\\{J_{\mu}^{B^{*}_{QQ}}(x)J_{\nu}^{B^{*\dagger}_{QQ}}(0)\\}|0\rangle_{F},$
(3)

where F is the EBGEM field and
$F_{\alpha\beta}=i(\varepsilon_{\alpha}q_{\beta}-\varepsilon_{\beta}q_{\alpha})$
with $\varepsilon_{\beta}$ and $q_{\alpha}$ being the polarization and four-
momentum of the EBGEM field, respectively. Since the EBGEM field can be made
arbitrarily small, the correlation function in Eq. (3) can be acquired by
expanding in powers of the EBGEM field,

$\Pi_{\mu\nu}(p,q)=\Pi_{\mu\nu}^{(0)}(p,q)+\Pi_{\mu\nu}^{(1)}(p,q)+....,$ (4)

and keeping only terms $\Pi_{\mu\nu}^{(1)}(p,q)$, which corresponds to the
single photon emission Ball:2002ps . The primary advantage of using the EBGEM
field approach relies on the fact that it separates the soft and hard photon
emissions in an explicitly gauge invariant way Ball:2002ps . The
$\Pi_{\mu\nu}^{(0)}(p,q)$ is the correlation function in the lack of the EBGEM
field, and leads to the two point sum rules of the hadrons, which is not
relevant for our case.

After these general comments, we can now move on deriving the LCSR for the MDM
of the DHBs. The correlation function given in Eq. (3) can be calculated with
regards to hadronic properties, known as hadronic representation. In addition
to this it can be obtained with regards to the quark-gluon properties in the
deep Euclidean region, known as QCD representation. By matching the results of
these representations using the dispersion relation and quark–hadron duality
ansatz, one can acquire the corresponding sum rules.

We begin to evaluate the correlation function with respect to hadronic degrees
of freedom comprising the physical properties of the particles under
investigation. For that purpose, we embed a complete set of $B^{*}_{QQ}$
baryons into the correlation function. So, we get

$\displaystyle\Pi^{Had}_{\mu\nu}(p,q)$ $\displaystyle=\frac{\langle 0\mid
J_{\mu}^{B^{*}_{QQ}}\mid{B^{*}_{QQ}}(p)\rangle}{[p^{2}-m_{{B^{*}_{QQ}}}^{2}]}$
$\displaystyle\langle{B^{*}_{QQ}}(p)\mid{B^{*}_{QQ}}(p+q)\rangle_{F}$
$\displaystyle\frac{\langle{B^{*}_{QQ}}(p+q)\mid\bar{J}_{\nu}^{B^{*}_{QQ}}\mid
0\rangle}{[(p+q)^{2}-m_{{B^{*}_{QQ}}}^{2}]}+...,$ (5)

where the dots stand for contributions of higher states and the continuum. The
matrix elements in Eq. (II) are defined as Pascalutsa:2006up ; Ramalho:2009vc
,

$\displaystyle\langle 0\mid J_{\mu}(0)\mid{B^{*}_{QQ}}(p,s)\rangle$
$\displaystyle=\lambda_{{B^{*}_{QQ}}}u_{\mu}(p,s),$ (6)
$\displaystyle\langle{B^{*}_{QQ}}(p)\mid{B^{*}_{QQ}}(p+q)\rangle_{F}$
$\displaystyle=-e\,\bar{u}_{\mu}(p)\Bigg{\\{}F_{1}(q^{2})g_{\mu\nu}\varepsilon\\!\\!\\!/$
$\displaystyle-\frac{1}{2m_{{B^{*}_{QQ}}}}\Big{[}F_{2}(q^{2})g_{\mu\nu}+F_{4}(q^{2})$
$\displaystyle\times\frac{q_{\mu}q_{\nu}}{(2m_{{B^{*}_{QQ}}})^{2}}\Big{]}\varepsilon\\!\\!\\!/q\\!\\!\\!/+\frac{F_{3}(q^{2})}{(2m_{{B^{*}_{QQ}}})^{2}}$
$\displaystyle\times
q_{\mu}q_{\nu}\varepsilon\\!\\!\\!/\Bigg{\\}}u_{\nu}(p+q),$ (7)

where $\lambda_{{B^{*}_{QQ}}}$ is the residue of ${B^{*}_{QQ}}$ baryon and
$u_{\mu}(p,s)$ is the Rarita-Schwinger spinor. Summation over spins of
${B^{*}_{QQ}}$ baryon is carried out as:

$\displaystyle\sum_{s}u_{\mu}(p,s)\bar{u}_{\nu}(p,s)$ $\displaystyle=$
$\displaystyle-\Big{(}p\\!\\!\\!/+m_{{B^{*}_{QQ}}}\Big{)}\Big{[}g_{\mu\nu}-\frac{1}{3}\gamma_{\mu}\gamma_{\nu}$
(8)
$\displaystyle-\frac{2\,p_{\mu}p_{\nu}}{3\,m^{2}_{{B^{*}_{QQ}}}}+\frac{p_{\mu}\gamma_{\nu}-p_{\nu}\gamma_{\mu}}{3\,m_{{B^{*}_{QQ}}}}\Big{]}.$

Substituting Eqs. (II)-(8) into Eq. (3) for hadronic side we obtain

$\displaystyle\Pi^{Had}_{\mu\nu}(p,q)$
$\displaystyle=-\frac{\lambda_{{}_{{B^{*}_{QQ}}}}^{2}\,\Big{(}p\\!\\!\\!/+m_{{B^{*}_{QQ}}}\Big{)}}{[(p+q)^{2}-m_{{}_{{B^{*}_{QQ}}}}^{2}][p^{2}-m_{{}_{{B^{*}_{QQ}}}}^{2}]}$
$\displaystyle\Big{[}g_{\mu\nu}-\frac{1}{3}\gamma_{\mu}\gamma_{\nu}-\frac{2\,p_{\mu}p_{\nu}}{3\,m^{2}_{{B^{*}_{QQ}}}}+\frac{p_{\mu}\gamma_{\nu}-p_{\nu}\gamma_{\mu}}{3\,m_{{B^{*}_{QQ}}}}\Big{]}$
$\displaystyle\times\Bigg{\\{}F_{1}(q^{2})g_{\mu\nu}\varepsilon\\!\\!\\!/-\frac{1}{2m_{{B^{*}_{QQ}}}}\Big{[}F_{2}(q^{2})g_{\mu\nu}$
$\displaystyle+F_{4}(q^{2})\frac{q_{\mu}q_{\nu}}{(2m_{{B^{*}_{QQ}}})^{2}}\Big{]}\varepsilon\\!\\!\\!/q\\!\\!\\!/+\frac{F_{3}(q^{2})}{(2m_{{B^{*}_{QQ}}})^{2}}q_{\mu}q_{\nu}\varepsilon\\!\\!\\!/\Bigg{\\}}.$
(9)

As a principle, we can achieve the last form of the hadronic representation of
the correlator using the Eq. (II), however we face with two problems. One of
them is related to the fact that not all Lorentz structures appearing in Eq.
(II) are independent. The second problem is the correlator can also get
contributions from spin-1/2 baryons, which should be removed. To eliminate the
spin one half contributions and acquire just independent structures in the
correlator, we perform the ordering for Dirac matrices as
$\gamma_{\mu}p\\!\\!\\!/\varepsilon\\!\\!\\!/q\\!\\!\\!/\gamma_{\nu}$ and
remove terms with $\gamma_{\mu}$ at the beginning, $\gamma_{\nu}$ at the end
and all those proportional to $p_{\mu}$ and $p_{\nu}$ Belyaev:1982cd . As a
result, for hadronic side we get,

$\displaystyle\Pi^{Had}_{\mu\nu}(p,q)=$
$\displaystyle-\frac{\lambda_{{}_{{B^{*}_{QQ}}}}^{2}}{[(p+q)^{2}-m_{{}_{{B^{*}_{QQ}}}}^{2}][p^{2}-m_{{}_{{B^{*}_{QQ}}}}^{2}]}$
$\displaystyle\times\Bigg{[}-g_{\mu\nu}p\\!\\!\\!/\varepsilon\\!\\!\\!/q\\!\\!\\!/\,F_{1}(q^{2})+m_{{B^{*}_{QQ}}}g_{\mu\nu}\varepsilon\\!\\!\\!/q\\!\\!\\!/\,F_{2}(q^{2})$
$\displaystyle+~{}\mbox{other independent structures}\Bigg{]}.$ (10)

The MDM form factor, $G_{M}(q^{2})$, is defined with respect to the form
factors $F_{i}(q^{2})$ in the following manner Pascalutsa:2006up ;
Ramalho:2009vc :

$\displaystyle G_{M}(q^{2})$ $\displaystyle=$
$\displaystyle\big{[}F_{1}(q^{2})+F_{2}(q^{2})\big{]}(1+\frac{4}{5}\tau)-\frac{2}{5}\big{[}F_{3}(q^{2})$
(11) $\displaystyle+F_{4}(q^{2})\big{]}\tau(1+\tau),$

where $\tau=-\frac{q^{2}}{4m^{2}_{{B^{*}_{QQ}}}}$. At $q^{2}=0$, the magnetic
dipole form factors are obtained with respect tothe functions $F_{i}(0)$ as:

$\displaystyle G_{M}(0)$ $\displaystyle=$ $\displaystyle F_{1}(0)+F_{2}(0).$
(12)

The MDM ($\mu_{{B^{*}_{QQ}}}$), is defined in the following way:

$\displaystyle\mu_{{B^{*}_{QQ}}}$ $\displaystyle=$
$\displaystyle\frac{e}{2m_{{B^{*}_{QQ}}}}G_{M}(0).$ (13)

In this study we achieve QCD sum rules for the form factors $F_{i}(q^{2})$ at
first, after that in numerical calculations we will make use of the above
equations to obtain the values of the MDMs using the QCD sum rules for the
form factors. The final form of the hadronic representation with respect to
the chosen structures in momentum space is:

$\displaystyle\Pi^{Had}_{\mu\nu}(p,q)$ $\displaystyle=$
$\displaystyle\Pi_{1}^{Had}g_{\mu\nu}p\\!\\!\\!/\varepsilon\\!\\!\\!/q\\!\\!\\!/\,+\Pi_{2}^{Had}g_{\mu\nu}\varepsilon\\!\\!\\!/q\\!\\!\\!/\,+...,$
(14)

where $\Pi_{i}^{Had}$ are functions of the form factors $F_{i}(q^{2})$ and
other hadronic parameters; and $...$ represents other independent structures.

To obtain the expression of the correlation function with respect to the
quark-gluon parameters, the explicit form for the interpolating current of the
$B^{*}_{QQ}$ baryons needs to be chosen. In this work, we consider the
$B^{*}_{QQ}$ baryons with the quantum numbers $J^{P}=\frac{3}{2}^{+}$. The
interpolating current is given as Aliev:2012iv ,

$\displaystyle J_{\mu}^{B^{*}_{QQ}}(x)$
$\displaystyle={\frac{1}{\sqrt{3}}}\epsilon^{abc}\Big{\\{}(q^{aT}C\gamma_{\mu}Q^{b})Q^{\prime
c}+(q^{aT}C\gamma_{\mu}Q^{\prime b})Q^{c}$
$\displaystyle+(Q^{aT}C\gamma_{\mu}Q^{\prime b})q^{c}\Big{\\}}~{},$ (15)

where $q$ is the light; and $Q$ and $Q^{\prime}$ are the two heavy quarks,
respectively. We give the quark content of the spin-3/2 DHBs in Table 1.

Table 1: The quark content of the spin-3/2 DHBs.
$\begin{array}[]{l|c|c|cc}\hline\cr\hline\cr\mbox{Baryon}&\mbox{ $q$}&\mbox{
$Q$}&\mbox{ $Q^{\prime}$}\\\
\hline\cr\hline\cr\Xi_{QQ}^{\ast}&u~{}\mbox{or}~{}d&b~{}\mbox{or}~{}c&b~{}\mbox{or}~{}c\\\
\Xi_{QQ^{\prime}}^{\ast}&u~{}\mbox{or}~{}d&b&c\\\
\Omega_{QQ}^{\ast}&s&b~{}\mbox{or}~{}c&b~{}\mbox{or}~{}c\\\
\Omega_{QQ^{\prime}}^{\ast}&s&b&c\\\ \hline\cr\hline\cr\end{array}$

After contracting pairs of quark fields and using the Wick’s theorem, the
correlation function becomes:

$\displaystyle\Pi^{QCD}_{\mu\nu}(p)$ $\displaystyle=$
$\displaystyle-\frac{i}{3}\,\varepsilon^{abc}\varepsilon^{a^{\prime}b^{\prime}c^{\prime}}\int
d^{4}xe^{ip\cdot x}\langle
0|\Bigg{\\{}S_{Q^{\prime}}^{cc^{\prime}}\mbox{Tr}\Big{[}S_{Q}^{ba^{\prime}}\gamma_{\nu}\widetilde{S}_{q}^{ab^{\prime}}\gamma_{\mu}\Big{]}+S_{q}^{cc^{\prime}}\mbox{Tr}\Big{[}S_{Q^{\prime}}^{ba^{\prime}}\gamma_{\nu}\widetilde{S}_{Q}^{ab^{\prime}}\gamma_{\mu}\Big{]}$
(16)
$\displaystyle+S_{Q}^{cc^{\prime}}\mbox{Tr}\Big{[}S_{q}^{ba^{\prime}}\gamma_{\nu}\widetilde{S}_{Q^{\prime}}^{ab^{\prime}}\gamma_{\mu}\Big{]}+S_{Q}^{cb^{\prime}}\gamma_{\nu}\widetilde{S}_{Q^{\prime}}^{aa^{\prime}}\gamma_{\mu}S_{q}^{bc^{\prime}}+S_{Q}^{ca^{\prime}}\gamma_{\nu}\widetilde{S}_{q}^{bb^{\prime}}\gamma_{\mu}S_{Q^{\prime}}^{ac^{\prime}}$
$\displaystyle+S_{Q^{\prime}}^{ca^{\prime}}\gamma_{\nu}\widetilde{S}_{Q}^{bb^{\prime}}\gamma_{\mu}S_{q}^{ac^{\prime}}+S_{Q^{\prime}}^{cb^{\prime}}\gamma_{\nu}\widetilde{S}_{q}^{aa^{\prime}}\gamma_{\mu}S_{Q}^{bc^{\prime}}+S_{q}^{ca^{\prime}}\gamma_{\nu}\widetilde{S}_{Q^{\prime}}^{bb^{\prime}}\gamma_{\mu}S_{Q}^{ac^{\prime}}$
$\displaystyle+S_{q}^{cb^{\prime}}\gamma_{\nu}\widetilde{S}_{Q}^{aa^{\prime}}\gamma_{\mu}S_{Q^{\prime}}^{bc^{\prime}}\Bigg{\\}}|0\rangle_{F},$

where $\tilde{S}_{Q(q)}^{ij}(x)=CS_{Q(q)}^{{ij}^{T}}(x)C$ and, $S_{q}^{ij}(x)$
and $S_{Q}^{ij}(x)$ are the light and heavy quark propagators, respectively.
The light and heavy quark propagators are given as Yang:1993bp ;
Belyaev:1985wza ,

$\displaystyle S_{q}(x)$
$\displaystyle=S_{q}^{free}-\frac{\bar{q}q}{12}\Big{(}1-i\frac{m_{q}x\\!\\!\\!/}{4}\Big{)}-\frac{\bar{q}\sigma.Gq}{192}x^{2}\Big{(}1-i\frac{m_{q}x\\!\\!\\!/}{6}\Big{)}$
$\displaystyle-\frac{ig_{s}}{32\pi^{2}x^{2}}~{}G^{\mu\nu}(x)\bigg{[}\hbox
to0.0pt{/\hss}{x}\sigma_{\mu\nu}+\sigma_{\mu\nu}\hbox
to0.0pt{/\hss}{x}\bigg{]},$ $\displaystyle S_{Q}(x)$
$\displaystyle=S_{Q}^{free}-\frac{g_{s}m_{Q}}{16\pi^{2}}\int_{0}^{1}dv~{}G^{\mu\nu}(vx)\bigg{[}\big{(}\sigma_{\mu\nu}{x\\!\\!\\!/}+{x\\!\\!\\!/}\sigma_{\mu\nu}\big{)}$
$\displaystyle\times\frac{K_{1}(m_{Q}\sqrt{-x^{2}})}{\sqrt{-x^{2}}}+2\sigma^{\mu\nu}K_{0}(m_{Q}\sqrt{-x^{2}})\bigg{]},$
(17)

where

$\displaystyle
S_{q}^{free}=\frac{1}{2\pi^{2}x^{2}}\Big{(}i\frac{{x\\!\\!\\!/}}{x^{2}}-\frac{m_{q}}{2}\Big{)},$
$\displaystyle
S_{Q}^{free}=\frac{m_{Q}^{2}}{4\pi^{2}}\bigg{[}\frac{K_{1}(m_{Q}\sqrt{-x^{2}})}{\sqrt{-x^{2}}}+i\frac{{x\\!\\!\\!/}~{}K_{2}(m_{Q}\sqrt{-x^{2}})}{(\sqrt{-x^{2}})^{2}}\bigg{]},$
(18)

with $G^{\mu\nu}$ is the gluon field strength tensor, $K_{i}$ are second kind
of the Bessel functions, $m_{q}$ and $m_{Q}$ are the light and heavy quark
mass respectively.

(a)

(b)

Figure 1:  Feynman diagrams for the MDMs of the spin-3/2 DHBs. The thick,
thin, wavy and curly lines denote the heavy quark, light quark, photon and
gluon propagators, respectively. Diagrams (a) corresponding to the
perturbative photon vertex and, diagrams (b) denote the contributions coming
from the distribution amplitudes of the photon.

The correlator in Eq. (16) contains various contributions: the photon can be
emitted both perturbatively or nonperturbatively. When the photon is emitted
perturbatively, one of the propagators in Eq. (16) is substituted by

$\displaystyle S^{free}(x)\rightarrow\int d^{4}y\,S^{free}(x-y)\,\hbox
to0.0pt{/\hss}{\\!A}(y)\,S^{free}(y)\,,$ (19)

and the surviving two propagators are substituted with the full quark
propagators including the free (perturbative) part as well as the interacting
parts (with gluon or QCD vacuum) as nonperturbative contributions. The total
perturbative photon emission is achieved by carrying out the substitution
mentioned above for the perturbatively interacting quark propagator with the
photon and employing the substitution of the surviving propagators by their
free parts.

In case of nonperturbative photon emission, the light quark propagator in Eq.
(16) is substituted by

$\displaystyle
S_{\alpha\beta}^{ab}\rightarrow-\frac{1}{4}(\bar{q}^{a}\Gamma_{i}q^{b})(\Gamma_{i})_{\alpha\beta},$
(20)

where $\Gamma_{i}$ represent the full set of Dirac matrices. Under this
approach, two surviving quark propagators are taken as the full propagators
comprising perturbative as well as nonperturbative contributions. Once Eq.
(20) is inserted into Eq. (16), there seem matrix elements such as
$\langle\gamma(q)\left|\bar{q}(x)\Gamma_{i}q(0)\right|0\rangle$ and
$\langle\gamma(q)\left|\bar{q}(x)\Gamma_{i}G_{\alpha\beta}q(0)\right|0\rangle$,
representing the nonperturbative contributions. Furthermore, nonlocal
operators such as $\bar{q}q\bar{q}q$ and $\bar{q}G^{2}q$ are anticipated to
appear. In this study, we take into account operators includes only one gluon
field and three particle nonlocal operators and to disregard terms with four
quarks $\bar{q}q\bar{q}q$, and two gluons $\bar{q}G^{2}q$. In order to
calculate the nonperturbative contributions, we need the matrix elements of
the nonlocal operators between the photon states and the vacuum and these
matrix elements are described with respect to the photon distribution
amplitudes with definite twists. Up to twist-4 the explicit expressions of the
photon distribution amplitudes are given in Ball:2002ps . Using these
expressions for the propagators and distribution amplitudes for the photon,
the correlation functions from the QCD side can be computed.

The QCD and hadronic representations of the correlation function are then
matched using dispersion relation. The next step in deriving the sum rules for
the MDMs of the spin-$\frac{3}{2}$ DHBs is applying double Borel
transformations (${\cal B}$) over the $p^{2}$ and $(p+q)^{2}$ on the both
representations of the correlation function so as to suppress the
contributions of higher states and continuum. As a result, we obtain

$\displaystyle{\cal B}\Pi^{Had}_{\mu\nu}(p,q)={\cal
B}\Pi^{QCD}_{\mu\nu}(p,q),$ (21)

which leads to

$\displaystyle{\cal B}\Pi_{1}^{Had}={\cal B}\Pi_{1}^{QCD},~{}~{}~{}~{}{\cal
B}\Pi_{2}^{Had}={\cal B}\Pi_{2}^{QCD},$ (22)

corresponding to the structures
$g_{\mu\nu}p\\!\\!\\!/\varepsilon\\!\\!\\!/q\\!\\!\\!/$ and
$g_{\mu\nu}\varepsilon\\!\\!\\!/q\\!\\!\\!/$. In this manner we extract the
QCD sum rules for the form factors $F_{1}$ and $F_{2}$. They are very lengthy
functions, therefore we do not give their explicit expressions here. The
interested readers can find details of the calculations such as Borel
transformations and continuum subtraction in Refs. Agaev:2016srl ;
Azizi:2018duk

##  III Numerical analysis

In this section, we achieve numerical computations for the spin-$\frac{3}{2}$
DHBs. We use $m_{u}$ = $m_{d}$ = 0, $m_{s}$ = $96^{+8}_{-4}$ MeV, $m_{c}$ =
$1.67\pm 0.07$ GeV, $m_{b}$ = $4.78\pm 0.06$ GeV, $f_{3\gamma}$ =
$-0.0039~{}\text{GeV}^{2}$ Ball:2002ps , $\langle\bar{q}q\rangle$ = $(-0.24\pm
0.01)^{3}\,\text{GeV}^{3}$ Ioffe:2005ym , $m_{0}^{2}$ = $0.8\pm
0.1~{}\text{GeV}^{2}$, $\langle g_{s}^{2}G^{2}\rangle=0.88~{}GeV^{4}$ and
$\chi$ = $-2.85\pm 0.5~{}\text{GeV}^{-2}$ Rohrwild:2007yt . The masses of the
$\Xi^{*}_{QQ}$, $\Xi^{*}_{QQ^{\prime}}$, $\Omega^{*}_{QQ}$ and
$\Omega^{*}_{QQ^{\prime}}$ baryons are borrowed from Ref. Aliev:2012iv , in
which the mass sum rules have been used to compute them. These masses are used
to have the following values: $M_{\Xi^{*}_{cc}}=3.72\pm 0.18~{}\text{GeV}$,
$M_{\Omega^{*}_{cc}}=3.78\pm 0.16~{}\text{GeV}$, $M_{\Xi^{*}_{bc}}=7.25\pm
0.20~{}\text{GeV}$, $M_{\Omega^{*}_{bc}}=7.30\pm 0.20~{}\text{GeV}$,
$M_{\Xi^{*}_{bb}}=10.40\pm 1.00~{}\text{GeV}$ and
$M_{\Omega^{*}_{bb}}=10.50\pm 0.20~{}\text{GeV}$. In order to specify the MDMs
of DHBs, the value of the residues are needed. The residues of the DHBs are
computed in Ref. Aliev:2012iv . These residues are calculated to have the
following values: $\lambda_{\Xi^{*}_{cc}}=0.12\pm 0.01~{}\text{GeV}^{3}$,
$\lambda_{\Omega^{*}_{cc}}=0.14\pm 0.02~{}\text{GeV}^{3}$,
$\lambda_{\Xi^{*}_{bc}}=0.15\pm 0.01~{}\text{GeV}^{3}$,
$\lambda_{\Omega^{*}_{bc}}=0.18\pm 0.02~{}\text{GeV}^{3}$,
$\lambda_{\Xi^{*}_{bb}}=0.22\pm 0.03~{}\text{GeV}^{3}$ and
$\lambda_{\Omega^{*}_{bb}}=0.25\pm 0.03~{}\text{GeV}^{3}$. The parameters used
in the photon distribution amplitudes are given in Ref. Ball:2002ps .

The QCD sum rule for the MDMs of the DHBs, besides the above mentioned input
parameters, include also two more extra parameters. These parameters are the
continuum threshold $s_{0}$ and the Borel mass parameter $M^{2}$. According to
the QCD sum rules philosophy we need to find the working regions of these
parameters, where the MDMs of the DHBs be insensitive to the variation of
these parameters in their working regions. The $s_{0}$ is not entirely
optional parameter, it is preferred as the point at which the continuum and
excited states begin to contribute to the calculations. To decide the working
interval of the $s_{0}$, we enforce the conditions of operator product
expansion (OPE) convergence and pole dominance. Therefore, it is expected that
$s_{0}$ varies in the interval $(M_{B^{*}_{QQ}}+0.3)^{2}$ GeV${}^{2}\leq
s_{0}\leq(M_{B^{*}_{QQ}}+0.7)^{2}$ GeV2. From this point of view, we prefer
the value of the $s_{0}$ within the interval $s_{0}=(16-20)~{}\text{GeV}^{2}$
for $\Xi^{*}_{cc}$, $s_{0}=(58-62)~{}\text{GeV}^{2}$ for $\Xi^{*}_{bc}$,
$s_{0}=(116-120)~{}\text{GeV}^{2}$ for $\Xi^{*}_{bb}$,
$s_{0}=(18-22)~{}\text{GeV}^{2}$ for $\Omega^{*}_{cc}$,
$s_{0}=(60-64)~{}\text{GeV}^{2}$ for $\Omega^{*}_{bc}$ and
$s_{0}=(118-122)~{}\text{GeV}^{2}$ for $\Omega^{*}_{bb}$ baryons. The working
region for $M^{2}$ is achieved by requiring that the series of OPE in QCD
representation is convergent and the contribution of higher states and
continuum is efficiently suppressed. In technique language, the upper limit on
$M^{2}$ is found demanding the maximum pole contribution. The lower limit is
obtained demanding that the contribution of the perturbative part exceeds the
nonperturbative one and series of the operator product expansion in the
obtained sum rules converge. Our numerical calculations indicates that both
conditions are satisfied when $M^{2}$ changes in the regions:
$4~{}\text{GeV}^{2}\leq M^{2}\leq 6~{}\text{GeV}^{2}$ for $\Xi^{*}_{cc}$,
$7~{}\text{GeV}^{2}\leq M^{2}\leq 9~{}\text{GeV}^{2}$ for $\Xi^{*}_{bc}$,
$10~{}\text{GeV}^{2}\leq M^{2}\leq 14~{}\text{GeV}^{2}$ for $\Xi^{*}_{bb}$,
$5~{}\text{GeV}^{2}\leq M^{2}\leq 7~{}\text{GeV}^{2}$ for $\Omega^{*}_{cc}$,
$8~{}\text{GeV}^{2}\leq M^{2}\leq 10~{}\text{GeV}^{2}$ for $\Omega^{*}_{bc}$
and $11~{}\text{GeV}^{2}\leq M^{2}\leq 15~{}\text{GeV}^{2}$ for
$\Omega^{*}_{bb}$ baryons. As an example in Fig. 2, we present the
dependencies of the MDMs of doubly charmed baryons on $M^{2}$ at several fixed
values of the $s_{0}$. As is seen from the figure, although being not entirely
insensitive, the MDMs exhibit acceptable dependency on the extra parameters,
$s_{0}$ and $M^{2}$ which is reasonable in the error limits of the QCD sum
rule formalism.

Our final results on the MDMs for the spin-$\frac{3}{2}$ DHBs are

$\displaystyle\mu_{\Xi_{cc}^{*++}}$ $\displaystyle=$ $\displaystyle 2.94\pm
0.95,$ $\displaystyle\mu_{\Xi_{bc}^{*+}}$ $\displaystyle=$ $\displaystyle
2.63\pm 0.82,$ $\displaystyle\mu_{\Xi_{bb}^{*0}}$ $\displaystyle=$
$\displaystyle 2.30\pm 0.55,$ $\displaystyle\mu_{\Xi_{cc}^{*+}}$
$\displaystyle=$ $\displaystyle-0.67\pm 0.11,$
$\displaystyle\mu_{\Xi_{bc}^{*0}}$ $\displaystyle=$ $\displaystyle-0.96\pm
0.32,$ $\displaystyle\mu_{\Xi_{bb}^{*-}}$ $\displaystyle=$
$\displaystyle-1.39\pm 0.32,$ $\displaystyle\mu_{\Omega_{cc}^{*+}}$
$\displaystyle=$ $\displaystyle-0.52\pm 0.07,$
$\displaystyle\mu_{\Omega_{bc}^{*+}}$ $\displaystyle=$ $\displaystyle-1.11\pm
0.33,$ $\displaystyle\mu_{\Omega_{bb}^{*-}}$ $\displaystyle=$
$\displaystyle-1.56\pm 0.33,$ (23)

where the quoted errors in the results are in connection with the
uncertainties in the values of the input parameters and the photon
distribution amplitudes, as well as the variations in the computations of the
working windows $M^{2}$ and $s_{0}$. We also need to emphasize that the main
source of uncertainties is the variations with respect to $s_{0}$ and the
results weakly depend on the choices of the $M^{2}$.

In Table 2, we present the our numerical results for the MDMs and comparison
of the with various other models such as, the relativistic harmonic
confinement model (RHM) Gadaria:2016omw , MIT bag model Bernotas:2012nz ;
Simonis:2018rld , nonrelativistic quark model (NRQM) Albertus:2006ya , hyper
central constituent model (HCQM) Shah:2016vmd ; Shah:2017liu , effective mass
(EMS) and screened charge scheme (SCS) Dhir:2009ax ; Dhir:2013nka  and heavy
baryon chiral perturbation theory (HBChBT) Meng:2017dni . From a comparison of
our values with the predictions of other models we observe from this table
that for the the $\Xi_{cc}^{*++}$ baryon, practically all methods give,
approximately, similar predictions. For the $\Xi_{cc}^{*+}$ and
$\Omega^{*+}_{cc}$ baryons, there are large discrepancy among results not only
the magnitude but also by the sign. The reason for such inconsistencies is
relatively easy to understand. The sign of the MDM depends on what is stronger
-two heavy quarks and one light quarks. In our analysis, the light quark
overcome two heavy quarks and give the dominant contributions. For the the
$\Xi_{bb}^{*0}$ baryon, our estimation is consistent within the errors with
Refs. Albertus:2006ya ; Gadaria:2016omw ; Meng:2017dni  and unlike other
approaches. For the the $\Xi_{bb}^{*-}$ baryon, nearly all models give,
approximately, similar predictions except the values of Refs. Bernotas:2012nz
; Simonis:2018rld , which are small. For the $\Omega^{*-}_{bb}$ baryon, our
estimation is consistent within the errors with Refs. Shah:2016vmd ;
Shah:2017liu ; Meng:2017dni  and different from other results. For the the
$\Xi_{bc}^{*+}$ baryon, we see that within errors our predictions in good
agreement with the Refs Albertus:2006ya ; Gadaria:2016omw ; Dhir:2009ax ;
Dhir:2013nka ; Meng:2017dni . For the $\Xi^{*0}_{bc}$ baryon, while the sign
of the MDM is correctly determined, there is a large discrepancy among
results. For the $\Omega^{*0}_{bc}$ baryon, nearly all methods give,
moderately, similar approximations except the results of Ref. Meng:2017dni
and this work, which are quite large.

Table 2:  MDMs of the spin-$\frac{3}{2}$ DHBs (in nuclear magnetons $\mu_{N}$). Approaches |  $\Xi_{cc}^{*++}$ |  $\Xi_{cc}^{*+}$ |  $\Omega_{cc}^{*+}$ |  $\Xi_{bb}^{*0}$ |  $\Xi_{bb}^{*-}$ |  $\Omega_{bb}^{*-}$ |  $\Xi_{bc}^{*+}$ |  $\Xi_{bc}^{*0}$ |  $\Omega_{bc}^{*0}$
---|---|---|---|---|---|---|---|---|---
NRQM Albertus:2006ya  | 2.67 | -0.31 | 0.14 | 1.87 | -1.11 | -0.66 | 2.27 | -0.71 | -0.26
HCQM Shah:2016vmd ; Shah:2017liu  | 2.22 | 0.07 | 0.29 | 1.61 | -1.74 | -1.24 | 1.56 | -0.38 | -0.18
HBChBT-I Meng:2017dni  | 3.51 | -0.27 | -0.64 | 2.83 | -1.33 | -1.54 | 3.22 | -0.84 | -1.09
HBChBT-II Meng:2017dni  | 3.63 | -0.37 | -0.65 | 2.87 | -1.38 | -1.55 | 3.27 | -0.89 | -1.10
EMS Dhir:2009ax ; Dhir:2013nka  | 2.41 | -0.11 | 0.16 | 1.60 | -0.98 | -0.70 | 2.01 | -0.55 | -0.28
SCS Dhir:2009ax ; Dhir:2013nka  | 2.52 | 0.04 | 0.21 | 1.50 | -1.02 | -0.80 | 2.02 | -0.50 | -0.30
MIT Bag model-I Bernotas:2012nz  | 2.00 | 0.16 | 0.33 | 0.92 | -0.65 | -0.52 | 1.41 | -0.25 | -0.11
MIT Bag model-II Simonis:2018rld  | 2.35 | -0.18 | -0.05 | 1.40 | -0.88 | -0.70 | 1.88 | -0.54 | -0.33
RHM Gadaria:2016omw  | 2.72 | -0.23 | 0.16 | 2.30 | -1.32 | -0.86 | 2.68 | -0.76 | -0.32
This work | 2.94 | -0.67 | -0.52 | 2.30 | -1.39 | -1.56 | 2.63 | -0.96 | -1.11

There is no experimental data for the MDMs of the DHBs. But, one can define
some useful splittings rules for the checking and comparing the results. In
the heavy quark limit, $\Xi_{QQq}^{*}$, $\Xi_{QQ^{\prime}q}^{*}$,
$\Omega_{QQq}^{*}$ and $\Omega_{QQ^{\prime}q}^{*}$ baryons have the same MDMs.
Because they have the same light degrees of freedom and the contribution of
heavy quarks will vanished in this limit. If one keep the heavy quark mass,
the MDMs are not equal anymore. However one can expect that they may satisfy
following relations as the $\Xi_{QQ^{\prime}q}^{*}$ and
$\Omega_{QQ^{\prime}q}^{*}$ are intermediate states:

$\displaystyle\Delta\mu_{u}\equiv\mu_{\Xi_{bbu}^{*0}}+\mu_{\Xi_{ccu}^{*++}}-2\,\mu_{\Xi_{bcu}^{*+}}=0,$
$\displaystyle\Delta\mu_{d}\equiv\mu_{\Xi_{bbd}^{*-}}+\mu_{\Xi_{ccd}^{*+}}-2\,\mu_{\Xi_{bcd}^{*0}}=0,$
(24)
$\displaystyle\Delta\mu_{s}\equiv\mu_{\Omega_{bbs}^{*-}}+\mu_{\Omega_{ccs}^{*+}}-2\,\mu_{\Omega_{bcs}^{*0}}=0.$

We observe from Table 3 that all baryons almost satisfy that condition except
the results are obtained in Refs. Shah:2016vmd ; Shah:2017liu .

Table 3:  Comparision of the first splitting. Approaches |  $\Delta\mu_{u}$ | $\Delta\mu_{d}$ | $\Delta\mu_{s}$
---|---|---|---
NRQM Albertus:2006ya  | 0.00 | 0.00 | 0.00
HCQM Shah:2016vmd ; Shah:2017liu  | 0.71 | -0.91 | -0.59
HBChBT-I Meng:2017dni  | -0.10 | 0.08 | 0.00
HBChBT-II Meng:2017dni  | -0.04 | 0.03 | 0.00
EMS Dhir:2009ax ; Dhir:2013nka  | -0.01 | 0.01 | 0.02
SCS Dhir:2009ax ; Dhir:2013nka  | -0.02 | 0.02 | 0.01
MIT Bag model-I Bernotas:2012nz  | 0.10 | 0.01 | 0.03
MIT Bag model-II Simonis:2018rld  | -0.01 | 0.02 | -0.09
RHM Gadaria:2016omw  | -0.34 | -0.03 | -0.06
This work | -0.02 | -0.14 | 0.14

As we mentioned beginning of this section, we work on $m_{u}=m_{d}=0$ limit.
The second splittings can be defined

$\displaystyle\Delta{\mu}^{QQ}_{u-d}\equiv\mu_{\Xi^{*}_{QQu}}-\mu_{\Xi^{*}_{QQd}},$
$\displaystyle\Delta{\mu}^{QQ^{\prime}}_{u-d}\equiv\mu_{\Xi^{*}_{QQ^{\prime}u}}-\mu_{\Xi^{*}_{QQd}},$
(25)

which gives the difference occur as a result of the different charge of $u$
and $d$ quark. The third splittings can also be defined,

$\displaystyle\Delta{\mu}^{QQ}_{d-s}\equiv\mu_{\Xi^{*}_{QQd}}-\mu_{\Omega^{*}_{QQs}},$
$\displaystyle\Delta{\mu}^{QQ^{\prime}}_{d-s}\equiv\mu_{\Xi^{*}_{QQ^{\prime}d}}-\mu_{\Omega^{*}_{QQ^{\prime}s}},$
(26)

which gives the difference occur as a result of the different quark mass of
$d$ and $s$. In the second and third splittings, the heavy quark contribution
has been canceled out. Therefore,
$\Delta\mu_{u-d/d-s}^{cc}=\Delta\mu_{u-d/d-s}^{bb}=\Delta\mu_{u-d/d-s}^{bc}$
is expected. We can see most results in the Table 4 satisfy this relation,
except the $\Delta\mu_{d-s}^{cc}$ of this work.

Table 4:  Comparision of the second and third splittings. Approaches |  $\Delta\mu_{u-d}^{cc}$ | $\Delta\mu_{u-d}^{bb}$ | $\Delta\mu_{u-d}^{bc}$ | $\Delta\mu_{d-s}^{cc}$ | $\Delta\mu_{d-s}^{bb}$ | $\Delta\mu_{d-s}^{bc}$
---|---|---|---|---|---|---
NRQM Albertus:2006ya  | 2.98 | 2.98 | 2.98 | -0.45 | -0.45 | -0.45
HCQM Shah:2016vmd ; Shah:2017liu  | 2.15 | 3.35 | 1.94 | -0.22 | -0.50 | -0.20
HBChBT-I Meng:2017dni  | 3.78 | 4.16 | 4.06 | 0.37 | 0.21 | 0.25
HBChBT-II Meng:2017dni  | 4.00 | 4.25 | 4.16 | 0.28 | 0.17 | 0.21
EMS Dhir:2009ax ; Dhir:2013nka  | 2.52 | 2.58 | 2.56 | -0.27 | -0.28 | -0.27
SCS Dhir:2009ax ; Dhir:2013nka  | 2.48 | 2.52 | 2.52 | -0.17 | -0.22 | -0.20
MIT Bag model-I Bernotas:2012nz  | 1.84 | 1.57 | 1.66 | -0.17 | -0.13 | -0.14
MIT Bag model-II Simonis:2018rld  | 2.53 | 2.28 | 2.42 | -0.13 | -0.18 | -0.21
RHM Gadaria:2016omw  | 2.95 | 3.62 | 3.44 | -0.39 | -0.46 | -0.44
This work | 3.61 | 3.69 | 3.59 | -0.15 | 0.17 | 0.15

##  IV Discussion and concluding remarks

In the presented paper we have evaluated the MDMs of the spin-$\frac{3}{2}$
DHBs by means of the light-cone QCD sum rule. The MDMs of the DHBs encodes key
knowledge of their internal structure and shape deformations. Measurement of
the MDMs of the spin-$\frac{3}{2}$ DHBs in future experiments can be very
helpful understanding the internal structure of these baryons. However, the
direct measurement of the MDMs of the spin-$\frac{3}{2}$ DHBs are unlikely in
the near future. Therefore, any unstraightforward approximation of the MDMs of
the spin-$\frac{3}{2}$ DHBs could be very helpful. Comparison of our results
with the estimation of other theoretical models is presented. As can be seen
from the MDM results of the DHBs given in Table 2, the results obtained using
different models lead to rather different estimations, which can be used to
distinguish these models. Obviously, more studies are needed to understand the
current situation.

Figure 2:  The dependence of the MDMs for the spin-$\frac{3}{2}$ doubly
charmed baryons on the $M^{2}$ at various fixed values of the $s_{0}$.

##  V Acknowledgements

We are grateful to V. S. Zamiralov for useful discussions, comments and
remarks.

## References

  * (1) M. Mattson, et al., Phys. Rev. Lett. 89 (2002) 112001.
  * (2) R. Chistov, et al., Phys. Rev. Lett. 97 (2006) 162001.
  * (3) S. P. Ratti, Nucl. Phys. Proc. Suppl. 115 (2003) 33–36.
  * (4) B. Aubert, et al., Phys. Rev. D74 (2006) 011103.
  * (5) R. Aaij, et al., Phys. Rev. Lett. 119 (11) (2017) 112001.
  * (6) R. Aaij, et al. (2019).  arXiv:1911.08594.
  * (7) C. Albertus, E. Hernandez, J. Nieves, J. M. Verde-Velasco, Eur. Phys. J. A32 (2007) 183–199, [Erratum: Eur. Phys. J.A36,119(2008)].
  * (8) R.-H. Li, C.-D. Lü, W. Wang, F.-S. Yu, Z.-T. Zou, Phys. Lett. B767 (2017) 232–235.
  * (9) W. Wang, F.-S. Yu, Z.-X. Zhao, Eur. Phys. J. C77 (11) (2017) 781.
  * (10) W. Wang, Z.-P. Xing, J. Xu, Eur. Phys. J. C77 (11) (2017) 800.
  * (11) Y.-J. Shi, W. Wang, Y. Xing, J. Xu, Eur. Phys. J. C78 (1) (2018) 56.
  * (12) Y.-J. Shi, W. Wang, Z.-X. Zhao (2019).  arXiv:1902.01092.
  * (13) Y.-J. Shi, Y. Xing, Z.-X. Zhao, Eur. Phys. J. C79 (6) (2019) 501.  arXiv:1903.03921.
  * (14) J. Hu, T. Mehen, Phys. Rev. D73 (2006) 054003.
  * (15) L.-Y. Xiao, K.-L. Wang, Q.-f. Lu, X.-H. Zhong, S.-L. Zhu, Phys. Rev. D96 (9) (2017) 094005.
  * (16) H.-S. Li, L. Meng, Z.-W. Liu, S.-L. Zhu, Phys. Lett. B777 (2018) 169–176.
  * (17) F.-S. Yu, H.-Y. Jiang, R.-H. Li, C.-D. Lü, W. Wang, Z.-X. Zhao, Chin. Phys. C42 (5) (2018) 051001.
  * (18) Q.-F. Lü, K.-L. Wang, L.-Y. Xiao, X.-H. Zhong, Phys. Rev. D96 (11) (2017) 114006\.
  * (19) E.-L. Cui, H.-X. Chen, W. Chen, X. Liu, S.-L. Zhu, Phys. Rev. D97 (3) (2018) 034018\.
  * (20) K. U. Can, G. Erkol, B. Isildak, M. Oka, T. T. Takahashi, Phys. Lett. B726 (2013) 703–709.
  * (21) T. Branz, A. Faessler, T. Gutsche, M. A. Ivanov, J. G. Korner, V. E. Lyubovitskij, B. Oexl, Phys. Rev. D81 (2010) 114036.
  * (22) S. K. Bose, L. P. Singh, Phys. Rev. D22 (1980) 773.
  * (23) B. Patel, A. K. Rai, P. C. Vinodkumar, 2008.  arXiv:0803.0221.
  * (24) B. Silvestre-Brac, Few Body Syst. 20 (1996) 1–25.
  * (25) B. Patel, A. K. Rai, P. C. Vinodkumar, J. Phys. G35 (2008) 065001.
  * (26) A. N. Gadaria, N. R. Soni, J. N. Pandya, DAE Symp. Nucl. Phys. 61 (2016) 698–699.
  * (27) B. Julia Diaz, D. O. Riska, Nucl. Phys. A739 (2004) 69 – 88.
  * (28) A. Faessler, T. Gutsche, M. A. Ivanov, J. G. Korner, V. E. Lyubovitskij, D. Nicmorus, K. Pumsa-ard, Phys. Rev. D73 (2006) 094013.
  * (29) K. U. Can, G. Erkol, B. Isildak, M. Oka, T. T. Takahashi, JHEP 05 (2014) 125.
  * (30) H.-S. Li, L. Meng, Z.-W. Liu, S.-L. Zhu, Phys. Rev. D96 (7) (2017) 076011.
  * (31) A. Bernotas, V. Simonis (2012).  arXiv:1209.2900.
  * (32) D. B. Lichtenberg, Phys. Rev. D15 (1977) 345.
  * (33) Y.-s. Oh, D.-P. Min, M. Rho, N. N. Scoccola, Nucl. Phys. A534 (1991) 493–512.
  * (34) V. Simonis (2018).  arXiv:1803.01809.
  * (35) M.-Z. Liu, Y. Xiao, L.-S. Geng, Phys. Rev. D98 (1) (2018) 014040.
  * (36) A. N. Hiller Blin, Z.-F. Sun, M. J. Vicente Vacas, Phys. Rev. D98 (5) (2018) 054025.
  * (37) L. Meng, H.-S. Li, Z.-W. Liu, S.-L. Zhu, Eur. Phys. J. C77 (12) (2017) 869.
  * (38) R. Dhir, R. C. Verma, Eur. Phys. J. A42 (2009) 243–249.
  * (39) U. Özdem, J. Phys. G46 (2019) 035003.
  * (40) E. Bagan, M. Chabab, S. Narison, Phys. Lett. B306 (1993) 350–356.
  * (41) R. Roncaglia, D. B. Lichtenberg, E. Predazzi, Phys. Rev. D52 (1995) 1722–1725.
  * (42) D. Ebert, R. N. Faustov, V. O. Galkin, A. P. Martynenko, V. A. Saleev, Z. Phys. C76 (1997) 111–115.
  * (43) S.-P. Tong, Y.-B. Ding, X.-H. Guo, H.-Y. Jin, X.-Q. Li, P.-N. Shen, R. Zhang, Phys. Rev. D62 (2000) 054024.
  * (44) C. Itoh, T. Minamikawa, K. Miura, T. Watanabe, Phys. Rev. D61 (2000) 057502.
  * (45) S. S. Gershtein, V. V. Kiselev, A. K. Likhoded, A. I. Onishchenko, Phys. Rev. D62 (2000) 054021.
  * (46) V. V. Kiselev, A. K. Likhoded, Phys. Usp. 45 (2002) 455–506.
  * (47) V. V. Kiselev, A. K. Likhoded, O. N. Pakhomova, V. A. Saleev, Phys. Rev. D66 (2002) 034030.
  * (48) I. M. Narodetskii, M. A. Trusov, Phys. Atom. Nucl. 65 (2002) 917–924, [Yad. Fiz.65,949(2002)].
  * (49) R. Lewis, N. Mathur, R. M. Woloshyn, Phys. Rev. D64 (2001) 094509.
  * (50) D. Ebert, R. N. Faustov, V. O. Galkin, A. P. Martynenko, Phys. Rev. D66 (2002) 014008\.
  * (51) N. Mathur, R. Lewis, R. M. Woloshyn, Phys. Rev. D66 (2002) 014502.
  * (52) J. M. Flynn, F. Mescia, A. S. B. Tariq, JHEP 07 (2003) 066.
  * (53) J. Vijande, H. Garcilazo, A. Valcarce, F. Fernandez, Phys. Rev. D70 (2004) 054022\.
  * (54) T.-W. Chiu, T.-H. Hsieh, Nucl. Phys. A755 (2005) 471–474.
  * (55) S. Migura, D. Merten, B. Metsch, H.-R. Petry, Eur. Phys. J. A28 (2006) 41.
  * (56) A. P. Martynenko, Phys. Lett. B663 (2008) 317–321.
  * (57) L. Tang, X.-H. Yuan, C.-F. Qiao, X.-Q. Li, Commun. Theor. Phys. 57 (2012) 435–444.
  * (58) X. Liu, H.-X. Chen, Y.-R. Liu, A. Hosaka, S.-L. Zhu, Phys. Rev. D77 (2008) 014031\.
  * (59) W. Roberts, M. Pervin, Int. J. Mod. Phys. A23 (2008) 2817–2860.
  * (60) A. Valcarce, H. Garcilazo, J. Vijande, Eur. Phys. J. A37 (2008) 217–225.
  * (61) L. Liu, H.-W. Lin, K. Orginos, A. Walker-Loud, Phys. Rev. D81 (2010) 094505.
  * (62) C. Alexandrou, J. Carbonell, D. Christaras, V. Drach, M. Gravina, M. Papinutto, Phys. Rev. D86 (2012) 114501.
  * (63) T. M. Aliev, K. Azizi, M. Savci, Nucl. Phys. A895 (2012) 59–70.
  * (64) T. M. Aliev, K. Azizi, M. Savci, J. Phys. G40 (2013) 065003.
  * (65) Y. Namekawa, et al., Phys. Rev. D87 (9) (2013) 094512.
  * (66) M. Karliner, J. L. Rosner, Phys. Rev. D90 (9) (2014) 094007.
  * (67) Z.-F. Sun, Z.-W. Liu, X. Liu, S.-L. Zhu, Phys. Rev. D91 (9) (2015) 094030.
  * (68) H.-X. Chen, W. Chen, Q. Mao, A. Hosaka, X. Liu, S.-L. Zhu, Phys. Rev. D91 (5) (2015) 054034.
  * (69) Z.-F. Sun, M. J. Vicente Vacas, Phys. Rev. (9) (2016) 094002.
  * (70) Z. Shah, K. Thakkar, A. K. Rai, Eur. Phys. J. C76 (10) (2016) 530.
  * (71) V. V. Kiselev, A. V. Berezhnoy, A. K. Likhoded, Phys. Atom. Nucl. 81 (3) (2018) 369–372, [Yad. Fiz.81, no.3, 356(2018)].
  * (72) H.-X. Chen, Q. Mao, W. Chen, X. Liu, S.-L. Zhu, Phys. Rev. D96 (3) (2017) 031501, [Erratum: Phys. Rev.D96,no.11,119902(2017)].
  * (73) L. Meng, N. Li, S.-L. Zhu, Phys. Rev. D95 (11) (2017) 114019.
  * (74) S. Narison, R. Albuquerque, Phys. Lett. B694 (2011) 217–225.
  * (75) J.-R. Zhang, M.-Q. Huang, Phys. Rev. D78 (2008) 094007.
  * (76) Z.-H. Guo, Phys. Rev. D96 (7) (2017) 074004.
  * (77) X.-Z. Weng, X.-L. Chen, W.-Z. Deng, Phys. Rev. D97 (5) (2018) 054008.
  * (78) V. L. Chernyak, I. R. Zhitnitsky, Nucl. Phys. B345 (1990) 137–172.
  * (79) V. M. Braun, I. E. Filyanov, Z. Phys. C44 (1989) 157, [Yad. Fiz.50,818(1989)].
  * (80) I. I. Balitsky, V. M. Braun, A. V. Kolesnichenko, Nucl. Phys. B312 (1989) 509–550.
  * (81) Z. Shah, A. K. Rai, Eur. Phys. J. C77 (2) (2017) 129.
  * (82) R. Dhir, C. S. Kim, R. C. Verma, Phys. Rev. D88 (2013) 094002.
  * (83) P. Ball, V. M. Braun, N. Kivel, Nucl. Phys. B649 (2003) 263–296.
  * (84) V. Pascalutsa, M. Vanderhaeghen, S. N. Yang, Phys. Rept. 437 (2007) 125–232.
  * (85) G. Ramalho, M. T. Pena, F. Gross, Phys. Lett. B678 (2009) 355–358.
  * (86) V. M. Belyaev, B. L. Ioffe, Sov. Phys. JETP 57 (1983) 716–721, [Zh. Eksp. Teor. Fiz.84,1236(1983)].
  * (87) K.-C. Yang, W. Y. P. Hwang, E. M. Henley, L. S. Kisslinger, Phys. Rev. D47 (1993) 3001–3012.
  * (88) V. M. Belyaev, B. Yu. Blok, Z. Phys. C30 (1986) 151.
  * (89) S. S. Agaev, K. Azizi, H. Sundu, Phys. Rev. D93 (11) (2016) 114036.
  * (90) K. Azizi, A. R. Olamaei, S. Rostami, Eur. Phys. J. A54 (9) (2018) 162.
  * (91) B. L. Ioffe, Prog. Part. Nucl. Phys. 56 (2006) 232–277.
  * (92) J. Rohrwild, JHEP 09 (2007) 073.